diff options
| -rwxr-xr-x | Abfind.f | 28 | ||||
| -rwxr-xr-x | Abpop.f | 26 | ||||
| -rwxr-xr-x | Abunplot.f | 20 | ||||
| -rwxr-xr-x | Begin.f | 22 | ||||
| -rwxr-xr-x | Binary.f | 28 | ||||
| -rwxr-xr-x | Binplot.f | 64 | ||||
| -rwxr-xr-x | Binplotprep.f | 18 | ||||
| -rwxr-xr-x | Blends.f | 66 | ||||
| -rwxr-xr-x | Boxit.f | 10 | ||||
| -rwxr-xr-x | Cdcalc.f | 8 | ||||
| -rwxr-xr-x | Chabund.f | 62 | ||||
| -rwxr-xr-x | Cog.f | 8 | ||||
| -rwxr-xr-x | Cogplot.f | 14 | ||||
| -rwxr-xr-x | Cogsyn.f | 6 | ||||
| -rwxr-xr-x | Correl.f | 14 | ||||
| -rwxr-xr-x | Crosscorr.f | 2 | ||||
| -rwxr-xr-x | Curve.f | 10 | ||||
| -rwxr-xr-x | Damping.f | 62 | ||||
| -rwxr-xr-x | Defcolor.f | 4 | ||||
| -rwxr-xr-x | Discov.f | 2 | ||||
| -rwxr-xr-x | Doflux.f | 10 | ||||
| -rwxr-xr-x | Drawcurs.f | 10 | ||||
| -rwxr-xr-x | Eqlib.f | 44 | ||||
| -rwxr-xr-x | Ewfind.f | 12 | ||||
| -rwxr-xr-x | Ewweighted.f | 8 | ||||
| -rwxr-xr-x | Fakeline.f | 6 | ||||
| -rwxr-xr-x | Findtic.f | 8 | ||||
| -rwxr-xr-x | Finish.f | 34 | ||||
| -rwxr-xr-x | Gammabark.f | 22 | ||||
| -rwxr-xr-x | Getasci.f | 4 | ||||
| -rwxr-xr-x | Getcount.f | 2 | ||||
| -rwxr-xr-x | Getnum.f | 4 | ||||
| -rwxr-xr-x | Getsyns.f | 18 | ||||
| -rwxr-xr-x | Gridplo.f | 14 | ||||
| -rwxr-xr-x | Gridsyn.f | 18 | ||||
| -rwxr-xr-x | Infile.f | 24 | ||||
| -rwxr-xr-x | Inlines.f | 78 | ||||
| -rwxr-xr-x | Inmodel.f | 76 | ||||
| -rwxr-xr-x | Invert.f | 10 | ||||
| -rwxr-xr-x | Jexpint.f | 4 | ||||
| -rwxr-xr-x | Lineabund.f | 12 | ||||
| -rwxr-xr-x | Lineinfo.f | 78 | ||||
| -rwxr-xr-x | Linlimit.f | 48 | ||||
| -rwxr-xr-x | Makeplot.f | 44 | ||||
| -rwxr-xr-x | Molquery.f | 18 | ||||
| -rwxr-xr-x | Moog.f | 34 | ||||
| -rwxr-xr-x | Moogsilent.f | 36 | ||||
| -rwxr-xr-x | Nansi.f | 20 | ||||
| -rwxr-xr-x | Nearly.f | 54 | ||||
| -rwxr-xr-x | Number.f | 4 | ||||
| -rwxr-xr-x | Obshead.f | 54 | ||||
| -rwxr-xr-x | Oneline.f | 44 | ||||
| -rwxr-xr-x | OpacHelium.f | 2 | ||||
| -rwxr-xr-x | OpacHydrogen.f | 34 | ||||
| -rwxr-xr-x | Opaccouls.f | 8 | ||||
| -rwxr-xr-x | Opacit.f | 12 | ||||
| -rwxr-xr-x | Opacmetals.f | 66 | ||||
| -rwxr-xr-x | Opacscat.f | 2 | ||||
| -rwxr-xr-x | Params.f | 176 | ||||
| -rwxr-xr-x | Partfn.f | 2 | ||||
| -rwxr-xr-x | Partnew.f | 2 | ||||
| -rwxr-xr-x | Plotit.f | 12 | ||||
| -rwxr-xr-x | Plotremember.f | 12 | ||||
| -rwxr-xr-x | Pltabun.f | 16 | ||||
| -rwxr-xr-x | Pltcog.f | 16 | ||||
| -rwxr-xr-x | Pltflux.f | 10 | ||||
| -rwxr-xr-x | Pltspec.f | 116 | ||||
| -rwxr-xr-x | Pointcurs.f | 4 | ||||
| -rwxr-xr-x | Prinfo.f | 10 | ||||
| -rwxr-xr-x | Readobs.f | 22 | ||||
| -rwxr-xr-x | Rinteg.f | 4 | ||||
| -rwxr-xr-x | Setmols.f | 12 | ||||
| -rwxr-xr-x | Smooth.f | 152 | ||||
| -rwxr-xr-x | Specplot.f | 82 | ||||
| -rwxr-xr-x | Stats.f | 24 | ||||
| -rwxr-xr-x | Sunder.f | 2 | ||||
| -rwxr-xr-x | Synpop.f | 34 | ||||
| -rwxr-xr-x | Synspec.f | 40 | ||||
| -rwxr-xr-x | Synth.f | 24 | ||||
| -rwxr-xr-x | Tablepop.f | 40 | ||||
| -rwxr-xr-x | Taukap.f | 2 | ||||
| -rwxr-xr-x | Total.f | 4 | ||||
| -rwxr-xr-x | Trudamp.f | 30 | ||||
| -rwxr-xr-x | Ucalc.f | 14 | ||||
| -rwxr-xr-x | Vargauss.f | 40 | ||||
| -rwxr-xr-x | Vmacro.f | 2 | ||||
| -rwxr-xr-x | Voigt.f | 14 | ||||
| -rwxr-xr-x | Wavecalc.f | 8 | ||||
| -rwxr-xr-x | Weedout.f | 26 | ||||
| -rwxr-xr-x | Writenumber.f | 4 |
90 files changed, 1203 insertions, 1203 deletions
@@ -72,7 +72,7 @@ c*****set some parameters c*****find the range of lines of a species 5 call linlimit - if (lim1line .lt. 0) then + if (lim1line < 0) then call finish (0) return endif @@ -86,7 +86,7 @@ c*****find out whether molecular equilibrium is involved in the species c*****force each abundance of a species member to predict the c line equivalent width; here is the code for ordinary species - if (molflag .eq. 0) then + if (molflag == 0) then abundin = dlog10(xabund(iabatom)) + 12.0 do lim1=lim1line,lim2line call lineabund (abundin) @@ -106,12 +106,12 @@ c agreement enddo call stats call lineinfo (3) - if (t(jtau5).lt.3800 .or. - . atom1(lim1line).gt.100.0 .or. - . int(atom1(lim1line)+0.0001).eq.6 .or. - . int(atom1(lim1line)+0.0001).eq.8) then - if (iternumber .lt. 6) then - if (dabs(average-abundin) .gt. 0.02) then + if (t(jtau5)<3800 .or. + . atom1(lim1line)>100.0 .or. + . int(atom1(lim1line)+0.0001)==6 .or. + . int(atom1(lim1line)+0.0001)==8) then + if (iternumber < 6) then + if (dabs(average-abundin) > 0.02) then xabund(iabatom) = 10.**(average-12.) iternumber = iternumber + 1 call eqlib @@ -137,15 +137,15 @@ c agreement c*****here a plot may be made on the terminal (and paper) if there c are enough lines; then the user will be prompted on some c options concerning what is seen on the plot - if (plotopt .ne. 0) then + if (plotopt /= 0) then call pltabun - if (choice.eq.'v') then + if (choice=='v') then rewind nf1out rewind nf2out write (nf2out,1002) linitle,moditle choice = ' ' go to 100 - elseif (choice .eq. 'm') then + elseif (choice == 'm') then close (unit=nfmodel) close (unit=nflines) rewind nf1out @@ -166,7 +166,7 @@ c options concerning what is seen on the plot c*****quit, or go on to another species? - if (silent .eq. 'y') then + if (silent == 'y') then choice = 'y' nchars = 0 else @@ -175,8 +175,8 @@ c*****quit, or go on to another species? call getasci (nchars,maxline) choice = chinfo(1:1) endif - if (choice.eq.'y' .or. nchars.le.0) then - if (mode .eq. 2) then + if (choice=='y' .or. nchars<=0) then + if (mode == 2) then go to 5 else call finish (0) @@ -40,7 +40,7 @@ c*****read in the model atmospheres and their summary output files line = synpre num = 80 call getcount (num,line) - if (mmod .lt. 10) then + if (mmod < 10) then write(line(num+1:num+1),1013) mmod else write(line(num+1:num+2),1014) mmod @@ -57,7 +57,7 @@ c*****read in the model atmospheres and their summary output files line = modpre num = 80 call getcount (num,line) - if (mmod .lt. 10) then + if (mmod < 10) then write(line(num+1:num+1),1013) mmod else write(line(num+1:num+2),1014) mmod @@ -83,7 +83,7 @@ c into 1-line lists for computational efficiency call infile ('input ',nflines,'formatted ',0,nchars, . flines,lscreen) call inlines (1) - if (nlines .gt. 1000) then + if (nlines > 1000) then write (*,1005) stop endif @@ -99,7 +99,7 @@ c*****do the curves of growth; store the results lim2 = lim1 waveold = 0. call curve - if (ncurve .gt. 50) then + if (ncurve > 50) then nmodcurve(mmod,lim1) = 50 else nmodcurve(mmod,lim1) = ncurve @@ -128,7 +128,7 @@ c*****set some parameters c*****now, inside the c*****define the range of lines for a species 5 call linlimit - if (lim1line .lt. 0) then + if (lim1line < 0) then call finish (0) return endif @@ -171,7 +171,7 @@ c compute a weighted <EW_calc> call ewweighted rwlgcal = dlog10(ewmod(lim1)/wave1(lim1)) do i=2,ntabtot - if (rwtab(i) .gt. rwlgcal) then + if (rwtab(i) > rwlgcal) then gflgcal = gftab(i-1) + (gftab(i)-gftab(i-1))* . (rwlgcal-rwtab(i-1))/(rwtab(i)-rwtab(i-1)) exit @@ -179,7 +179,7 @@ c compute a weighted <EW_calc> enddo rwlgobs = dlog10(width(lim1)/wave1(lim1)) do i=2,ntabtot - if (rwtab(i) .gt. rwlgobs) then + if (rwtab(i) > rwlgobs) then gflgobs = gftab(i-1) + (gftab(i)-gftab(i-1))* . (rwlgobs-rwtab(i-1))/(rwtab(i)-rwtab(i-1)) exit @@ -188,8 +188,8 @@ c compute a weighted <EW_calc> rwlgerror = rwlgobs - rwlgcal diffngf = gflgobs - gflgcal deltangf = deltangf + diffngf - if (dabs(rwlgerror) .lt. 0.01) exit - if (k .eq. 30) then + if (dabs(rwlgerror) < 0.01) exit + if (k == 30) then write (*,1004) stop endif @@ -210,7 +210,7 @@ c*****here a plot may be made on the terminal (and paper) if there c are enough lines; then the user will be prompted on some c options concerning what is seen on the plot - if (plotopt .ne. 0) then + if (plotopt /= 0) then call blankstring (moditle) moditle(1:70) = popitle(1:70) moditle(57:80) = 'EW-POPULATION ' @@ -220,15 +220,15 @@ c options concerning what is seen on the plot *****quit, or go on to another species? array = 'DO ANOTHER SPECIES ([y]/n)? ' - if (silent .eq. 'n') then + if (silent == 'n') then nchars = 28 call getasci (nchars,maxline) choice = chinfo(1:1) else choice = 'y' endif - if (choice.eq.'y' .or. nchars.le.0) then - if (mode .eq. 2) then + if (choice=='y' .or. nchars<=0) then + if (mode == 2) then go to 5 else call finish (0) @@ -19,7 +19,7 @@ c****************************************************************************** c*****dump the data into working arrays j = 0 do l=lim1obs,lim2obs - if (abundout(l) .ne. 999.99) then + if (abundout(l) /= 999.99) then j = j + 1 ep(j) = e(l,1) abb(j) = abundout(l) @@ -39,7 +39,7 @@ c*****find the plot boundaries for the excitation potential plot do j=1,kount xhi = amax1(xhi,ep(j)) enddo - if (xhi-xlo .lt. 5.) then + if (xhi-xlo < 5.) then xlo = amax1((xlo+xhi)/2.-2.5,-0.2) xhi = xlo + 5.0 else @@ -54,10 +54,10 @@ c*****find the plot boundaries for the excitation potential plot do j=1,kount yhi = amax1(yhi,abb(j)) enddo - if (yhi-ylo .lt. 0.5) then + if (yhi-ylo < 0.5) then ylo = (ylo+yhi)/2. - 0.30 yhi = ylo + 0.60 - elseif (yhi-ylo .lt. 1.0) then + elseif (yhi-ylo < 1.0) then ylo = (ylo+yhi)/2. - 0.55 yhi = ylo + 1.10 else @@ -104,7 +104,7 @@ c*****make the excitation potential plot call sm_ltype (2) call sm_relocate (xlo,ymed) call sm_draw (xhi,ymed) - if (kount .gt. 2 .and. deltaep .gt. 1.5) then + if (kount > 2 .and. deltaep > 1.5) then call sm_ltype (3) call defcolor (3) call sm_relocate (xlo,real(xxm1*xlo+xxb1)) @@ -116,11 +116,11 @@ c*****make the excitation potential plot call sm_relocate (xlo+0.05*(xhi-xlo),ylo+0.15*(yhi-ylo)) call sm_expand (0.8) ich = idint(charge(lim1obs) + 0.1) - if (ich .eq. 1) then + if (ich == 1) then ion = ' I ' - elseif (ich .eq. 2) then + elseif (ich == 2) then ion = ' II ' - elseif (ich .eq. 3) then + elseif (ich == 3) then ion = ' III' endif iatom = idint(atom1(lim1obs)) @@ -178,7 +178,7 @@ c*****make the equivalent width plot call sm_ltype (2) call sm_relocate (xlo,ymed) call sm_draw (xhi,ymed) - if (kount .gt. 2 .and. deltarw .gt. 0.5) then + if (kount > 2 .and. deltarw > 0.5) then call sm_ltype (3) call defcolor (3) call sm_relocate (xlo,real(xxm2*xlo+xxb2)) @@ -242,7 +242,7 @@ c*****make the wavelength plot, and exit normally call sm_lweight (4.0) call sm_ltype (2) call sm_draw (xhi,ymed) - if (kount .gt. 2 .and. deltawv .gt. 500.) then + if (kount > 2 .and. deltawv > 500.) then call sm_ltype (3) call defcolor (3) call sm_relocate (xlo,real(xxm3*xlo+xxb3)) @@ -16,7 +16,7 @@ c*************************************************************************** c*****define the number of text screen lines for silent mode; c this number is hardwired, since it is not really needed at run time. - if (silent .eq. 'y') then + if (silent == 'y') then maxline = 24 write (*,*) 'maxline', maxline pause @@ -42,12 +42,12 @@ c system. open (99,file='/tmp/moog.tmpsize') 5 read (99,1010,end=15) line do i=1,77 - if (line(i:i+3) .eq. 'rows') then - if (machine .eq. 'Linux') then + if (line(i:i+3) == 'rows') then + if (machine == 'Linux') then read (line(i+4:i+6),1011) maxline - elseif (machine .eq. 'Darwin') then + elseif (machine == 'Darwin') then read (line(i-4:i-2),1011) maxline - elseif (machine .eq. 'Solaris') then + elseif (machine == 'Solaris') then read (line(i+6:i+8),1011) maxline endif go to 10 @@ -59,7 +59,7 @@ c system. ikount = 2 call getasci (nchars,ikount) choice = chinfo(1:1) - if (choice.eq.'y' .or. nchars.le.0) then + if (choice=='y' .or. nchars<=0) then go to 10 else call finish (0) @@ -67,7 +67,7 @@ c system. 10 close (99,status='delete') write (systemcall,*) 'rm -f /tmp/moog.tmpsize' call system (systemcall) - if (maxline .lt. 10) then + if (maxline < 10) then maxline = 24 else maxline = maxline - 2 @@ -83,7 +83,7 @@ c*****open data files carried with the source code: Barklem damping nfbarklem = 35 num = 60 call getcount (num,moogpath) - if (moogpath(num:num) .ne. '/') then + if (moogpath(num:num) /= '/') then num = num + 1 moogpath(num:num) = '/' endif @@ -96,7 +96,7 @@ c*****open data files carried with the source code: Barklem UV damping nfbarklemUV = 36 num = 60 call getcount (num,moogpath) - if (moogpath(num:num) .ne. '/') then + if (moogpath(num:num) /= '/') then num = num + 1 moogpath(num:num) = '/' endif @@ -116,10 +116,10 @@ c write a header and find the appropriate parameter file, and exit normally lscreen = 4 nargs = command_argument_count() - if (nargs .gt. 0) then + if (nargs > 0) then call get_command_argument(1, fparam) else - if (silent .eq. 'y') then + if (silent == 'y') then fparam = 'batch.par' else fparam = 'no_filename_given' @@ -22,8 +22,8 @@ c*****examine the parameter file ncall = 1 1 call params linprintopt = linprintalt - if (begin .eq. 0) then - if (numpecatom .gt. 0) then + if (begin == 0) then + if (numpecatom > 0) then do i=3,95 binpec(syncount,i) = pec(i) do j=1,numatomsyn @@ -32,7 +32,7 @@ c*****examine the parameter file enddo endif else - if (numpecatom .gt. 0) then + if (numpecatom > 0) then do i=3,95 pec(i) = binpec(syncount,i) do j=1,numatomsyn @@ -58,12 +58,12 @@ c spectra, and (if desired) IRAF-style smoothed spectra nchars = 20 call infile ('output ',nf2out,'formatted ',0,nchars, . f2out,lscreen) - if (syncount .eq. 1) then + if (syncount == 1) then f7out = f2out else f8out = f2out endif - if (iraf .ne. 0) then + if (iraf /= 0) then nf4out = 23 lscreen = lscreen + 2 array = 'IRAF ("rtext") OUTPUT' @@ -90,7 +90,7 @@ c*****open the line list file and the strong line list file nchars = 13 call infile ('input ',nflines,'formatted ',0,nchars, . flines,lscreen) - if (dostrong .gt. 0) then + if (dostrong > 0) then nfslines = 32 lscreen = lscreen + 2 array = 'THE STRONG LINE LIST' @@ -101,7 +101,7 @@ c*****open the line list file and the strong line list file c*****do the syntheses - if (numpecatom .eq. 0 .or. numatomsyn .eq. 0) then + if (numpecatom == 0 .or. numatomsyn == 0) then isorun = 1 nlines = 0 mode = 3 @@ -123,7 +123,7 @@ c*****do the syntheses linprintopt = 0 enddo endif - if (syncount .eq. 1) then + if (syncount == 1) then fluxprimary = flux else fluxsecondary = flux @@ -133,7 +133,7 @@ c*****do the syntheses c*****finish the syntheses call finish (1) istat = ivcleof(4,1) - if (control .ne. 'gridend') go to 1 + if (control /= 'gridend') go to 1 c*****combine the synthetic spectra for plotting @@ -141,12 +141,12 @@ c*****combine the synthetic spectra for plotting c*****now plot the spectrum, maybe iterating abundances, and end the program - if (plotopt.eq.2 .and. specfileopt.gt.0) then + if (plotopt==2 .and. specfileopt>0) then nfobs = 33 lscreen = lscreen + 2 array = 'THE OBSERVED SPECTRUM' nchars = 21 - if (specfileopt.eq.1 .or. specfileopt.eq.3) then + if (specfileopt==1 .or. specfileopt==3) then call infile ('input ',nfobs,'unformatted',2880,nchars, . fobs,lscreen) else @@ -154,14 +154,14 @@ c*****now plot the spectrum, maybe iterating abundances, and end the program . fobs,lscreen) endif endif - if (plotopt .ne. 0) then + if (plotopt /= 0) then control = 'binary ' nf2out = nf9out nf3out = nf10out call pltspec (lscreen,ncall) - if (choice .eq. 'n') then + if (choice == 'n') then do syncount=1,2 - if (numpecatom .gt. 0) then + if (numpecatom > 0) then do i=3,95 pec(i) = binpec(syncount,i) do j=1,numatomsyn @@ -18,7 +18,7 @@ c****************************************************************************** c*****for grid syntheses, dump out relevant information to a file - if (choice .eq. 'g') then + if (choice == 'g') then write (nf6out,3001) syncount write (nf6out,3002) obsitle, moditle, smitle endif @@ -43,7 +43,7 @@ c*****write smoothing information at the top of the plot c*****write the non-varying isotopic information at the top of the plot do i=1,120 - if (isoitle(i:i) .ne. ' ') go to 112 + if (isoitle(i:i) /= ' ') go to 112 enddo isoitle(1:16) = 'no isotopic data' 112 call sm_relocate (-0.080,1.065) @@ -51,7 +51,7 @@ c*****write the non-varying isotopic information at the top of the plot c*****define the real plot limits - if (xlo .lt. xhi) then + if (xlo < xhi) then call sm_limits (xlo,xhi,ylo,yhi) iflip = 0 else @@ -65,7 +65,7 @@ c*****define the real plot limits c*****draw and label the box for the spectra call defcolor (1) - if (whichwin .eq. '1of1') then + if (whichwin == '1of1') then idev = 1 call sm_window (1,1,1,1,1,1) else @@ -79,7 +79,7 @@ c*****draw and label the box for the spectra call sm_lweight (2.0) call sm_expand (0.8) call sm_box (1,2,4,4) - if (iflip .eq. 1) then + if (iflip == 1) then array = 'Wavenumber' else array = 'Wavelength' @@ -93,15 +93,15 @@ c*****plot the synthetic spectra call sm_lweight (2.2) call sm_expand (0.7) do i=1,100 - if (pec(i) .ne. 0) go to 111 + if (pec(i) /= 0) go to 111 enddo 111 do j=1,nsyn - if (choice.eq.'h' .or. choice.eq.'f' .or. - . choice.eq.'g') then + if (choice=='h' .or. choice=='f' .or. + . choice=='g') then call defcolor (8) call sm_ltype (j-1) else - if (smterm(1:3) .eq. 'x11') then + if (smterm(1:3) == 'x11') then call defcolor (j+1) call sm_ltype (0) else @@ -110,7 +110,7 @@ c*****plot the synthetic spectra endif endif call sm_connect (xsyn,chunk(1,j),kount) - if (iflip .eq. 1) then + if (iflip == 1) then call sm_relocate (xhi+0.045*(xlo-xhi), . ylo+(0.12+0.06*j)*(yhi-ylo)) call sm_draw (xhi+0.005*(xlo-xhi), @@ -134,9 +134,9 @@ c*****plot the synthetic spectra c*****plot the observed spectrum - if (plotopt .eq. 2) then + if (plotopt == 2) then call defcolor (1) - if (choice.eq.'h' .or. choice.eq.'f') then + if (choice=='h' .or. choice=='f') then call sm_lweight (4.0) else call sm_lweight (2.2) @@ -146,10 +146,10 @@ c*****plot the observed spectrum style(1) = 43.5 call sm_ptype (style,1) mount = lim2obs - lim1obs + 1 - if (mount .lt. 500) then + if (mount < 500) then call sm_points (xobs(lim1obs),yobs(lim1obs),mount) else - if (histoyes .eq. 1) then + if (histoyes == 1) then call sm_histogram (xobs(lim1obs),yobs(lim1obs),mount) else call sm_connect (xobs(lim1obs),yobs(lim1obs),mount) @@ -157,7 +157,7 @@ c*****plot the observed spectrum endif call sm_lweight (2.2) call sm_expand (0.7) - if (iflip .eq. 1) then + if (iflip == 1) then call sm_relocate (xhi+0.05*(xlo-xhi),ylo+0.12*(yhi-ylo)) else call sm_relocate (xlo+0.05*(xhi-xlo),ylo+0.12*(yhi-ylo)) @@ -165,7 +165,7 @@ c*****plot the observed spectrum call sm_label (obsitle) endif do i=1,2 - if (iflip .eq. 1) then + if (iflip == 1) then call sm_relocate (xhi+0.05*(xlo-xhi), . ylo+0.01+0.06*(2-i)*(yhi-ylo)) else @@ -174,41 +174,41 @@ c*****plot the observed spectrum endif call sm_label (modbin(i)) enddo - if (whichwin.eq.'1of1' .or. plotopt.ne.2) then + if (whichwin=='1of1' .or. plotopt/=2) then return endif c*****this section of code is executed only if a deviations plot is desired; c find the starting and stopping points in the arrays for the deviations - if (xsyn(kount) .le. xobs(lim1obs)) go to 1000 - if (xsyn(1) .gt. xobs(lim2obs)) go to 1000 - if (xsyn(1) .gt. xobs(lim1obs)) go to 150 + if (xsyn(kount) <= xobs(lim1obs)) go to 1000 + if (xsyn(1) > xobs(lim2obs)) go to 1000 + if (xsyn(1) > xobs(lim1obs)) go to 150 lim3obs = lim1obs do k=2,kount - if (xsyn(k) .gt. xobs(lim3obs)) then + if (xsyn(k) > xobs(lim3obs)) then lim1syn = k - 1 go to 155 endif enddo 150 lim1syn = 1 do l=lim1obs,lim2obs - if (xsyn(lim1syn) .le. xobs(l)) then + if (xsyn(lim1syn) <= xobs(l)) then lim3obs = l go to 155 endif enddo -155 if (xsyn(kount) .lt. xobs(lim2obs)) go to 160 +155 if (xsyn(kount) < xobs(lim2obs)) go to 160 lim4obs = lim2obs do k=lim1syn,kount - if (xsyn(k) .gt. xobs(lim4obs)) then + if (xsyn(k) > xobs(lim4obs)) then lim2syn = k go to 165 endif enddo 160 lim2syn = kount do l=lim3obs,lim2obs - if (xsyn(lim2syn) .lt. xobs(l)) then + if (xsyn(lim2syn) < xobs(l)) then lim4obs = l - 1 go to 165 endif @@ -221,7 +221,7 @@ c of the synthetic spectra is considered sufficient lpoint = lim1syn devsigma = 0. do i=lim3obs,lim4obs -170 if (xsyn(lpoint+1) .lt. xobs(i)) then +170 if (xsyn(lpoint+1) < xobs(i)) then lpoint = lpoint + 1 go to 170 endif @@ -235,7 +235,7 @@ c of the synthetic spectra is considered sufficient c from first set of deviations, define the plot limits, draw and label box - if (j .eq. 1) then + if (j == 1) then yup = -1000. ydown = +1000. do i=lim3obs,lim4obs @@ -265,12 +265,12 @@ c from first set of deviations, define the plot limits, draw and label box c plot the array of deviations - if (choice.eq.'h' .or. choice.eq.'f' .or. - . choice.eq.'g') then + if (choice=='h' .or. choice=='f' .or. + . choice=='g') then call defcolor (8) call sm_ltype (j-1) else - if (smterm(1:3) .eq. 'x11') then + if (smterm(1:3) == 'x11') then call defcolor (j+1) call sm_ltype (0) else @@ -289,7 +289,7 @@ c plot the array of deviations call sm_draw(xhi-0.215*(xhi-xlo), . ydown+(0.10+0.06*j)*(yup-ydown)) call sm_label (array) - if (choice .eq. 'g') then + if (choice == 'g') then noff = 80*(j-1) write (nf6out,3002) abitle(noff+1:noff+80) write (nf6out,3003) devsigma, velsh @@ -298,7 +298,7 @@ c plot the array of deviations c reset the spectrum plot boundaries before exiting - if(xlo .lt. xhi) then + if(xlo < xhi) then call sm_limits (xlo,xhi,ylo,yhi) iflip = 0 else diff --git a/Binplotprep.f b/Binplotprep.f index 78c8ce5..c9ca462 100755 --- a/Binplotprep.f +++ b/Binplotprep.f @@ -35,7 +35,7 @@ c*****open the files nchars = 40 call infile ('output ',nf9out,'formatted ',0,nchars, . f9out,lscreen) - if (plotopt .ne. 0) then + if (plotopt /= 0) then nf10out = 24 lscreen = lscreen + 2 array = 'SMOOTHED SYNTHESIS OUTPUT FOR COMBINED BINARY' @@ -49,7 +49,7 @@ c*****open the files call infile ('output ',nf5out,'formatted ',0,nchars, . f5out,lscreen) endif - if (plotopt .gt. 1) then + if (plotopt > 1) then nf6out = 27 lscreen = lscreen + 2 array = 'SPECTRUM COMPARISON OUTPUT' @@ -86,19 +86,19 @@ c*****read back the header information from the individual star raw c synthetic spectra 99 read (nf7out,1004,end=100) array read (nf8out,1004) chinfo - if (array(1:7).eq.'Isotopi') then + if (array(1:7)=='Isotopi') then write (nf9out,1004) array go to 99 - elseif (array(1:7).eq.'ALL abu') then + elseif (array(1:7)=='ALL abu') then write (nf9out,1010) array(1:59), chinfo(54:59) go to 99 - elseif (array(1:7).eq.'Changin') then + elseif (array(1:7)=='Changin') then write (nf9out,1011) array(1:37), chinfo(31:37) go to 99 - elseif (array(1:7).eq.'element') then + elseif (array(1:7)=='element') then write (nf9out,1012) array(1:30), chinfo(25:30) go to 99 - elseif (array(1:7).eq.'MODEL: ') then + elseif (array(1:7)=='MODEL: ') then write (nf9out,1013) array(8:43), chinfo(8:43) modbin(1) = array(1:80) modbin(2) = chinfo(1:80) @@ -120,14 +120,14 @@ c and dump to combined raw synthetic spectrum file z(i) = (1.0 - z(i))/lumratio enddo normalization = 1.0 + 1.0/lumratio - if (pshift .gt. 0) then + if (pshift > 0) then do i=1,pshift dev(i) = 1.0 enddo do i=pshift+1,kount dev(i) = (y(i)+z(i-pshift))/normalization enddo - elseif (pshift .eq. 0) then + elseif (pshift == 0) then do i=1,kount dev(i) = (y(i)+z(i))/normalization enddo @@ -68,7 +68,7 @@ c*****start the large loop that will go through each blended feature ewsynthopt = -1 mode = 4 30 do lll=1,1000 - if (lim2line .eq. nlines) exit + if (lim2line == nlines) exit c*****define the set of lines responsible for a blended feature @@ -82,26 +82,26 @@ c*****make sure that the element whose abundance is to be fit has c a representative line of the blend ifind = 0 do j=lim1,lim2 - if (atom1(j) .lt. 100.) then - if (iatom .eq. int(atom1(j))) then + if (atom1(j) < 100.) then + if (iatom == int(atom1(j))) then ifind = 1 exit endif else call sunder (atom1(j),ia,ib) - if (iatom.eq.ia .or. iatom.eq.ib) then + if (iatom==ia .or. iatom==ib) then ifind = 1 exit endif endif enddo - if (ifind .eq. 0) then + if (ifind == 0) then do j=lim1,lim2 abundout(j) = 999.99 enddo write (nf1out,1002) lim1line = lim2line + 1 - if (lim1line .le. nlines+nstrong) cycle + if (lim1line <= nlines+nstrong) cycle endif @@ -127,31 +127,31 @@ c arbitrarily c*****here we go for another iteration - if (dabs(error) .ge. 0.0075) then + if (dabs(error) >= 0.0075) then rwlcomp = dlog10(w(ncurve-1)/wave1(lim1)) - if (rwlcomp.lt.-5.2 .and. rwlgobs.lt.-5.2) then + if (rwlcomp<-5.2 .and. rwlgobs<-5.2) then ratio = ratio - elseif (rwlcomp.ge.-5.2 .and. rwlgobs.ge.-5.2) then + elseif (rwlcomp>=-5.2 .and. rwlgobs>=-5.2) then ratio = ratio**2.0 else ratio = ratio**1.5 endif gf1(ncurve) = gf1(ncurve-1)*ratio do j=lim1,lim2 - if (atom1(j) .gt. 100.) then + if (atom1(j) > 100.) then call sunder (atom1(j),ia,ib) - if (ia.eq.iatom .or. ib.eq.iatom) then + if (ia==iatom .or. ib==iatom) then do i=1,ntau kapnu0(j,i) = kapnu0(j,i)*ratio enddo endif - elseif (int(atom1(j)) .eq. iatom) then + elseif (int(atom1(j)) == iatom) then do i=1,ntau kapnu0(j,i) = kapnu0(j,i)*ratio enddo endif enddo - if (k .eq. 20) then + if (k == 20) then write (*,1008) stop endif @@ -164,14 +164,14 @@ c*****here we go for another iteration c*****here we do the final calculation when the predicted and observed are close gf1(ncurve) = gf1(ncurve-1)*ratio do j=lim1,lim2 - if (atom1(j) .gt. 100.) then + if (atom1(j) > 100.) then call sunder (atom1(j),ia,ib) - if (ia.eq.iatom .or. ib.eq.iatom) then + if (ia==iatom .or. ib==iatom) then do i=1,ntau kapnu0(j,i) = kapnu0(j,i)*ratio enddo endif - elseif (int(atom1(j)) .eq. iatom) then + elseif (int(atom1(j)) == iatom) then do i=1,ntau kapnu0(j,i) = kapnu0(j,i)*ratio enddo @@ -182,13 +182,13 @@ c*****here we do the final calculation when the predicted and observed are close widout(lim1) = w(ncurve) diff = dlog10(gf1(ncurve)) abundout(lim1) = dlog10(xabund(iatom)) + 12.0 + diff - if (ncurve .ne. 1) then + if (ncurve /= 1) then write (nf1out,1001) ncurve endif c*****here is where some auxiliary things like mean depth are computed - if (lim1.eq.lim2 .and. linprintopt.ge.3) then + if (lim1==lim2 .and. linprintopt>=3) then wave = wave1(lim1) call taukap call cdcalc(2) @@ -201,11 +201,11 @@ c*****here is where some auxiliary things like mean depth are computed cdmean = rinteg(xref,dummy1,dummy2,ntau,first)/ . rinteg(xref,cd,dummy2,ntau,first) do i=1,ntau - if (cdmean .lt. cd(i)) exit + if (cdmean < cd(i)) exit enddo write (nf1out,1005) lim1, cdmean, i, xref(i) do i=1,ntau - if (taunu(i)+taulam(i) .ge. 1.) exit + if (taunu(i)+taulam(i) >= 1.) exit enddo write (nf1out,1006) lim1, i, dlog10(tauref(i)), . dlog10(taulam(i)), dlog10(taunu(i)) @@ -215,7 +215,7 @@ c*****here is where some auxiliary things like mean depth are computed c*****assign the abundance to the strongest line of the blend that c contains cogatom; put 999.99's as the abundances of all but c this line; go back for another blended feature - if (lim2 .gt. lim1) then + if (lim2 > lim1) then abunblend = abundout(lim1) widblend = widout(lim1) strongest = 0. @@ -224,17 +224,17 @@ c this line; go back for another blended feature abundout(j) = 999.99 enddo do j=lim1,lim2 - if (atom1(j) .lt. 100.) then - if (dint(atom1(j)) .eq. cogatom) then - if (kapnu0(j,jtau5) .gt. strongest) then + if (atom1(j) < 100.) then + if (dint(atom1(j)) == cogatom) then + if (kapnu0(j,jtau5) > strongest) then strongest = kapnu0(j,jtau5) linstrongest = j endif endif else call sunder (atom1(j),ia,ib) - if (dble(ia).eq.cogatom .or. dble(ib).eq.cogatom) then - if (kapnu0(j,jtau5) .gt. strongest) then + if (dble(ia)==cogatom .or. dble(ib)==cogatom) then + if (kapnu0(j,jtau5) > strongest) then strongest = kapnu0(j,jtau5) linstrongest = j endif @@ -245,7 +245,7 @@ c this line; go back for another blended feature widout(linstrongest) = widblend endif lim1line = lim2line + 1 - if (lim1line .le. nlines+nstrong) cycle + if (lim1line <= nlines+nstrong) cycle enddo @@ -260,23 +260,23 @@ c*****do abundance statistics; print out a summary of the abundances c*****here a plot may be made on the terminal (and paper) if there c are enough lines; then the user will be prompted on some c options concerning what is seen on the plot - if (plotopt .ne. 0) then + if (plotopt /= 0) then call pltabun - if (choice.eq.'v' .or. choice.eq.'m') go to 100 + if (choice=='v' .or. choice=='m') go to 100 endif c*****now the option will be to redo the molecular equilibrium and redo c the last species, at the user's option. - if (neq .ne. 0) then + if (neq /= 0) then do n=1,neq - if (iatom .eq. iorder(n)) then + if (iatom == iorder(n)) then write (array,1004) ikount = min0(nlines+11,maxline) nchars = 70 35 call getasci (nchars,ikount) choice = chinfo(1:1) - if (choice.eq.'y' .or. nchars.le.0) then + if (choice=='y' .or. nchars<=0) then xabund(iatom) = 10.**(average-12.) call eqlib call nearly (1) @@ -284,7 +284,7 @@ c the last species, at the user's option. lim2line = 0 rewind nf2out go to 30 - elseif (choice .eq. 'n') then + elseif (choice == 'n') then call finish (0) else go to 35 @@ -15,16 +15,16 @@ c****************************************************************************** c*****if only a synthetic spectrum is being plotted, use beginning and c end points of the synthetic spectrum instead of those of the observed c spectrum - if (plotopt .lt. 2) then + if (plotopt < 2) then do i=1,kount - if (xlo .le. xsyn(i)) then + if (xlo <= xsyn(i)) then lim1obs = i go to 10 endif enddo lim1obs = kount 10 do i=lim1obs,kount - if (xhi .lt. xsyn(i)) then + if (xhi < xsyn(i)) then lim2obs = i -1 return endif @@ -36,14 +36,14 @@ c spectrum c*****and here is the same logic when an observed spectrum exists else do i=1,lount - if (xlo .le. xobs(i)) then + if (xlo <= xobs(i)) then lim1obs = i go to 20 endif enddo lim1obs = lount 20 do i=lim1obs,lount - if (xhi .lt. xobs(i)) then + if (xhi < xobs(i)) then lim2obs = i -1 return endif @@ -10,11 +10,11 @@ c****************************************************************************** include 'Linex.com' c*****continuum "contribution curve" calculation - if (number .eq. 1) then + if (number == 1) then do i=1,ntau scont(i) = ((1.19089d+25/wave**2)*1.0d+10)/(wave**3* . (dexp(1.43879d+08/(wave*t(i)))-1.0d+00)) - if (fluxintopt .eq. 1) then + if (fluxintopt == 1) then cd(i) = kaplam(i)*tauref(i)*scont(i)* . dexp(-taulam(i))/(0.4343*kapref(i)) else @@ -27,8 +27,8 @@ c*****line plus continuum "contribution curve" calculation else do i=1,ntau sline(i) = scont(i) - if (fluxintopt .eq. 1) then - if (taulam(i)+taunu(i) .le. 50.) then + if (fluxintopt == 1) then + if (taulam(i)+taunu(i) <= 50.) then exptau = dexp(-taulam(i)-taunu(i)) else exptau = 0. @@ -22,7 +22,7 @@ c*****set the local temporary parameters to what they were for the last synth newnumisosyn = numisosyn newnumpecatom = numpecatom newnumatomsyn = numatomsyn - if (newnumatomsyn .eq. 0) newnumatomsyn = 1 + if (newnumatomsyn == 0) newnumatomsyn = 1 do i=3,95 newpec(i) = pec(i) do j=1,newnumatomsyn @@ -30,23 +30,23 @@ c*****set the local temporary parameters to what they were for the last synth enddo enddo do i=1,newnumiso - if(newnumiso .ge. 1) then + if(newnumiso >= 1) then newisotope(i)=isotope(i) do j=1,newnumisosyn newisoabund(i,j)=isoabund(i,j) enddo endif enddo - if (newnumatomsyn .le. 0) newnumatomsyn = 1 + if (newnumatomsyn <= 0) newnumatomsyn = 1 c*****SPECIAL CASE c for the "binary" driver, present a limited set of options; only c the abundances of the elements already chosen for variation in c the parameter file can be altered -11 if (control .eq. 'binary ') then +11 if (control == 'binary ') then istat = ivcleof(4,1) - if (syncount .eq. 1) then + if (syncount == 1) then array = 'Abundance Alterations for the PRIMARY STAR: ' else array = 'Abundance Alterations for the SECONDARY STAR: ' @@ -56,7 +56,7 @@ c the parameter file can be altered istat = ivwrite(6,1,array,49) line = 6 do i=3,95 - if(newpec(i) .eq. 1) then + if(newpec(i) == 1) then line = line + 1 write (array,1001) i,(newpecabund(i,j),j=1,newnumatomsyn) istat = ivwrite(line,1,array,57) @@ -70,16 +70,16 @@ c the parameter file can be altered nchars = 21 call getasci (nchars,line+4) choice = chinfo(1:1) - if (choice.ne.'c' .and. choice.ne.'x' .and. - . choice.ne.'q') go to 11 + if (choice/='c' .and. choice/='x' .and. + . choice/='q') go to 11 endif c*****for the "synth" driver, present the user with the current abundance c alterations, and the options; find out what is desired -1 if (control .eq. 'synth' .or. - . control .eq. 'isotop' .or. - . control .eq. 'synpop') then +1 if (control == 'synth' .or. + . control == 'isotop' .or. + . control == 'synpop') then line = 4 istat = ivcleof(line,1) write (array,1006) @@ -87,13 +87,13 @@ c alterations, and the options; find out what is desired istat = ivwrite(6,1,array,72) line = line + 1 do i=3,95 - if(newpec(i) .eq. 1) then + if(newpec(i) == 1) then line = line + 1 write (array,1001) i,(newpecabund(i,j),j=1,newnumatomsyn) istat = ivwrite(line,1,array,57) endif enddo - if (newnumiso .gt. 0) then + if (newnumiso > 0) then do i=1,newnumiso line = line + 1 write (array,1002) i, newisotope(i), @@ -109,32 +109,32 @@ c alterations, and the options; find out what is desired nchars = 21 call getasci (nchars,line+3) choice = chinfo(1:1) - if (choice.ne.'c' .and. choice.ne.'i' .and. - . choice.ne.'n' .and. choice.ne.'x' .and. - . choice.ne.'q') go to 1 + if (choice/='c' .and. choice/='i' .and. + . choice/='n' .and. choice/='x' .and. + . choice/='q') go to 1 endif c*****for "synth", change elemental abundances - if (choice .eq. 'c') then + if (choice == 'c') then 20 istat = ivcleof(line+1,1) array = 'Which element to change? ' nchars = 25 call getnum (nchars,line+1,xnum,shortnum) - if (xnum.lt.2.0 .or. xnum.gt.95.) go to 20 + if (xnum<2.0 .or. xnum>95.) go to 20 j = nint(xnum) - if (newpec(j) .eq. 1) then + if (newpec(j) == 1) then 25 array = 'n = new abundances, or z = zero offsets? ' nchars = 41 call getasci (nchars,line+2) choice2 = chinfo(1:1) - if (choice2 .eq. 'z') then + if (choice2 == 'z') then newpec(j) = 0 do i=1,5 newpecabund(j,i) = 0.0 enddo newnumpecatom = newnumpecatom - 1 - elseif (choice2 .eq. 'n') then + elseif (choice2 == 'n') then write (array,1004) istat = ivwrite(line+3,1,array,40) read (*,*) (newpecabund(j,i),i=1,newnumatomsyn) @@ -151,7 +151,7 @@ c*****for "synth", change elemental abundances c*****for "synth", change isotopic factors - elseif (choice .eq. 'i') then + elseif (choice == 'i') then istat = ivcleof(line+1,1) array = 'Options: c = change an isotopic factor' istat = ivwrite(line+1,1,array,49) @@ -161,18 +161,18 @@ c*****for "synth", change isotopic factors nchars = 22 call getasci (nchars,line+3) choice2 = chinfo(1:1) - if (choice2 .eq. 'c') then + if (choice2 == 'c') then 35 istat = ivcleof(line+1,1) array = 'Which isotope number from the list? ' nchars = 36 call getnum (nchars,line+1,xnum,shortnum) j = nint(xnum) - if (j.lt.1 .or. j.gt.newnumiso) go to 35 + if (j<1 .or. j>newnumiso) go to 35 istat = ivcleof(line+2,1) array = 'What are the new division factors? ' istat = ivwrite(line+2,1,array,35) read (*,*) (newisoabund(j,i),i=1,newnumisosyn) - elseif (choice2 .eq. 'n') then + elseif (choice2 == 'n') then newnumiso = newnumiso + 1 istat = ivcleof(line+1,1) array = 'What is the new isotope designation? ' @@ -188,24 +188,24 @@ c*****for "synth", change isotopic factors c*****for "synth", change the number of syntheses - elseif (choice .eq. 'n') then + elseif (choice == 'n') then 55 array = 'How many synths? ' nchars = 17 call getnum (nchars,line+5,xnum,shortnum) - if (xnum .gt. 5.) go to 55 + if (xnum > 5.) go to 55 newnumatomsyn = nint(xnum) go to 1 c*****for "synth", exit the routine without changing anything - elseif (choice .eq. 'x') then + elseif (choice == 'x') then return c*****for "synth", make the proposed alterations permanent; then c return to the calling routine, which allegedly will c redo the syntheses. - elseif (choice .eq. 'q') then + elseif (choice == 'q') then numiso = newnumiso numisosyn = newnumisosyn numpecatom = newnumpecatom @@ -217,7 +217,7 @@ c redo the syntheses. enddo enddo do i=1,numiso - if(numiso .ge. 1) then + if(numiso >= 1) then isotope(i) = newisotope(i) do j=1,numisosyn isoabund(i,j)=newisoabund(i,j) @@ -229,7 +229,7 @@ c redo the syntheses. c*****loop back and print out the main menu again - if (control .eq. 'synth ') then + if (control == 'synth ') then go to 1 else go to 11 @@ -61,7 +61,7 @@ c*****open and read the line list file; get ready for the line calculations c*****define the range of lines (the whole list, in this case) mode = 1 call linlimit - if (lim1line .lt. 0) then + if (lim1line < 0) then call finish (0) return endif @@ -72,7 +72,7 @@ c*****do the curves of growth, making plots if desired lim2 = lim1 call curve call pltcog - if (choice .eq. 'm') then + if (choice == 'm') then close (unit=nfmodel) close (unit=nflines) rewind nf1out @@ -83,13 +83,13 @@ c*****do the curves of growth, making plots if desired fmodel = 'no_filename_given' go to 102 endif - if (choice .eq. 'v') go to 101 + if (choice == 'v') go to 101 array = 'DO ANOTHER CURVE-OF-GROWTH ([y]/n)? ' nchars = 36 lscreen = 16 call getasci (nchars,lscreen) choice = chinfo(1:1) - if (choice .eq. 'n') then + if (choice == 'n') then call finish (0) return endif @@ -25,8 +25,8 @@ c*****define the plot boundaries ylo = real(nint(rwlow*10))/10. - 0.1 yhi = real(nint(rwhigh*10))/10. + 0.1 do i=1,ncurve - if (rwplot(i) .gt. rwlow) then - if (gfplot(i) .gt. 0) then + if (rwplot(i) > rwlow) then + if (gfplot(i) > 0) then xlo = real(int(gfplot(i)*10))/10 else xlo = real(int(gfplot(i)*10))/10 - 0.1 @@ -35,8 +35,8 @@ c*****define the plot boundaries endif enddo 10 do i=1,ncurve - if (rwplot(i) .gt. rwhigh) then - if (gfplot(i) .gt. 0) then + if (rwplot(i) > rwhigh) then + if (gfplot(i) > 0) then xhi = real(int(gfplot(i)*10))/10 + 0.1 else xhi = real(int(gfplot(i)*10))/10 @@ -83,11 +83,11 @@ c*****plot the computed curve-of-growth points; exit normally call sm_points (gfplot,rwplot,ncurve) call defcolor (1) ich = idint(charge(lim1) + 0.1) - if (ich .eq. 1) then + if (ich == 1) then ion = ' I ' - elseif (ich .eq. 2) then + elseif (ich == 2) then ion = ' II ' - elseif (ich .eq. 3) then + elseif (ich == 3) then ion = ' III' endif iatom = idint(atom1(lim1)) @@ -67,13 +67,13 @@ c*****do the syntheses nlines = 0 call synspec call total - if (ncurve .eq. 1) then + if (ncurve == 1) then wstart = 10.**(rwlow)*wave1(lim1) wstop = 10.**(rwhigh)*wave1(lim1) endif molopt = 1 linprintopt = 0 - if (w(ncurve) .gt. wstart) then + if (w(ncurve) > wstart) then pecabund(iatom,1) = pecabund(iatom,1) - rwstep go to 10 endif @@ -84,7 +84,7 @@ c*****do the syntheses call total molopt = 1 linprintopt = 0 - if (w(ncurve) .lt. wstop) then + if (w(ncurve) < wstop) then pecabund(iatom,1) = pecabund(iatom,1) + rwstep go to 20 endif @@ -26,13 +26,13 @@ c****************************************************************************** c*****find the array subscripts for the minimum and maximum c wavelengths in the synthetic spectrum. do i=1,kount - if (xsyn(i) .ge. wavemin) then + if (xsyn(i) >= wavemin) then ilosyn = i go to 5 endif enddo 5 do i=ilosyn,kount - if (xsyn(i) .gt. wavemax) then + if (xsyn(i) > wavemax) then ihisyn = i - 1 go to 10 endif @@ -51,7 +51,7 @@ c*****dump the synthetic spectrum wavelengths and fluxes into working arrays c*****find the array subscript for the minimum wavelength in the c observed spectrum do j=1,lount - if (xobs(j) .ge. xgood(1)) then + if (xobs(j) >= xgood(1)) then jloobs = j go to 15 endif @@ -63,7 +63,7 @@ c spectrum that matches the wavelength step size of the c synthetic spectrum 15 j = jloobs do i=1,itotsyn - if (xgood(i) .gt. (xobs(j+1)+xobs(j))/2.) j = j + 1 + if (xgood(i) > (xobs(j+1)+xobs(j))/2.) j = j + 1 q = xgood(i) zgood(i) = yobs(j-1)*(q-xobs(j))*(q-xobs(j+1))/ . ((xobs(j-1)-xobs(j))*(xobs(j-1)-xobs(j+1))) + @@ -89,12 +89,12 @@ c*****interpolate to find the maximum of the correlation function corrmax = spy(1) mc = 1 do i=2,imax - if (spy(i) .gt. corrmax) then + if (spy(i) > corrmax) then corrmax = spy(i) mc = i endif enddo - if (mc.eq.1 .or. mc.eq.imax) then + if (mc==1 .or. mc==imax) then deltawave = 0. return endif @@ -113,7 +113,7 @@ c*****interpolate to find the maximum of the correlation function jmax = 1 corrmax = yinterp(1) do j=2,21 - if (yinterp(j) .gt. corrmax) then + if (yinterp(j) > corrmax) then jmax = j corrmax = yinterp(j) shiftmax = xinterp(j) diff --git a/Crosscorr.f b/Crosscorr.f index 5e33d12..457c6e8 100755 --- a/Crosscorr.f +++ b/Crosscorr.f @@ -17,7 +17,7 @@ c****************************************************************************** xysum = 0. x2sum = 0. y2sum = 0. - if (ishift .ge. 0) then + if (ishift >= 0) then minpt = 1 maxpt = npx - ishift else @@ -22,7 +22,7 @@ c*****set up the parameters c*****increment the log(gf) value backward by rwstep and redo the calculations c until the end (rwlow) is reached 31 call oneline (2) - if (w(ncurve) .gt. wstart) then + if (w(ncurve) > wstart) then gf1(ncurve) = gf1(ncurve)/dec do i=1,ntau kapnu0(lim1,i) = kapnu0(lim1,i)/dec @@ -39,26 +39,26 @@ c until the end (rwhigh) is reached kapnu0(lim1,i) = kapnu0(lim1,i)*dec enddo 61 call oneline (2) - if (w(ncurve) .lt. wstop) then + if (w(ncurve) < wstop) then ncurve = ncurve + 1 go to 60 endif c*****end the computations with a summary print - if (nf2out .ne. 0 .and. lim1 .eq. 1) + if (nf2out /= 0 .and. lim1 == 1) . write (nf2out,1001) moditle do i=1,ncurve w(i) = dlog10(w(i)/wave1(lim1)) gf1(i) = dlog10(gf1(i)) enddo iatom = atom1(lim1) - if(iatom .ge. 100) iatom = 1 + if(iatom >= 100) iatom = 1 abund = dlog10(xabund(iatom)) + 12. write (nf1out,1002) wave1(lim1),atom1(lim1),e(lim1,1), . abund,ncurve write (nf1out,1003) (gf1(i),w(i),i=1,ncurve) - if (nf2out .eq. 0) return + if (nf2out == 0) return write (nf2out,1002) wave1(lim1),atom1(lim1),e(lim1,1), . abund,ncurve write (nf2out,1003) (gf1(i),w(i),i=1,ncurve) @@ -13,43 +13,43 @@ c****************************************************************************** j = linnumber iwave = int(wave1(j)) iatom10 = nint(10.*atom1(j)) - if (dampnum(j) .lt. 0.) dampnum(j) = 10.**dampnum(j) + if (dampnum(j) < 0.) dampnum(j) = 10.**dampnum(j) c*****for a few lines, explicit detailed broadening terms have c appeared in the literature, and so do these lines with a c sepaarate subroutine - if (itru .eq. 0) then + if (itru == 0) then c Ca II - if (iatom10 .eq. 201) then - if (iwave .eq. 8498 .or. - . iwave .eq. 8542 .or. - . iwave .eq. 8662 .or. - . iwave .eq. 3933) then + if (iatom10 == 201) then + if (iwave == 8498 .or. + . iwave == 8542 .or. + . iwave == 8662 .or. + . iwave == 3933) then call trudamp (j) damptype(j) = 'TRUEgam' return endif c CH - elseif(iatom10 .eq. 1060) then - if (iwave .eq. 3693) then + elseif(iatom10 == 1060) then + if (iwave == 3693) then call trudamp (j) damptype(j) = 'TRUEgam' return endif c Ca I - elseif (iatom10 .eq. 200) then - if (iwave.eq.6717 .or. iwave.eq.6318 - . .or. iwave.eq.6343 .or. iwave.eq.6361) then + elseif (iatom10 == 200) then + if (iwave==6717 .or. iwave==6318 + . .or. iwave==6343 .or. iwave==6361) then call trudamp (j) damptype(j) = 'TRUEgam' return endif c Ca I autoionization - elseif (iatom10 .eq. 200) then - if (iwave.eq.6318 .or. - . iwave.eq.6343 .or. - . iwave.eq.6361) then + elseif (iatom10 == 200) then + if (iwave==6318 .or. + . iwave==6343 .or. + . iwave==6361) then call trudamp (j) damptype(j) = 'TRUEgam' return @@ -82,19 +82,19 @@ c c6 done as in dampingopt = 0 c*****these damping calculations are done at each atmosphere level - if (linprintopt .gt. 2) write (nf1out,1001) j, wave1(j) + if (linprintopt > 2) write (nf1out,1001) j, wave1(j) do i=1,ntau ich = nint(charge(j)) v1 = dsqrt(2.1175d8*t(i)*(1.0/amass(j)+1.008)) c*****first calculate an Unsold approximation to gamma_VanderWaals - if (atom1(j) .gt. 100.) then + if (atom1(j) > 100.) then ebreakup = 7.0 else ebreakup = chi(j,ich) endif - if (e(j,1).ge.ebreakup .or. e(j,2).ge.ebreakup) then + if (e(j,1)>=ebreakup .or. e(j,2)>=ebreakup) then unsold = 1.0e-33 else unsold = dabs(1.61d-33*(13.598*charge(j)/(ebreakup - @@ -105,15 +105,15 @@ c*****first calculate an Unsold approximation to gamma_VanderWaals c*****dampingopt = 0 or c*****dampingopt = 1 and no Barklem data - if (dampingopt .eq. 0 .or. - . (dampingopt.eq.1 .and. gambark(j).lt.0)) then - if (dampnum(j) .eq. 0.0) then + if (dampingopt == 0 .or. + . (dampingopt==1 .and. gambark(j)<0)) then + if (dampnum(j) == 0.0) then damptype(j) = 'UNSLDc6' gammav = 17.0*unsold**0.4*v1**0.6*numdens(1,1,i) - elseif (dampnum(j) .lt. 1.0d-15) then + elseif (dampnum(j) < 1.0d-15) then damptype(j) = ' MYc6' gammav = 17.0*dampnum(j)**0.4*v1**0.6*numdens(1,1,i) - elseif (dampnum(j) .lt. 1.0d-04) then + elseif (dampnum(j) < 1.0d-04) then damptype(j) = 'MYgamma' gammav = dampnum(j)*(t(i)/10000.)**0.3*numdens(1,1,i) else @@ -124,23 +124,23 @@ c*****dampingopt = 1 and no Barklem data c*****dampingopt = 1 with extant Barklem data - elseif (dampingopt.eq.1 .and. gambark(j).gt.0.) then + elseif (dampingopt==1 .and. gambark(j)>0.) then damptype(j) = 'BKgamma' gammav = . gambark(j)*(t(i)/10000.)**alpbark(j)*numdens(1,1,i) c*****dampingopt = 2 - elseif (dampingopt .eq. 2) then + elseif (dampingopt == 2) then damptype(j) = 'BLKWLc6' gammav = 17.0*((1.0 + 0.67*e(j,1))*unsold)**0.4* . v1**0.6*numdens(1,1,i) c*****dampingopt = 3 - elseif (dampingopt .eq. 3) then + elseif (dampingopt == 3) then damptype(j) = 'NXTGNc6' - if (dampnum(j) .le. 1.0d-10) dampnum(j) = 1.0 + if (dampnum(j) <= 1.0d-10) dampnum(j) = 1.0 c6h = dabs(1.01d-32*(charge(j)**2)* . (13.598/(ebreakup - e(j,1)))**2 - 1.61d-33* . (13.598/(ebreakup-e(j,2)))**2) @@ -160,7 +160,7 @@ c*****dampingopt = 3 c*****compute radiative broadening either by an approximate formula or c*****the value in Barklem.dat) - if (gamrad(j).ne.0.0 .and. dampingopt .eq. 1) then + if (gamrad(j)/=0.0 .and. dampingopt == 1) then gammar = gamrad(j) else gammar = 2.223d15/wave1(j)**2 @@ -169,7 +169,7 @@ c*****the value in Barklem.dat) c*****now Stark broadening (approximate formulae) excdiff = chi(j,nint(charge(j))) - e(j,2) - if (excdiff .gt. 0.0 .and. atom1(j).lt.100.) then + if (excdiff > 0.0 .and. atom1(j)<100.) then effn2 = 13.6*charge(j)**2/excdiff else effn2 = 25. @@ -180,7 +180,7 @@ c*****now Stark broadening (approximate formulae) c*****now finish by summing the gammas and computing the Voigt *a* values gammatot = gammar + gammas + gammav a(j,i) = gammatot*wave1(j)*1.0d-8/(12.56636*dopp(j,i)) - if (linprintopt .gt. 2) write (nf1out,1002) i, gammar, + if (linprintopt > 2) write (nf1out,1002) i, gammar, . gammas, gammav, gammatot, a(j,i) enddo return @@ -20,8 +20,8 @@ c*****assign colors to character arrays colors(8) = 'black ' - if (choice.eq.'h' .or. choice.eq.'f' .or. - . choice.eq.'g') then + if (choice=='h' .or. choice=='f' .or. + . choice=='g') then call sm_ctype (colors(8)) else call sm_ctype (colors(icolor)) @@ -12,7 +12,7 @@ c****************************************************************************** im = amol do i=1,5 i3 = im/itest(i) - if(i3 .eq. i1) i2 = i2 + 1 + if(i3 == i1) i2 = i2 + 1 im = im - i3*itest(i) enddo return @@ -48,26 +48,26 @@ c*****open and read the model atmosphere c*****compute the flux curve wave = start 1 call opacit (2,wave) - if (modprintopt .ge. 2) + if (modprintopt >= 2) . write(nf1out,1002) wave,(kaplam(i),i=1,ntau) call cdcalc (1) first = 0.4343*cd(1) flux = rinteg(xref,cd,dummy1,ntau,first) - if (flux .le. 0.1) flux = 0. - if (iunits .eq. 1) then + if (flux <= 0.1) flux = 0. + if (iunits == 1) then write (nf1out,1003) 1.d-4*wave,flux else write (nf1out,1004) wave,flux endif waveinv = 1.0d4/wave - if (flux .gt. 0.) then + if (flux > 0.) then fluxlog = dlog10(flux) else fluxlog = -1.0 endif write (nf2out,1001) wave, flux, waveinv, fluxlog wave = wave + step - if (wave .le. sstop) go to 1 + if (wave <= sstop) go to 1 call pltflux @@ -10,9 +10,9 @@ c****************************************************************************** call sm_graphics - if (whichwin .eq. '1of1') then + if (whichwin == '1of1') then call sm_window (1,1,1,1,1,1) - elseif (whichwin .eq. '2of2') then + elseif (whichwin == '2of2') then call sm_defvar ('y_gutter','0.0') call sm_window (1,2,1,1,1,1) endif @@ -29,11 +29,11 @@ c****************************************************************************** call sm_relocate (xplotpos,yplotpos-0.11*(yhi-ylo)) call sm_putlabel (5,array) call writenumber (yplotpos) - if (whichwin(4:4) .eq. '1') then + if (whichwin(4:4) == '1') then call sm_relocate (xplotpos,yplotpos-0.15*(yhi-ylo)) - elseif (whichwin(4:4) .eq. '2') then + elseif (whichwin(4:4) == '2') then call sm_relocate (xplotpos,yplotpos-0.18*(yhi-ylo)) - elseif (whichwin(4:4) .eq. '3') then + elseif (whichwin(4:4) == '3') then call sm_relocate (xplotpos,yplotpos-0.21*(yhi-ylo)) endif call sm_putlabel (5,array) @@ -33,9 +33,9 @@ c*****clear the arrays c*****the number of species to be considered has been defined in "Inmodel"; c*****either read in the dissociation data for a molecular species do jmol=1,nmol - if (amol(jmol) .ge. 100.) then + if (amol(jmol) >= 100.) then do k=1,110 - if (datmol(1,k) .eq. amol(jmol)) go to 11 + if (datmol(1,k) == amol(jmol)) go to 11 enddo write (nf1out,1001) amol(jmol) stop @@ -54,7 +54,7 @@ c*****or read the ionization data for an atomic species it = ti do jj=1,2 att = atom+0.1*(jj-1) - if (partflag(iatom1,jj) .gt. 0) then + if (partflag(iatom1,jj) > 0) then uu(jj) = partnew(att,jj,it) else uu(jj) = ucalc(att,it) @@ -64,7 +64,7 @@ c*****or read the ionization data for an atomic species enddo endif enddo - if (molopt .ge. 2) + if (molopt >= 2) . write (nf1out,1002) (amol(i),(const(j,i),j=1,6),i=1,nmol) @@ -77,24 +77,24 @@ c*****understood, but is not explicitly contained in 'ident'. atom = amol(jmol) 2 call sunder(atom,iatom1,iatom2) do k=1,30 - if (iatom1 .eq. iorder(k)) go to 4 + if (iatom1 == iorder(k)) go to 4 enddo neq = neq + 1 iorder(neq) = iatom1 ident(neq,1) = jmol go to 6 4 do kk=1,nmax - if (ident(k,kk).eq.0 .or. ident(k,kk).eq. jmol) go to 7 + if (ident(k,kk)==0 .or. ident(k,kk)== jmol) go to 7 enddo nmax = nmax + 1 kk = nmax 7 ident(k,kk) = jmol -6 if (iatom2 .ne. 0) then +6 if (iatom2 /= 0) then atom = iatom2 go to 2 endif enddo - if (molopt .ge. 2) then + if (molopt >= 2) then do i=1,neq dummy1(i) = iorder(i) enddo @@ -114,7 +114,7 @@ c*****calculate *xfic* and make a first guess at *xatom* korder = iorder(k) xfic(k) = xabund(korder)*nhtot(i) enddo - if (i .lt. ntau) then + if (i < ntau) then do k=1,neq xatom(k) = xatom(k)*nhtot(i)/nhtot(i+1) enddo @@ -149,8 +149,8 @@ c (1 eV = 1.60219E-12 ergs). Also, k = 1.38065E-16 erg/K, c h = 6.626076E-27 erg s, and pi = 3.1415926536. 27 do jmol=1,nmol atom = amol(jmol) - if (atom .ge. 100.) then - if (t(i) .gt. 12000.) then + if (atom >= 100.) then + if (t(i) > 12000.) then xmol(jmol,i) = 1.0d-20 else xmol(jmol,i) = 1. @@ -158,10 +158,10 @@ c h = 6.626076E-27 erg s, and pi = 3.1415926536. 37 call sunder(atom,iatom1,iatom2) count = count + 1. do k=1,neq - if (iorder(k) .eq. iatom1) + if (iorder(k) == iatom1) . xmol(jmol,i) = xmol(jmol,i)*xatom(k) enddo - if (iatom2 .ne. 0) then + if (iatom2 /= 0) then atom = iatom2 go to 37 endif @@ -185,7 +185,7 @@ c*****compute the number of ions: . (const(m+1,jmol)-const(m,jmol))*delt iatom1 = atom do k=1,neq - if (iorder(k) .eq. iatom1) xmol(jmol,i) = + if (iorder(k) == iatom1) xmol(jmol,i) = . 4.825d15*u1*t(i)**1.5/ne(i)*dexp(-1.1605d4* . const(1,jmol)/t(i))*xatom(k) enddo @@ -206,12 +206,12 @@ c*****respect to each atom. kderiv = iorder(kk) do 28 j=1,nmax jmol = ident(k,j) - if (jmol .eq. 0) go to 28 + if (jmol == 0) go to 28 call discov(amol(jmol),kderiv,num2) - if (num2 .eq. 0) go to 28 + if (num2 == 0) go to 28 call discov(amol(jmol),korder,num1) c(k,kk) = c(k,kk) + xmol(jmol,i)*num1*num2/xatom(kk) - if (k .eq. kk) deltax(k) = deltax(k) + num1*xmol(jmol,i) + if (k == kk) deltax(k) = deltax(k) + num1*xmol(jmol,i) 28 continue enddo enddo @@ -233,17 +233,17 @@ c*****decide if another iteration is needed iflag = 0 do k=1,neq c*****here oscillations are damped out - if (x1*xcorr(k) .lt. -0.5*x1**2) xcorr(k) = 0.5*xcorr(k) + if (x1*xcorr(k) < -0.5*x1**2) xcorr(k) = 0.5*xcorr(k) x1 = xatom(k) - if (dabs(xcorr(k)/xatom(k)) .gt. 0.005) iflag = 1 + if (dabs(xcorr(k)/xatom(k)) > 0.005) iflag = 1 xatom(k) = xatom(k) - xcorr(k) c*****fix element number densities which are ridiculous - if (xatom(k).le.0.0 .or. xatom(k).ge.1.001*xfic(k)) then + if (xatom(k)<=0.0 .or. xatom(k)>=1.001*xfic(k)) then iflag = 1 xatom(k) = 1.0d-2*dabs(xatom(k)+xcorr(k)) endif enddo - if (iflag .ne. 0) go to 27 + if (iflag /= 0) go to 27 c*****print out atomic and molecular partial pressures. *xamol* is the @@ -255,7 +255,7 @@ c*****number density for each neutral atom do jmol=1,nmol pmol(jmol) = dlog10(xmol(jmol,i)*tk) enddo - if (molopt .ge. 2) then + if (molopt >= 2) then pglog = dlog10(pgas(lev)) write (nf1out,1004) lev,int(t(lev)+0.001),pglog, . (patom(i),i=1,neq), (pmol(i),i=1,nmol) @@ -79,7 +79,7 @@ c call prinfo (lscreen) iatom = atom1(lim1) xab = dlog10(xabund(iatom)) + 12. ich = idint(charge(lim1) + 0.1) - if (iatom .lt. 100) then + if (iatom < 100) then write (array,1003) wave1(lim1), e(lim1,1), . dlog10(gf(lim1)), names(iatom), . ion(ich), xab, 1000.*widout(lim1) @@ -123,7 +123,7 @@ c*****(re)compute the line optical depth at line center and the C_d curve c*****compute layer where continuum optical depth > 1 do i=1,ntau - if (taulam(i) .ge. 1.) then + if (taulam(i) >= 1.) then xdepthlam1 = xdepth(i-1) + (1.-taulam(i-1))* . (xdepth(i)-xdepth(i-1))/(taulam(i)-taulam(i-1)) write (nf2out,1013) int(xdepthlam1), i @@ -133,12 +133,12 @@ c*****compute layer where continuum optical depth > 1 c compute layer where line center optical depth > 1 -10 if (taunu0(ntau) .lt. 1.) then +10 if (taunu0(ntau) < 1.) then write (nf2out,1016) go to 20 endif do i=1,ntau - if (taunu0(i) .ge. 1.) then + if (taunu0(i) >= 1.) then xdepthnu01 = xdepth(i-1) + (1.-taunu0(i-1))* . (xdepth(i)-xdepth(i-1))/(taunu0(i)-taunu0(i-1)) write (nf2out,1014) int(xdepthnu01), i @@ -149,7 +149,7 @@ c compute layer where line center optical depth > 1 c compute layer where line center plus continuum optical depth > 1 20 do i=1,ntau - if (taunu0(i)+taulam(i) .ge. 1.) then + if (taunu0(i)+taulam(i) >= 1.) then tautot1 = taulam(i-1) + taunu0(i-1) tautot2 = taulam(i) + taunu0(i) xdepthtot1 = xdepth(i-1) + (1.-tautot1)* @@ -171,7 +171,7 @@ c compute layer where line center plus continuum optical depth > 1 cdinteg = rinteg(xref,dummy1,dummy2,ntau,0.) xrefmean = xrefcdinteg/cdinteg do i=1,ntau - if (xrefmean .le. xref(i)) then + if (xrefmean <= xref(i)) then xdepthxrefmean = xdepth(i-1) + (xrefmean-xref(i-1))* . (xdepth(i)-xdepth(i-1))/(xref(i)-xref(i-1)) write (nf2out,1017) int(xdepthxrefmean), i, tauref(i), diff --git a/Ewweighted.f b/Ewweighted.f index bbbb397..bbf7aac 100755 --- a/Ewweighted.f +++ b/Ewweighted.f @@ -12,7 +12,7 @@ c****************************************************************************** c*****start the computations for a line - if (dabs(rwlgerror) .lt. 0.01) write (nf7out,1001) + if (dabs(rwlgerror) < 0.01) write (nf7out,1001) ewweighttot = 0. weighttot = 0 modcount = 0 @@ -25,7 +25,7 @@ c for the assumed abundance do mmod=1,modtot ncurvetot = nmodcurve(mmod,lim1) do icurve=3,ncurvetot-2 - if (gfmodtab(mmod,lim1,icurve) .gt. xngf) then + if (gfmodtab(mmod,lim1,icurve) > xngf) then ic = icurve - 1 pp = (xngf-gfmodtab(mmod,lim1,ic))/0.15 rw = rwmodtab(mmod,lim1,ic-1)*(-pp)*(pp-1.)*(pp-2.)/6. + @@ -41,7 +41,7 @@ c*****add this EW to the total, weighting it by flux*radius^2*relcount 10 ewweight = ew*weightmod(mmod,lim1) weighttot = weighttot + weightmod(mmod,lim1) ewweighttot = ewweighttot + ewweight - if (dabs(rwlgerror) .lt. 0.01) + if (dabs(rwlgerror) < 0.01) . write (nf7out,1005) fmodinput(mmod), . fmodoutput(mmod), radius(mmod), relcount(mmod), . fluxmod(mmod,lim1), weightmod(mmod,lim1), @@ -52,7 +52,7 @@ c*****add this EW to the total, weighting it by flux*radius^2*relcount c*****write out the mean EW ewmod(lim1) = ewweighttot/weighttot abundout(lim1) =dlog10(xabund(iatom))+deltangf+12. - if (dabs(rwlgerror) .lt. 0.01) then + if (dabs(rwlgerror) < 0.01) then write(nf7out,1006) atom1(lim1), wave1(lim1), e(lim1,1), . dlog10(gf(lim1)), 1000.*ewmod(lim1), . abundout(lim1) @@ -30,11 +30,11 @@ c*****here we use a real Fe I line because it has Barklem damping data alphabk(1) = 0.238 gambark(1) = 10.**gammabk(1) alpbark(1) = (1.-alphabk(1))/2. - if (dampingopt .eq. 0) then + if (dampingopt == 0) then damptype(1) = 'UNSLDc6' - elseif (dampingopt .eq. 1) then + elseif (dampingopt == 1) then damptype(1) = 'BKgamma' - elseif (dampingopt .eq. 2) then + elseif (dampingopt == 2) then damptype(1) = 'BLKWLc6' else damptype(1) = 'NEXTGEN' @@ -9,15 +9,15 @@ c****************************************************************************** span = end2 - end1 - if (span .lt. 0.) span = - span + if (span < 0.) span = - span spanlog = alog10(span) size = ifix(spanlog) - if (spanlog .lt. 0.) size = size - 1 + if (spanlog < 0.) size = size - 1 chop = spanlog - size - if (chop .lt. 0.31) then + if (chop < 0.31) then bigtic = 10.**(size)/2 smltic = bigtic/5 - elseif (chop .lt. 0.71) then + elseif (chop < 0.71) then bigtic = 10.**(size) smltic = bigtic/5 else @@ -10,26 +10,26 @@ c****************************************************************************** c close the files - if (nfmodel .ne. 0) close (unit=nfmodel) - if (nflines .ne. 0) close (unit=nflines) - if (nfslines .ne. 0) close (unit=nfslines) - if (nftable .ne. 0) close (unit=nftable) - if (nfobs .ne. 0) close (unit=nfobs) - if (nf1out .ne. 0) close (unit=nf1out) - if (nf2out .ne. 0) close (unit=nf2out) - if (nf3out .ne. 0) close (unit=nf3out) - if (nf4out .ne. 0) close (unit=nf4out) - if (nf5out .ne. 0) close (unit=nf5out) - if (control .ne. 'gridsyn' .and. control .ne. 'gridplo') then - if (nf6out .ne. 0) close (unit=nf6out) - if (nf7out .ne. 0) close (unit=nf7out) - if (nf8out .ne. 0) close (unit=nf8out) - if (nf9out .ne. 0) close (unit=nf9out) - if (nf10out .ne. 0) close (unit=nf10out) + if (nfmodel /= 0) close (unit=nfmodel) + if (nflines /= 0) close (unit=nflines) + if (nfslines /= 0) close (unit=nfslines) + if (nftable /= 0) close (unit=nftable) + if (nfobs /= 0) close (unit=nfobs) + if (nf1out /= 0) close (unit=nf1out) + if (nf2out /= 0) close (unit=nf2out) + if (nf3out /= 0) close (unit=nf3out) + if (nf4out /= 0) close (unit=nf4out) + if (nf5out /= 0) close (unit=nf5out) + if (control /= 'gridsyn' .and. control /= 'gridplo') then + if (nf6out /= 0) close (unit=nf6out) + if (nf7out /= 0) close (unit=nf7out) + if (nf8out /= 0) close (unit=nf8out) + if (nf9out /= 0) close (unit=nf9out) + if (nf10out /= 0) close (unit=nf10out) endif c write the closing message - if (number .eq. 0) then + if (number == 0) then istat = ivcleof (4,1) write (array,1001) istat = ivwrite (5,1,array,79) diff --git a/Gammabark.f b/Gammabark.f index 94a697e..40794d4 100755 --- a/Gammabark.f +++ b/Gammabark.f @@ -23,8 +23,8 @@ c*****on first entry to this routine, read damping data from either c 'Barklem.dat' or 'BarklemUV.dat', depending on the wavelength region c of the linelist; read Van der Waals params, if radiative damping c data are present - if (firstread .eq. 0) then - if (wave1(nlines) .gt. 3000.) then + if (firstread == 0) then + if (wave1(nlines) > 3000.) then nwant = 35 else nwant = 36 @@ -33,7 +33,7 @@ c data are present call blankstring (line) read (nwant,1001,end=10) line read (line,*) wavebk(k), idbk(k), gammabk(k), alphabk(k) - if (line(34:) .ne. ' ') then + if (line(34:) /= ' ') then read(line(34:),*) gammarad(k) else gammarad(k) = 0.0 @@ -48,20 +48,20 @@ c*****identify the Barklem list positions of the wavelength limits of c the input line list wavemin = 10000000. do j=1,nlines+nstrong - if (wave1(j) .lt. wavemin) wavemin = wave1(j) + if (wave1(j) < wavemin) wavemin = wave1(j) enddo wavemax = 0. do j=1,nlines+nstrong - if (wave1(j) .gt. wavemax) wavemax = wave1(j) + if (wave1(j) > wavemax) wavemax = wave1(j) enddo do k=1,numbark - if (wavemin-wavebk(k) .lt. 1.0) then + if (wavemin-wavebk(k) < 1.0) then nummin = k exit endif enddo do k=nummin,numbark - if (wavebk(k)-wavemax .gt. 1.0) then + if (wavebk(k)-wavemax > 1.0) then nummax = k exit endif @@ -73,19 +73,19 @@ c*****search for Barklem data gambark(j) = -1. alpbark(j) = -1. gamrad(j) = -1. - if (atom1(j) .gt. 100.) cycle + if (atom1(j) > 100.) cycle iatom10 = nint(10.*atom1(j)) do k=nummin,nummax waveerror = -(wave1(j) - wavebk(k))/wavebk(k) iii = nint(10.*idbk(k)) - if (dabs(waveerror).lt.5.0d-06 .and. - . iii .eq. iatom10) then + if (dabs(waveerror)<5.0d-06 .and. + . iii == iatom10) then gamrad(j) = gammarad(k) gambark(j) = 10.**gammabk(k) alpbark(j) = (1.-alphabk(k))/2. exit endif - if (waveerror .gt. 5.0d-06) exit + if (waveerror > 5.0d-06) exit enddo enddo @@ -16,7 +16,7 @@ c****************************************************************************** istat = ivmove(line-1,1) istat = ivcleol() istat = ivmove(line-1,1) - if (num .lt. 10) then + if (num < 10) then write (errmess,1001) num 1001 format ('(a',i1,'$)') else @@ -25,7 +25,7 @@ c****************************************************************************** endif write (*,errmess) array num = 80 - num - if (num .lt. 10) then + if (num < 10) then write (errmess,1003) num 1003 format ('(a',i1,')') else @@ -11,7 +11,7 @@ c*************************************************************************** do i=num,1,-1 - if (linechars(i:i) .ne. ' ') go to 11 + if (linechars(i:i) /= ' ') go to 11 enddo num = -1 return @@ -14,9 +14,9 @@ c****************************************************************************** xnum = -9999. 1 call getasci (nchars,line) - if (nchars .lt. 0) return + if (nchars < 0) return call number (nchars,line,xnum) - if (xnum .eq. -9999.) go to 1 + if (xnum == -9999.) go to 1 xnumsngl = sngl(xnum) return @@ -15,30 +15,30 @@ c****************************************************************************** c*****if the syntheses need to be redone: first rewind the output files, c then close/reopen line list(s), then rewrite model atmosphere output - if (choice .eq. 'n') then + if (choice == 'n') then call chabund - if (choice .eq. 'x') call pltspec (lscreen,ncall) + if (choice == 'x') call pltspec (lscreen,ncall) rewind nf1out rewind nf2out - if (nflines .ne. 0) then + if (nflines /= 0) then close (unit=nflines) open (unit=nflines,file=flines,access='sequential', . form='formatted',blank='null',status='old', . iostat=jstat,err=10) endif - if (nfslines .ne. 0) then + if (nfslines /= 0) then close (unit=nfslines) open (unit=nfslines,file=fslines,access='sequential', . form='formatted',blank='null',status='old', . iostat=jstat,err=10) endif - if (plotopt .ne. 0) then + if (plotopt /= 0) then rewind nf3out endif write (nf1out,1002) modtype - if (modprintopt .ge. 1) then - if (modtype .eq. 'begn ' .or. - . modtype .eq. 'BEGN ') write (nf1out,1003) + if (modprintopt >= 1) then + if (modtype == 'begn ' .or. + . modtype == 'BEGN ') write (nf1out,1003) write (nf1out,1102) moditle do i=1,ntau dummy1(i) = dlog10(pgas(i)) @@ -62,7 +62,7 @@ c then close/reopen line list(s), then rewrite model atmosphere output c*****now do the syntheses - if (numpecatom .eq. 0 .or. numatomsyn .eq. 0) then + if (numpecatom == 0 .or. numatomsyn == 0) then isynth = 1 isorun = 1 nlines = 0 @@ -21,7 +21,7 @@ c spectra (output), and (if desired) IRAF-style smoothed spectra (output) nchars = 19 call infile ('input ',nf2out,'formatted ',0,nchars, . f2out,lscreen) - if (plotopt .ne. 0) then + if (plotopt /= 0) then nf3out = 22 lscreen = lscreen + 2 array = 'SMOOTHED SYNTHESES OUTPUT' @@ -35,7 +35,7 @@ c spectra (output), and (if desired) IRAF-style smoothed spectra (output) call infile ('output ',nf5out,'formatted ',0,nchars, . f5out,lscreen) endif - if (plotopt .gt. 1) then + if (plotopt > 1) then nf6out = 27 lscreen = lscreen + 2 array = 'SPECTRUM COMPARISON OUTPUT' @@ -43,7 +43,7 @@ c spectra (output), and (if desired) IRAF-style smoothed spectra (output) call infile ('output ',nf6out,'formatted ',0,nchars, . f6out,lscreen) endif - if (iraf .ne. 0) then + if (iraf /= 0) then nf4out = 23 lscreen = lscreen + 2 array = 'IRAF ("rtext") OUTPUT' @@ -54,12 +54,12 @@ c spectra (output), and (if desired) IRAF-style smoothed spectra (output) c*****now plot the spectrum - if (plotopt.eq.2 .and. specfileopt.gt.0) then + if (plotopt==2 .and. specfileopt>0) then nfobs = 33 lscreen = lscreen + 2 array = 'THE OBSERVED SPECTRUM' nchars = 21 - if (specfileopt.eq.1 .or. specfileopt.eq.3) then + if (specfileopt==1 .or. specfileopt==3) then call infile ('input ',nfobs,'unformatted',2880,nchars, . fobs,lscreen) else @@ -67,7 +67,7 @@ c*****now plot the spectrum . fobs,lscreen) endif endif - if (plotopt .ne. 0) then + if (plotopt /= 0) then line = 10 ncall = 1 call pltspec (line,ncall) @@ -75,7 +75,7 @@ c*****now plot the spectrum c*****finish - if (control .ne. 'gridend') then + if (control /= 'gridend') then call finish (1) go to 1 else @@ -32,7 +32,7 @@ c spectra, and (if desired) IRAF-style smoothed spectra nchars = 20 call infile ('output ',nf2out,'formatted ',0,nchars, . f2out,lscreen) - if (plotopt .gt. 0) then + if (plotopt > 0) then nf3out = 22 lscreen = lscreen + 2 array = 'SMOOTHED SYNTHESES OUTPUT' @@ -46,7 +46,7 @@ c spectra, and (if desired) IRAF-style smoothed spectra call infile ('output ',nf5out,'formatted ',0,nchars, . f5out,lscreen) endif - if (plotopt .gt. 1) then + if (plotopt > 1) then nf6out = 27 lscreen = lscreen + 2 array = 'SPECTRUM COMPARISON OUTPUT' @@ -54,7 +54,7 @@ c spectra, and (if desired) IRAF-style smoothed spectra call infile ('output ',nf6out,'formatted ',0,nchars, . f6out,lscreen) endif - if (iraf .ne. 0) then + if (iraf /= 0) then nf4out = 23 lscreen = lscreen + 2 array = 'IRAF ("rtext") OUTPUT' @@ -81,7 +81,7 @@ c*****open the line list file and the strong line list file nchars = 13 call infile ('input ',nflines,'formatted ',0,nchars, . flines,lscreen) - if (dostrong .gt. 0) then + if (dostrong > 0) then nfslines = 32 lscreen = lscreen + 2 array = 'THE STRONG LINE LIST' @@ -92,7 +92,7 @@ c*****open the line list file and the strong line list file c*****do the syntheses - if (numpecatom .eq. 0 .or. numatomsyn .eq. 0) then + if (numpecatom == 0 .or. numatomsyn == 0) then isorun = 1 isynth = 1 nlines = 0 @@ -119,12 +119,12 @@ c*****do the syntheses c*****now plot the spectrum - if (plotopt.eq.2 .and. specfileopt.gt.0) then + if (plotopt==2 .and. specfileopt>0) then nfobs = 33 lscreen = lscreen + 2 array = 'THE OBSERVED SPECTRUM' nchars = 21 - if (plotopt.eq.1 .or. specfileopt.eq.3) then + if (plotopt==1 .or. specfileopt==3) then call infile ('input ',nfobs,'unformatted',2880,nchars, . fobs,lscreen) else @@ -132,14 +132,14 @@ c*****now plot the spectrum . fobs,lscreen) endif endif - if (plotopt .ne. 0) then + if (plotopt /= 0) then ncall = 1 call pltspec (lscreen,ncall) endif c*****finish - if (control .ne. 'gridend') then + if (control /= 'gridend') then call finish (1) go to 1 else @@ -13,19 +13,19 @@ c****************************************************************************** c decide on the file status desired jstat = 0 - if (type .eq. 'input ') then + if (type == 'input ') then kstat = 'old ' - elseif (type .eq. 'output ') then + elseif (type == 'output ') then kstat = 'new ' - elseif (type .eq. 'overout') then + elseif (type == 'overout') then kstat = 'unknown' endif c write out the appropriate message about this file 5 nchars = charcount - if (fname .eq. 'optional_output_file') then + if (fname == 'optional_output_file') then return - elseif (fname .eq. 'no_filename_given') then + elseif (fname == 'no_filename_given') then array(charcount+1:charcount+24) ='; what is the filename? ' charcount = charcount + 24 call getasci (charcount,line) @@ -35,11 +35,11 @@ c write out the appropriate message about this file array(charcount+25:79) = fname charcount = 79 call putasci (charcount,line) - if (type .ne. 'input ') kstat = 'unknown' + if (type /= 'input ') kstat = 'unknown' endif c open the file specified by the user, earlier or now -6 if (mode .eq. 'formatted ') then +6 if (mode == 'formatted ') then open (unit=iunit,file=fname,access='sequential', . form=mode,blank='null',status=kstat, . iostat=jstat,err=10) @@ -55,8 +55,8 @@ c open the file specified by the user, earlier or now c here are the file reading error messages; c if an expected file is not found, 118 is the error code for SunOS, 1018 c is for Solaris, and 2 is for Redhat Linux operating systems. -10 if (jstat .eq. 118 .or. jstat .eq. 1018 .or. - . jstat .eq. 2) then +10 if (jstat == 118 .or. jstat == 1018 .or. + . jstat == 2) then write (errmess,1001) jstat istat = ivwrite (line+2,3,errmess,44) fname = 'no_filename_given' @@ -64,12 +64,12 @@ c is for Solaris, and 2 is for Redhat Linux operating systems. go to 5 c if a file is in danger of being over-written, 117 is the error code for c SunOS, 1017 is for Solaris, and 128 is for Redhat Linux operating systems. - elseif (jstat .eq. 117 .or. jstat .eq. 1017 .or. - . jstat .eq. 128) then + elseif (jstat == 117 .or. jstat == 1017 .or. + . jstat == 128) then write (errmess,1002) jstat istat = ivwrite (line+2,3,errmess,41) read (*,*) yesno - if (yesno .eq. 'y') then + if (yesno == 'y') then kstat = 'unknown' go to 6 else @@ -17,38 +17,38 @@ c****************************************************************************** integer n2 - if (num .eq. 2) go to 4 - if (num .eq. 6) go to 340 + if (num == 2) go to 4 + if (num == 6) go to 340 n1marker = 1 n2 = 0 c*****decide if certain element abundances need to be modified. - if (numpecatom .gt. 0) then + if (numpecatom > 0) then do iatom=3,95 xabund(iatom) = 10.**pecabund(iatom,isynth)* . 10.**abfactor(isynth)*xabu(iatom) enddo endif - if (num .ne. 5) then + if (num /= 5) then write (nf1out,1004) xmetals = abscale + abfactor(isynth) - if (ninetynineflag .eq. 1) then + if (ninetynineflag == 1) then write (nf1out,1005) xmetals - if (nf2out .gt. 0) write (nf2out,1005) xmetals + if (nf2out > 0) write (nf2out,1005) xmetals else - if (nf2out .gt. 0) write (nf2out,1006) abscale + if (nf2out > 0) write (nf2out,1006) abscale endif do j=1,93 - if (pec(j) .gt. 0 ) then + if (pec(j) > 0 ) then dummy1(j) = dlog10(xabund(j)) + 12.0 - if (dummy1(j) .le. -10.) then + if (dummy1(j) <= -10.) then write (nf1out,1008) names(j),dummy1(j) - if (nf2out .gt. 0) + if (nf2out > 0) . write (nf2out,1008) names(j),dummy1(j) else write (nf1out,1007) names(j),dummy1(j) - if (nf2out .gt. 0) + if (nf2out > 0) . write (nf2out,1007) names(j),dummy1(j) endif endif @@ -57,12 +57,12 @@ c*****decide if certain element abundances need to be modified. c*****output information about the isotopic ratios - if (numiso .gt. 0) then + if (numiso > 0) then write (nf1out,1014) do i=1,numiso iiso = isotope(i) write (nf1out,1015) iiso, isotope(i), isoabund(i,isorun) - if (nf2out .gt. 0) write (nf2out,1015) + if (nf2out > 0) write (nf2out,1015) . iiso, isotope(i), isoabund(i,isorun) enddo endif @@ -74,7 +74,7 @@ c if 'printstrong' gt 0 then the strong lines have c been printed printstrong = -1 - if (num .ne. 4) then + if (num /= 4) then rewind nflines wave = start read (nflines,1001) linitle @@ -83,10 +83,10 @@ c been printed c*****read in the strong lines if needed 302 nstrong = 0 - if (dostrong .gt. 0 ) then + if (dostrong > 0 ) then rewind nfslines do j=1,41 - if (linfileopt .eq. 0) then + if (linfileopt == 0) then read (nfslines,1002,end=340) swave1(j),satom1(j),se(j), . sgf(j),sdampnum(j),sd0(j),swidth(j) else @@ -97,12 +97,12 @@ c*****read in the strong lines if needed iatom = satom1(j) scharge(j) = 1.0 + dble(int(10.0*(satom1(j) - iatom) . +0.0001)) - if (scharge(j) .gt. 3.) then + if (scharge(j) > 3.) then write (*,1003) swave1(i), satom1(i) stop endif enddo - if (nstrong .gt. 40) then + if (nstrong > 40) then write(*,*) 'STRONG LINE LIST HAS MORE THAN 40 LINES. THIS' write(*,*) 'IS NOT ALLOWED. I QUIT!' stop @@ -111,7 +111,7 @@ c*****read in the strong lines if needed 340 nlines = 2500 - nstrong j = 1 -333 if (linfileopt .eq. 0) then +333 if (linfileopt == 0) then read (nflines,1002,end=311) wave1(j),atom1(j),e(j,1),gf(j), . dampnum(j),d0(j),width(j) else @@ -120,26 +120,26 @@ c*****read in the strong lines if needed endif iatom = atom1(j) charge(j) = 1.0 + dble(int(10.0*(atom1(j) - iatom)+0.0001)) - if (charge(j) .gt. 3.) then + if (charge(j) > 3.) then write (*,1003) wave1(j), atom1(j) stop endif - if (width(j) .lt. 0.) then - if (control .eq. 'blends ') then + if (width(j) < 0.) then + if (control == 'blends ') then write (*,*) 'BLENDS cannot have negative EWs! I QUIT!' stop else go to 333 endif endif - if (iunits .eq. 1) wave1(j) = 1.d+4*wave1(j) + if (iunits == 1) wave1(j) = 1.d+4*wave1(j) j = j + 1 - if (j .le. nlines) go to 333 + if (j <= nlines) go to 333 311 nlines = j - 1 c*****append the strong lines here if necessary - if (dostrong .gt. 0) then + if (dostrong > 0) then do k=1,nstrong wave1(nlines+k) = swave1(k) atom1(nlines+k) = satom1(k) @@ -155,7 +155,7 @@ c*****append the strong lines here if necessary c*****here groups of lines for blended features are defined do j=1,nlines+nstrong - if (wave1(j) .lt. 0.) then + if (wave1(j) < 0.) then group(j) = 1 wave1(j) = dabs(wave1(j)) width(j) = width(j-1) @@ -167,7 +167,7 @@ c*****here groups of lines for blended features are defined c*****here excitation potentials are changed from cm^-1 to eV, if needed do j=1,nlines+nstrong - if (e(j,1) .gt. 50.) then + if (e(j,1) > 50.) then do jj=1,nlines+nstrong e(jj,1) = 1.2389e-4*e(jj,1) enddo @@ -178,7 +178,7 @@ c*****here excitation potentials are changed from cm^-1 to eV, if needed c*****here log(gf) values are turned into gf values, if needed do j=1,nlines+nstrong - if (gfstyle.eq.0 .or. gf(j) .lt. 0) then + if (gfstyle==0 .or. gf(j) < 0) then do jj=1,nlines+nstrong gf(jj) = 10.**gf(jj) enddo @@ -190,7 +190,7 @@ c*****here log(gf) values are turned into gf values, if needed c*****turn log(RW) values and EW values in mA into EW values in A. Stuff c duplicate EW values of the first line of a blend into all blend members. do j=1,nlines+nstrong - if (width(j) .lt. 0.) then + if (width(j) < 0.) then width(j) = 10.**width(j)*wave1(j) else width(j) = width(j)/1000. @@ -208,13 +208,13 @@ c and a different one for atomic lines c*****here are the calculations specific to molecular lines - if (iatom .ge. 100) then + if (iatom >= 100) then call sunder (atom1(j),ia,ib) - if (ia .gt. ib) then + if (ia > ib) then write (*,1010) ia,ib stop endif - if (atom10-int(atom10) .le. 0.0) then + if (atom10-int(atom10) <= 0.0) then amass(j) = xam(ia) + xam(ib) mas1 = xam(ia) + 0.0000001 mas2 = xam(ib) + 0.0000001 @@ -223,8 +223,8 @@ c*****here are the calculations specific to molecular lines mas1 = jat100 - 100*int(atom10) jat10000 = int(10000.*(atom10+0.00001)) mas2 = jat10000 - 100*jat100 - if (mas1.gt.mas2 .or. mas1.le.0.0 .or. - . mas2.le.0.0) then + if (mas1>mas2 .or. mas1<=0.0 .or. + . mas2<=0.0) then write (*,1011) mas1, mas2 stop endif @@ -232,9 +232,9 @@ c*****here are the calculations specific to molecular lines endif c*****use an internal dissociation energy for molecules if the user c does not read one in - if (d0(j) .eq. 0.) then + if (d0(j) == 0.) then do k=1,110 - if (int(datmol(1,k)+0.01) .eq. + if (int(datmol(1,k)+0.01) == . int(atom1(j)+0.01)) then d0(j) = datmol(2,k) go to 390 @@ -251,7 +251,7 @@ c does not read one in c*****here are the calculations specific to atomic lines else - if (atom10-int(atom10) .le. 0.0) then + if (atom10-int(atom10) <= 0.0) then amass(j) = xam(iatom) else atom10 = atom10 + 0.00001 @@ -266,8 +266,8 @@ c*****here are the calculations specific to atomic lines c*****quit the routine normally - if (nlines+nstrong .lt. 2500) then - if (sstop .gt. wave1(nlines)+10.) sstop = wave1(nlines)+10. + if (nlines+nstrong < 2500) then + if (sstop > wave1(nlines)+10.) sstop = wave1(nlines)+10. endif lim1line = 1 return @@ -26,7 +26,7 @@ c*****Read in the key word to define the model type rewind nfmodel read (nfmodel,2001) modtype write (nf1out,1010) modtype - if (modtype .eq. 'begn ' .or. modtype .eq. 'BEGN ') + if (modtype == 'begn ' .or. modtype == 'BEGN ') . write (nf1out,1011) @@ -38,7 +38,7 @@ c*****Read the number of depth points read (nfmodel,2002) list list2 = list(11:) read (list2,*) ntau - if (ntau .gt. 100) then + if (ntau > 100) then write (array,1012) call prinfo (10) stop @@ -50,7 +50,7 @@ c MARCS code. This modtype is called "NEWMARCS". On each line c the numbers are: c tau(5000), t, pe, pgas, rho, model microtrubulent velocity, c and mean opacity (cm^2/gm) at the reference wavelength (5000A). - if (modtype .eq. 'NEWMARCS ') then + if (modtype == 'NEWMARCS ') then read (nfmodel,*) wavref do i=1,ntau read (nfmodel,*) tauref(i),t(i),ne(i),pgas(i),rho(i), @@ -62,7 +62,7 @@ c the numbers are: c layer number (not needed), log{tau(Rosseland)} (not needed), c log{tau(5000)}, depth, t, pe, pgas, prad (not read in) and c pturb (not read in) - elseif (modtype .eq. 'WEBMARCS') then + elseif (modtype == 'WEBMARCS') then read (nfmodel,*) wavref do i=1,ntau read (nfmodel,*) k, dummy1(k), tauref(i), dummy2(k), t(i), @@ -73,7 +73,7 @@ c MARCS code. This modtype is called "WEB2MARC". On each line c the numbers are: c atmospheric layer number (not needed), log{tau(5000)}, t, c log(Pe), log(Pgas), rhox - elseif (modtype .eq. 'WEB2MARC') then + elseif (modtype == 'WEB2MARC') then read (nfmodel,*) wavref do i=1,ntau read (nfmodel,*) k,tauref(i),t(i),ne(i),pgas(i),rhox(i) @@ -82,7 +82,7 @@ c OR: Read in a model from the output of the ATLAS code. This c modtype is called "KURUCZ". On each line the numbers are: c rhox, t, pgas, ne, and Rosseland mean opacity (cm^2/gm), and c two numbers not used by MOOG. - elseif (modtype .eq. 'KURUCZ ') then + elseif (modtype == 'KURUCZ ') then do i=1,ntau read (nfmodel,*) rhox(i),t(i),pgas(i),ne(i),kaprefmass(i) enddo @@ -92,7 +92,7 @@ c specific wavelength that is read in before the model. MOOG will c need to generate the opacities internally.On each line the numbers c are: tau, t, pgas, pe, density, mean molecular weight, two numbers c not used by MOOG, and Rosseland mean opacity (cm^2/gm). - elseif (modtype .eq. 'NEXTGEN ') then + elseif (modtype == 'NEXTGEN ') then read (nfmodel,*) wavref do i=1,ntau read (nfmodel,*) tauref(i),t(i),pgas(i),ne(i), rho(i), @@ -101,7 +101,7 @@ c not used by MOOG, and Rosseland mean opacity (cm^2/gm). c OR: Read in a model from the output of the MARCS code. This modtype c type is called "BEGN". On each line the numbers are: c tauross, t, log(pg), log(pe), mol weight, and kappaross. - elseif (modtype .eq. 'BEGN ') then + elseif (modtype == 'BEGN ') then do i=1,ntau read (nfmodel,*) tauref(i),t(i),pgas(i),ne(i), . molweight(i), kaprefmass(i) @@ -109,14 +109,14 @@ c tauross, t, log(pg), log(pe), mol weight, and kappaross. c OR: Read in a model generated from ATLAS, but without accompanying c opacities. MOOG will need to generate the opacities internally, c using a reference wavelength that it reads in before the model. - elseif (modtype .eq. 'KURTYPE') then + elseif (modtype == 'KURTYPE') then read (nfmodel,*) wavref do i=1,ntau read (nfmodel,*) rhox(i),t(i),pgas(i),ne(i) enddo c OR: Read in a model generated from ATLAS, with output generated c in Padova. The columns are in somewhat different order than normal - elseif (modtype .eq. 'KUR-PADOVA') then + elseif (modtype == 'KUR-PADOVA') then read (nfmodel,*) wavref do i=1,ntau read (nfmodel,*) tauref(i),t(i),kaprefmass(i), @@ -125,7 +125,7 @@ c in Padova. The columns are in somewhat different order than normal c OR: Read in a generic model that has a tau scale at a specific c wavelength that is read in before the model. c MOOG will need to generate the opacities internally. - elseif (modtype .eq. 'GENERIC ') then + elseif (modtype == 'GENERIC ') then read (nfmodel,*) wavref do i=1,ntau read (nfmodel,*) tauref(i),t(i),pgas(i),ne(i) @@ -146,7 +146,7 @@ c*****Compute other convenient forms of the temperatures c*****Convert from logarithmic Pgas scales, if needed - if (pgas(ntau)/pgas(1) .lt. 10.) then + if (pgas(ntau)/pgas(1) < 10.) then do i=1,ntau pgas(i) = 10.0**pgas(i) enddo @@ -154,7 +154,7 @@ c*****Convert from logarithmic Pgas scales, if needed c*****Convert from logarithmic Ne scales, if needed - if(ne(ntau)/ne(1) .lt. 20.) then + if(ne(ntau)/ne(1) < 20.) then do i=1,ntau ne(i) = 10.0**ne(i) enddo @@ -162,7 +162,7 @@ c*****Convert from logarithmic Ne scales, if needed c*****Convert from Pe to Ne, if needed - if(ne(ntau) .lt. 1.0e7) then + if(ne(ntau) < 1.0e7) then do i=1,ntau ne(i) = ne(i)/1.38054d-16/t(i) enddo @@ -180,14 +180,14 @@ c*****Read the microturbulence (either a single value to apply to c all layers, or a value for each of the ntau layers). c Conversion to cm/sec from km/sec is done if needed read (nfmodel,2003) (vturb(i),i=1,6) - if (vturb(2) .ne. 0.) then + if (vturb(2) /= 0.) then read (nfmodel,2003) (vturb(i),i=7,ntau) else do i=2,ntau vturb(i) = vturb(1) enddo endif - if (vturb(1) .lt. 100.) then + if (vturb(1) < 100.) then write (moditle(55:62),1008) vturb(1) do i=1,ntau vturb(i) = 1.0e5*vturb(i) @@ -204,7 +204,7 @@ c*****solar ones contained in array xsolar. list2 = list(11:) read (list2,*) natoms,abscale write (moditle(63:73),1009) abscale - if(natoms .ne. 0) + if(natoms /= 0) . read (nfmodel,*) (element(i),logepsilon(i),i=1,natoms) xhyd = 10.0**xsolar(1) xabund(1) = 1.0 @@ -213,7 +213,7 @@ c*****solar ones contained in array xsolar. xabund(i) = 10.0**(xsolar(i)+abscale)/xhyd xabu(i) = xabund(i) enddo - if (natoms .ne. 0) then + if (natoms /= 0) then do i=1,natoms xabund(idint(element(i))) = 10.0**logepsilon(i)/xhyd xabu(idint(element(i))) = 10.0**logepsilon(i)/xhyd @@ -230,17 +230,17 @@ c in this approximation (maybe make more general some day?) enddo wtmol = wtnum/(xam(1)*wtden) nomolweight = 0 - if (modtype .eq. 'BEGN ' .or. modtype .eq. 'NEXTGEN ') then + if (modtype == 'BEGN ' .or. modtype == 'NEXTGEN ') then nomolweight = 1 endif - if (nomolweight .ne. 1) then + if (nomolweight /= 1) then do i=1,ntau molweight(i) = wtmol enddo endif c*****Compute the density - if (modtype .ne. 'NEXTGEN ') then + if (modtype /= 'NEXTGEN ') then do i=1,ntau rho(i) = pgas(i)*molweight(i)*1.6606d-24/(1.38054d-16*t(i)) enddo @@ -273,16 +273,16 @@ c saved. c*****Set up the default molecule list - if (molset .eq. 0) then + if (molset == 0) then nmol = 30 else nmol = 59 endif - if (molset .eq. 0) then + if (molset == 0) then do i=1,110 amol(i) = smallmollist(i) enddo - elseif (molset .eq. 1) then + elseif (molset == 1) then do i=1,110 amol(i) = largemollist(i) enddo @@ -298,15 +298,15 @@ c molecular equilibrium if needed. read (nfmodel,2002,end=101) list list2 = list(11:) read (list2,*) moremol - if (moremol .ne. 0) then + if (moremol /= 0) then read (nfmodel,*) (bmol(i),i=1,moremol) append = 1 do k=1,moremol do l=1,nmol - if (nint(bmol(k)) .eq. nint(amol(l))) + if (nint(bmol(k)) == nint(amol(l))) . append = 0 enddo - if (append .eq. 1) then + if (append == 1) then nmol = nmol + 1 amol(nmol) = bmol(k) endif @@ -320,13 +320,13 @@ c*****do the general molecular equilibrium c*****SPECIAL NEEDS: for NEWMARCS models, to convert kaprefs to our units - if (modtype .eq. 'NEWMARCS ') then + if (modtype == 'NEWMARCS ') then do i=1,ntau kapref(i) = kaprefmass(i)*rho(i) enddo c SPECIAL NEEDS: for KURUCZ models, to create the optical depth array, c and to convert kaprefs to our units - elseif (modtype .eq. 'KURUCZ ') then + elseif (modtype == 'KURUCZ ') then first = rhox(1)*kaprefmass(1) tottau = rinteg(rhox,kaprefmass,tauref,ntau,first) tauref(1) = first @@ -337,18 +337,18 @@ c and to convert kaprefs to our units kapref(i) = kaprefmass(i)*rho(i) enddo c SPECIAL NEEDS: for NEXTGEN models, to convert kaprefs to our units - elseif (modtype .eq. 'NEXTGEN ') then + elseif (modtype == 'NEXTGEN ') then do i=1,ntau kapref(i) = kaprefmass(i)*rho(i) enddo c SPECIAL NEEDS: for BEGN models, to convert kaprefs to our units - elseif (modtype .eq. 'BEGN ') then + elseif (modtype == 'BEGN ') then do i=1,ntau kapref(i) = kaprefmass(i)*rho(i) enddo c SPECIAL NEEDS: for KURTYPE models, to create internal kaprefs, c and to compute taurefs from the kaprefs converted to mass units - elseif (modtype .eq. 'KURTYPE ') then + elseif (modtype == 'KURTYPE ') then call opacit (1,wavref) do i=1,ntau kaprefmass(i) = kapref(i)/rho(i) @@ -360,21 +360,21 @@ c and to compute taurefs from the kaprefs converted to mass units tauref(i) = tauref(i-1) + tauref(i) enddo c SPECIAL NEEDS: for NEWMARCS models, to convert kaprefs to our units - elseif (modtype .eq. 'KUR-PADOVA') then + elseif (modtype == 'KUR-PADOVA') then do i=1,ntau kapref(i) = kaprefmass(i)*rho(i) enddo c SPECIAL NEEDS: for generic models, to create internal kaprefs, - elseif (modtype .eq. 'GENERIC ' .or. - . modtype .eq. 'WEBMARCS ' .or. - . modtype .eq. 'WEB2MARC ') then + elseif (modtype == 'GENERIC ' .or. + . modtype == 'WEBMARCS ' .or. + . modtype == 'WEB2MARC ') then call opacit (1,wavref) endif c*****Convert from logarithmic optical depth scales, or vice versa. c xref will contain the log of the tauref - if(tauref(1) .lt. 0.) then + if(tauref(1) < 0.) then do i=1,ntau xref(i) = tauref(i) tauref(i) = 10.0**xref(i) @@ -387,7 +387,7 @@ c xref will contain the log of the tauref c*****Write information to output files - if (modprintopt .lt. 1) return + if (modprintopt < 1) return write (nf1out,1002) moditle do i=1,ntau dummy1(i) = dlog10(pgas(i)) @@ -29,18 +29,18 @@ c****************************************************************************** c here begins the inversion do n=1,order - if (n .eq. order) go to 244 + if (n == order) go to 244 g = a(n,n) n1 = n+1 ng = n do m=n1,order f = a(m,n) - if (f .gt. g) then + if (f > g) then g = f ng = m endif enddo - if (ng .eq. n) go to 244 + if (ng == n) go to 244 do k=1,order d = i(n,k) e = i(ng,k) @@ -52,7 +52,7 @@ c here begins the inversion a(n,k) = g enddo 244 g = a(n,n) - if (g .eq. 0.0) then + if (g == 0.0) then write (nf1out,1001) return endif @@ -60,7 +60,7 @@ c here begins the inversion a(n,k) = a(n,k)/g i(n,k) = i(n,k)/g enddo - if (n .eq. order) go to 27 + if (n == order) go to 27 do k=n1,order f = -a(k,n) do j=1,order @@ -26,8 +26,8 @@ c****************************************************************************** 1.0435694088/ data xsave /0./ c - if (x.ge.100.) goto 800 - if (x.le.-100.) goto 800 + if (x>=100.) goto 800 + if (x<=-100.) goto 800 xu=x if(xu)603,602,603 602 rex=1.d+00 diff --git a/Lineabund.f b/Lineabund.f index 9c709ea..c7d807d 100755 --- a/Lineabund.f +++ b/Lineabund.f @@ -23,7 +23,7 @@ c*****find the c-o-g gf that matches the observed RW rwlgobs = dlog10(width(lim1)/wave1(lim1)) gfobs = gftab(ntabtot) do i=2,ntabtot - if (rwtab(i) .gt. rwlgobs) then + if (rwtab(i) > rwlgobs) then gfobs = gftab(i-1) + (gftab(i)-gftab(i-1))* . (rwlgobs-rwtab(i-1))/(rwtab(i)-rwtab(i-1)) go to 15 @@ -37,7 +37,7 @@ c*****the computed RW rwlgcal = dlog10(w(ncurve)/wave1(lim1)) gfcal = gftab(ntabtot) do i=2,ntabtot - if (rwtab(i) .gt. rwlgcal) then + if (rwtab(i) > rwlgcal) then gfcal = gftab(i-1) + (gftab(i)-gftab(i-1))* . (rwlgcal-rwtab(i-1))/(rwtab(i)-rwtab(i-1)) go to 20 @@ -52,9 +52,9 @@ c*****square root of the proposed gf shift. 20 error = (w(ncurve)-width(lim1))/width(lim1) ratio = 10.**(gfobs-gfcal) ncurve = ncurve + 1 - if (dabs(error) .ge. 0.0015 .and. ncurve .lt. 20) then + if (dabs(error) >= 0.0015 .and. ncurve < 20) then rwlcomp = dlog10(w(ncurve-1)/wave1(lim1)) - if (rwlcomp .gt. -4.7) then + if (rwlcomp > -4.7) then gf1(ncurve) = gf1(ncurve-1)*dsqrt(ratio) do i=1,ntau kapnu0(lim1,i) = kapnu0(lim1,i)*dsqrt(ratio) @@ -71,7 +71,7 @@ c*****square root of the proposed gf shift. c*****if the observed and computed RWs are close, do a final gf adjustment c*****and finish with one more line recomputation - if (ncurve .eq. 30) then + if (ncurve == 30) then write (nf1out,1001) write (nf2out,1001) endif @@ -84,7 +84,7 @@ c*****and finish with one more line recomputation wid1comp(lim1) = w(1) diff = dlog10(gf1(ncurve)/gf(lim1)) abundout(lim1) = abundin + diff - if (ncurve .ne. 1) then + if (ncurve /= 1) then write (nf1out,1002) ncurve endif @@ -30,7 +30,7 @@ c*****here the line data are output to "standard_out"; all relevant c drivers use this c if you don't want any line output, linprintopt=0 will exit the routine -1 if (linprintopt .lt. 1) return +1 if (linprintopt < 1) return c if you want standard output, linprintopt=1 is chosen c linprintopt>=2 outputs ionization potentials, charges, masses, @@ -38,14 +38,14 @@ c reduced masses for molecules, c linprintopt>=3 outputs partition functions c lineprintop =4 outputs line-center opacities write (nf1out,1001) nlines - if (linprintopt .ge. 2) write (nf1out,1002) + if (linprintopt >= 2) write (nf1out,1002) do j=1,nlines ich = idint(charge(j) + 0.1) iatom = idint(atom1(j)) loggf = dlog10(gf(j)) logstrength = dlog10(strength(j)) - if (iatom .lt. 100) then - if (iunits .eq. 1) then + if (iatom < 100) then + if (iunits == 1) then write (nf1out,1003) j, 1.d-4*wave1(j), names(iatom), . ion(ich), atom1(j), e(j,1), loggf, damptype(j), . logstrength, 1000.*width(j) @@ -54,23 +54,23 @@ c lineprintop =4 outputs line-center opacities . ion(ich), atom1(j), e(j,1), loggf, damptype(j), . logstrength, 1000.*width(j) endif - if (linprintopt .ge. 2) write (nf1out,1005) + if (linprintopt >= 2) write (nf1out,1005) . (chi(j,k),k=1,3), charge(j), amass(j), rdmass(j) - elseif (iatom .lt. 10000) then + elseif (iatom < 10000) then call sunder (atom1(j),i1,i2) - if (i1 .eq. 1) then + if (i1 == 1) then l = i1 i1 = i2 i2 = l endif leftovr = idint(10000.*(atom1(j)-iatom)+0.1) - if (i1 .lt. 10) then + if (i1 < 10) then read (names(i1),1006) name write (molname,1007) name,names(i2),leftovr else write (molname,1008) names(i1),names(i2),leftovr endif - if (iunits .eq. 1) then + if (iunits == 1) then write (nf1out,1009) j, 1.d-4*wave1(j), molname, . atom1(j), e(j,1), loggf, damptype(j), . logstrength, 1000.*width(j) @@ -79,20 +79,20 @@ c lineprintop =4 outputs line-center opacities . atom1(j), e(j,1), loggf, damptype(j), . logstrength, 1000.*width(j) endif - if (linprintopt .ge. 2) + if (linprintopt >= 2) . write (nf1out,1005) . d0(j), (chi(j,k),k=1,2), charge(j), amass(j), . rdmass(j) - elseif (iatom .lt. 1000000) then + elseif (iatom < 1000000) then call sunder (atom1(j),i1,i2) xia = dble(i2) call sunder (xia,i2,i3) - if (iatom .eq. 10108) then + if (iatom == 10108) then molname = 'H_2O ' else molname = 'CO_2 ' endif - if (iunits .eq. 1) then + if (iunits == 1) then write (nf1out,1009) j, 1.d-4*wave1(j), molname, . atom1(j), e(j,1), loggf, damptype(j), . logstrength, 1000.*width(j) @@ -101,24 +101,24 @@ c lineprintop =4 outputs line-center opacities . atom1(j), e(j,1), loggf, damptype(j), . logstrength, 1000.*width(j) endif - if (linprintopt .ge. 2) + if (linprintopt >= 2) . write (nf1out,1005) . d0(j), (chi(j,k),k=1,2), charge(j), amass(j), . rdmass(j) endif enddo - if (start.ne.0.0 .or. sstop.ne.0.0) then - if (iunits .eq. 1) then + if (start/=0.0 .or. sstop/=0.0) then + if (iunits == 1) then write (nf1out,1011) oldstart,oldstop,oldstep,olddelta else write (nf1out,1012) start,sstop,step,delta endif - if (rwlow .ne. 0.) write (nf1out,1013) rwlow, rwhigh, rwstep + if (rwlow /= 0.) write (nf1out,1013) rwlow, rwhigh, rwstep endif - if (linprintopt .ge. 3) then + if (linprintopt >= 3) then write (nf1out,1014) do j=1,95 - if (elem(j) .ne. 0.) then + if (elem(j) /= 0.) then iatom = int(elem(j)) write (nf1out,1015) iatom, names(iatom), xam(j), . xchi1(j), xchi2(j), xchi3(j) @@ -128,7 +128,7 @@ c lineprintop =4 outputs line-center opacities endif enddo endif - if (linprintopt .ge. 4) then + if (linprintopt >= 4) then write (nf1out,1001) do j=1,nlines write (nf1out,1002) j,(kapnu0(j,i),i=1,ntau) @@ -145,8 +145,8 @@ c molecular line can possibly be in this category iatom = idint(atom1(j)) loggf = dlog10(gf(j)) logstrength = dlog10(strength(j)) - if (iatom .lt. 100) then - if (iunits .eq. 1) then + if (iatom < 100) then + if (iunits == 1) then write (nf1out,1003) j-nlines,1.d-4*wave1(j),names(iatom), . ion(ich), atom1(j), e(j,1), loggf, . damptype(j), logstrength @@ -166,18 +166,18 @@ c molecular line can possibly be in this category c*****results of force-fitting EW to yield abundances are output here c look here also for the calls to the trend line calculations -3 if (ifresh .eq.0) then +3 if (ifresh ==0) then write (nf2out,3001) linitle,moditle ifresh = 1 endif - if (cogatom .eq. 0.) then + if (cogatom == 0.) then iatom = iabatom else iatom = idint(cogatom) endif xab = dlog10(xabund(iatom)) + 12. ich = idint(charge(lim1obs) + 0.1) - if (atom1(lim1obs) .lt. 100.) then + if (atom1(lim1obs) < 100.) then write (array,3002) names(iatom), ion(ich) ,xab line = 1 call prinfo (line) @@ -189,13 +189,13 @@ c look here also for the calls to the trend line calculations write (nf2out,3003) else call sunder (atom1(lim1obs),ia,ib) - if (ia .eq. 1) then + if (ia == 1) then l = ia ia = ib ib = l endif leftovr = idint(10000.*(atom1(lim1obs)-iatom)+0.1) - if (ia .lt. 10) then + if (ia < 10) then read (names(ia),1006) name write (molname,1007) name,names(ib) else @@ -216,7 +216,7 @@ c look here also for the calls to the trend line calculations write (nf2out,3006) endif do l=lim1obs,lim2obs - if (abundout(l) .ne. 999.99) then + if (abundout(l) /= 999.99) then diff = abundout(l) - average else diff = 999.99 @@ -226,7 +226,7 @@ c look here also for the calls to the trend line calculations loggf = dlog10(gf(l)) write (array,3007) wave1(l), atom1(l), e(l,1), loggf, . ew, rw, abundout(l), diff - if (errmess(1:9) .ne. 'stopinfo!') then + if (errmess(1:9) /= 'stopinfo!') then line = line + 1 call prinfo (line) endif @@ -235,48 +235,48 @@ c look here also for the calls to the trend line calculations enddo write (array,3008) average, deviate, kount line = line + 1 - if (errmess(1:9) .ne. 'stopinfo!') call prinfo (line) + if (errmess(1:9) /= 'stopinfo!') call prinfo (line) write (nf2out,3008) average, deviate, kount - if (kount .gt. 2 .and. deltaep .gt. 1.5) then + if (kount > 2 .and. deltaep > 1.5) then write (array,3009) xxm1, xxb1, xxr1 - if (errmess(1:9) .ne. 'stopinfo!') then + if (errmess(1:9) /= 'stopinfo!') then line = line + 1 call prinfo (line) endif write (nf2out,3009) xxm1, xxb1, xxr1 else write (array,*) 'No statistics done for E.P. trends' - if (errmess(1:9) .ne. 'stopinfo!') then + if (errmess(1:9) /= 'stopinfo!') then line = line + 1 call prinfo (line) endif write (nf2out,*) 'No statistics done for E.P. trends' endif - if (kount .gt. 2 .and. deltarw .gt. 0.5) then + if (kount > 2 .and. deltarw > 0.5) then write (array,3010) xxm2, xxb2, xxr2 - if (errmess(1:9) .ne. 'stopinfo!') then + if (errmess(1:9) /= 'stopinfo!') then line = line + 1 call prinfo (line) endif write (nf2out,3010) xxm2, xxb2, xxr2 else write (array,*) 'No statistics done for R.W. trends' - if (errmess(1:9) .ne. 'stopinfo!') then + if (errmess(1:9) /= 'stopinfo!') then line = line + 1 call prinfo (line) endif write (nf2out,*) 'No statistics done for R.W. trends' endif - if (kount .gt. 2 .and. deltawv .gt. 500.) then + if (kount > 2 .and. deltawv > 500.) then write (array,3011) xxm3, xxb3, xxr3 - if (errmess(1:9) .ne. 'stopinfo!') then + if (errmess(1:9) /= 'stopinfo!') then line = line + 1 call prinfo (line) endif write (nf2out,3011) xxm3, xxb3, xxr3 else write (array,*) 'No statistics done for wavelength trends' - if (errmess(1:9) .ne. 'stopinfo!') then + if (errmess(1:9) /= 'stopinfo!') then line = line + 1 call prinfo (line) endif @@ -14,7 +14,7 @@ c****************************************************************************** lineflag = 0 c*****for single-line computations, the line rage is the whole line set; c this will be called from "ewfind" - if (mode .eq. 1) then + if (mode == 1) then lim1line = 1 lim2line = nlines return @@ -22,9 +22,9 @@ c this will be called from "ewfind" c*****for deriving abundances from single of lines of one species, delimit c the lines of that species as the range; called from "abfind" - elseif (mode .eq. 2) then - if (lim2line .eq. nlines) then - if (nlines .eq. 1) then + elseif (mode == 2) then + if (lim2line == nlines) then + if (nlines == 1) then lim1line = 1 lim2line = 1 mode = -1 @@ -33,18 +33,18 @@ c the lines of that species as the range; called from "abfind" endif return endif - if (lim1line .eq. 0) then + if (lim1line == 0) then lim1line = 1 else lim1line = lim2line + 1 endif - if (lim1line .eq. nlines) then + if (lim1line == nlines) then lim2line = lim1line return else oldatom = atom1(lim1line) do j=lim1line+1,nlines - if (atom1(j) .ne. oldatom) then + if (atom1(j) /= oldatom) then lim2line = j - 1 return endif @@ -58,25 +58,25 @@ c*****for spectrum synthesis, find the range of lines to include at each c wavelength step; called from "synspec"; if requested synthesis c begins more than 10A blueward or goes on more than 105A redward of the c linelist limits, the synthesis aborts with a message; - elseif (mode .eq. 3) then + elseif (mode == 3) then wavelo = wave - delta wavehi = wave + delta c requested synthesis too far from linelist limits - if (wavehi .lt. wave1(1)-10.0) then + if (wavehi < wave1(1)-10.0) then write (*,1004) stop - elseif (wavelo .gt. wave1(nlines)+10.0 .and. - . nlines+nstrong .lt. 2500) then + elseif (wavelo > wave1(nlines)+10.0 .and. + . nlines+nstrong < 2500) then write (*,1005) stop endif c blank synthesis at start or end of requested wavelength range - if (wavehi .lt. wave1(1)) then + if (wavehi < wave1(1)) then lim1line = 1 lim2line = 1 lineflag = -1 return - elseif (wavelo .gt. wave1(nlines)) then + elseif (wavelo > wave1(nlines)) then lim1line = nlines lim2line = nlines lineflag = -1 @@ -85,19 +85,19 @@ c blank synthesis at start or end of requested wavelength range c requested synthesis region is at least partially within the line list c first set the lower synthesis limit; at the beginning of a synthesis c of with a new chunk of lines, lim1line = 0 - if (lim1line.eq.0 .or. wavelo .lt. wave1(1)) lim1line = 1 + if (lim1line==0 .or. wavelo < wave1(1)) lim1line = 1 do j=lim1line,nlines - if (wavelo .lt. wave1(j)) then + if (wavelo < wave1(j)) then lim1line = j exit endif enddo c now set the upper synthesis limit do j=lim1line,nlines - if (wavehi .lt. wave1(j)) then + if (wavehi < wave1(j)) then lim2line = j - 1 - if (lim1line.eq.lim2line .and. - . wavelo.gt. wave1(lim1line)) then + if (lim1line==lim2line .and. + . wavelo> wave1(lim1line)) then lineflag = -1 endif return @@ -105,7 +105,7 @@ c now set the upper synthesis limit enddo c here the end of the line list has been reached; decide whether c to read in more lines - if (nlines+nstrong .eq. 2500) then + if (nlines+nstrong == 2500) then lim2line = -1 else lim2line = nlines @@ -115,20 +115,20 @@ c to read in more lines c*****for blended line force fits to the EWs, the range is a set of c lines in a particular blend - elseif (mode .eq. 4) then - if (lim1line .eq. 0) lim1line = 1 - if (group(lim1line) .ne. 0.) then + elseif (mode == 4) then + if (lim1line == 0) lim1line = 1 + if (group(lim1line) /= 0.) then write (array,1001) call prinfo (10) write (array,1002) call prinfo (11) stop endif - if (lim1line .eq. nlines) then + if (lim1line == nlines) then lim2line = lim1line else do j=lim1line+1,nlines - if (group(j) .ne. 1) then + if (group(j) /= 1) then lim2line = j - 1 return endif @@ -13,8 +13,8 @@ c****************************************************************************** c open the plot device: screen terminal - if (plotroutine(1:4) .eq. 'term') then - if (sm_device(smterm) .lt. 0) then + if (plotroutine(1:4) == 'term') then + if (sm_device(smterm) < 0) then write (array,1001) smterm istat = ivwrite(lscreen+1,1,array,79) write (nf1out,1007) array(1:79) @@ -24,16 +24,16 @@ c open the plot device: screen terminal c open the plot device: hardcopy sent to printer - if (plotroutine(1:4) .eq. 'hard') then - if (plotroutine(6:9) .eq. 'land') then - if (sm_device('postland') .lt. 0) then + if (plotroutine(1:4) == 'hard') then + if (plotroutine(6:9) == 'land') then + if (sm_device('postland') < 0) then write (array,1002) istat = ivwrite(lscreen+1,1,array,34) write (nf1out,1007) array(1:34) stop endif - elseif (plotroutine(6:9) .eq. 'port') then - if (sm_device('postport') .lt. 0) then + elseif (plotroutine(6:9) == 'port') then + if (sm_device('postport') < 0) then write (array,1009) istat = ivwrite(lscreen+1,1,array,34) write (nf1out,1007) array(1:34) @@ -45,8 +45,8 @@ c open the plot device: hardcopy sent to printer c open the plot device: postscript file - if (plotroutine(1:4) .eq. 'file') then - if (f5out .eq. 'optional_output_file') then + if (plotroutine(1:4) == 'file') then + if (f5out == 'optional_output_file') then array = 'Give the file name for the POSTSRIPT plot image: ' nchars = 49 call getasci (nchars,maxline) @@ -55,21 +55,21 @@ c open the plot device: postscript file nchars = 80 call getcount (nchars,f5out) endif - if (plotroutine(6:9) .eq. 'land') then - if (nchars .lt. 10) then + if (plotroutine(6:9) == 'land') then + if (nchars < 10) then write (errmess,1003) nchars else write (errmess,1004) nchars endif - elseif (plotroutine(6:9) .eq. 'port') then - if (nchars .lt. 10) then + elseif (plotroutine(6:9) == 'port') then + if (nchars < 10) then write (errmess,1005) nchars else write (errmess,1006) nchars endif endif write (array,errmess) f5out(1:nchars) - if (sm_device(array(1:nchars+13)) .lt. 0) then + if (sm_device(array(1:nchars+13)) < 0) then write (nf1out,1007) array(1:nchars+9) istat = ivwrite(lscreen+1,1,array,nchars+9) stop @@ -83,26 +83,26 @@ c issue standard beginning commands c call the routine that makes the desired plot - if (plotroutine(11:14) .eq. 'cog ') then + if (plotroutine(11:14) == 'cog ') then call cogplot - elseif (plotroutine(11:14) .eq. 'abun') then + elseif (plotroutine(11:14) == 'abun') then call abunplot - elseif (plotroutine(11:14) .eq. 'spec') then + elseif (plotroutine(11:14) == 'spec') then call specplot - elseif (plotroutine(11:14) .eq. 'bin ') then + elseif (plotroutine(11:14) == 'bin ') then call binplot - elseif (plotroutine(11:14) .eq. 'flux') then + elseif (plotroutine(11:14) == 'flux') then call fluxplot endif c issue standard ending commands; exit normally - if (plotroutine(1:4) .eq. 'file') then + if (plotroutine(1:4) == 'file') then f5out = 'optional_output_file' endif call sm_gflush - if (plotroutine(1:4).eq.'hard' .or. - . plotroutine(1:4).eq.'file') call sm_hardcopy + if (plotroutine(1:4)=='hard' .or. + . plotroutine(1:4)=='file') call sm_hardcopy call sm_alpha return @@ -17,12 +17,12 @@ c****************************************************************************** c*****the species is an atom: c*****search for it in the list of elements done in M.E. - if (atom1(lim1line) .lt. 100.) then - if (neq .eq. 0) then + if (atom1(lim1line) < 100.) then + if (neq == 0) then return else do n=1,neq - if (iabatom .eq. iorder(n)) molflag = 1 + if (iabatom == iorder(n)) molflag = 1 return enddo endif @@ -31,7 +31,7 @@ c*****search for it in the list of elements done in M.E. c*****the species is a molecule: c*****halt if M.E. wasn't done or didn't include this species else - if (neq .eq. 0) then + if (neq == 0) then lscreen = lscreen + 2 write (array,1001) iabatom call prinfo (lscreen) @@ -41,9 +41,9 @@ c*****halt if M.E. wasn't done or didn't include this species iaa = ia ibb = ib do n=1,neq - if (ia.eq.iorder(n) .or. ib.eq.iorder(n)) molflag = ia + if (ia==iorder(n) .or. ib==iorder(n)) molflag = ia enddo - if (molflag .eq. 0) then + if (molflag == 0) then lscreen = lscreen + 2 write (array,1002) iabatom call prinfo (lscreen) @@ -52,8 +52,8 @@ c*****halt if M.E. wasn't done or didn't include this species molflag = 1 c*****if molecule is a hydride, the non-H element will be varied - if (ia.eq.1 .or. ib.eq.1) then - if (ia .eq. 1) then + if (ia==1 .or. ib==1) then + if (ia == 1) then iabatom = ib else iabatom = ia @@ -67,7 +67,7 @@ c*****for other molecules, the user specifies which element will be varied call getnum (nchars,ikount+1,xnum,shortnum) iabatom = int(xnum+0.0001) - if (iabatom.ne.ia .and. iabatom.ne.ib) then + if (iabatom/=ia .and. iabatom/=ib) then write (array,1003) stop endif @@ -53,40 +53,40 @@ c*****invoke the overall starting routine c*****use one of the standard driver routines - if (control .eq. 'synplot') then + if (control == 'synplot') then call plotit - elseif (control .eq. 'synth ') then + elseif (control == 'synth ') then call synth - elseif (control .eq. 'cogsyn ') then + elseif (control == 'cogsyn ') then call cogsyn - elseif (control .eq. 'blends ') then + elseif (control == 'blends ') then call blends - elseif (control .eq. 'abfind ') then + elseif (control == 'abfind ') then call abfind - elseif (control .eq. 'ewfind ') then + elseif (control == 'ewfind ') then call ewfind - elseif (control .eq. 'cog ') then + elseif (control == 'cog ') then call cog - elseif (control .eq. 'calmod ') then + elseif (control == 'calmod ') then call calmod - elseif (control .eq. 'doflux ') then + elseif (control == 'doflux ') then call doflux - elseif (control .eq. 'weedout') then + elseif (control == 'weedout') then call weedout - elseif (control .eq. 'gridsyn') then + elseif (control == 'gridsyn') then call gridsyn - elseif (control .eq. 'gridplo') then + elseif (control == 'gridplo') then call gridplo - elseif (control .eq. 'binary ') then + elseif (control == 'binary ') then call binary - elseif (control .eq. 'abpop ') then + elseif (control == 'abpop ') then call abpop - elseif (control .eq. 'synpop ') then + elseif (control == 'synpop ') then call synpop c*****or, put in your own drivers in the form below.... - elseif (control .eq. 'mine ') then + elseif (control == 'mine ') then call mydriver @@ -96,7 +96,7 @@ c*****or else you are out of luck! istat = ivwrite (4,3,array,49) istat = ivmove (3,1) read (*,*) yesno - if (yesno .eq. 'y') then + if (yesno == 'y') then go to 1 else call finish (0) diff --git a/Moogsilent.f b/Moogsilent.f index 61e9bef..d8b5b64 100755 --- a/Moogsilent.f +++ b/Moogsilent.f @@ -52,45 +52,45 @@ c*****invoke the overall starting routine c*****use one of the standard driver routines ("isotop" is obsolete): - if (control .eq. 'synplot') then + if (control == 'synplot') then call plotit - elseif (control .eq. 'isoplot') then + elseif (control == 'isoplot') then call plotit - elseif (control .eq. 'synth ') then + elseif (control == 'synth ') then call synth - elseif (control .eq. 'cogsyn ') then + elseif (control == 'cogsyn ') then call cogsyn - elseif (control .eq. 'blends ') then + elseif (control == 'blends ') then call blends - elseif (control .eq. 'abfind ') then + elseif (control == 'abfind ') then call abfind - elseif (control .eq. 'ewfind ') then + elseif (control == 'ewfind ') then call ewfind - elseif (control .eq. 'cog ') then + elseif (control == 'cog ') then call cog - elseif (control .eq. 'calmod ') then + elseif (control == 'calmod ') then call calmod - elseif (control .eq. 'isotop ') then + elseif (control == 'isotop ') then control = 'synth ' call synth - elseif (control .eq. 'doflux ') then + elseif (control == 'doflux ') then call doflux - elseif (control .eq. 'weedout') then + elseif (control == 'weedout') then call weedout - elseif (control .eq. 'gridsyn') then + elseif (control == 'gridsyn') then call gridsyn - elseif (control .eq. 'gridplo') then + elseif (control == 'gridplo') then call gridplo - elseif (control .eq. 'binary ') then + elseif (control == 'binary ') then call binary - elseif (control .eq. 'abpop ') then + elseif (control == 'abpop ') then call abpop - elseif (control .eq. 'synpop ') then + elseif (control == 'synpop ') then call synpop c*****or, put in your own drivers in the form below.... - elseif (control .eq. 'mine ') then + elseif (control == 'mine ') then call mydriver @@ -42,7 +42,7 @@ c equivalence (dummy,string(1)) c count = ccount - if (y .gt. maxline .or. x .gt. 79) then + if (y > maxline .or. x > 79) then ivwrite = -1 return endif @@ -51,8 +51,8 @@ c dummy = arr count = min0(80-x,count) c - if (x .lt. 10) then - if (y .lt. 10) then + if (x < 10) then + if (y < 10) then write (*,1007) esc,y,x,(string(i),i=1,count) 1007 format(1x,a1,'[',i1,';',i1,'H',80a1) else @@ -60,7 +60,7 @@ c 1006 format(1x,a1,'[',i2,';',i1,'H',80a1) endif else - if (y .lt. 10) then + if (y < 10) then write (*,1005) esc,y,x,(string(i),i=1,count) 1005 format(1x,a1,'[',i1,';',i2,'H',80a1) else @@ -88,15 +88,15 @@ c integer y,x character esc*1 c - if (y .gt. maxline .or. x .gt. 79) then + if (y > maxline .or. x > 79) then ivmove = -1 return endif c esc = char(27) c - if (x .lt. 10) then - if (y .lt. 10) then + if (x < 10) then + if (y < 10) then write (*,1007) esc,y,x 1007 format(1x,a1,'[',i1,';',i1,'H') else @@ -104,7 +104,7 @@ c 1006 format(1x,a1,'[',i2,';',i1,'H') endif else - if (y .lt. 10) then + if (y < 10) then write (*,1005) esc,y,x 1005 format(1x,a1,'[',i1,';',i2,'H') else @@ -148,7 +148,7 @@ c character*1 esc c esc = char(27) - if(on.eq.1) then + if(on==1) then write(*,'(1x,a1,a)') esc,'[1m' else write(*,'(1x,a1,a)') esc,'[0m' @@ -171,7 +171,7 @@ c character*1 esc c esc = char(27) - if(on.eq.1) then + if(on==1) then write(*,'(1x,a1,a)') esc,'[4m' else write(*,'(1x,a1,a)') esc,'[0m' @@ -23,14 +23,14 @@ c****************************************************************************** c*****load in data for damping factors to be computed from the data c of Barklem, if desired - if (numpass.eq.1 .and. dampingopt.eq.1) call gammabark + if (numpass==1 .and. dampingopt==1) call gammabark c*****Locate the atmosphere level where taulam is near 0.5 - if (numpass.eq.1 .or. numpass.eq.3) then + if (numpass==1 .or. numpass==3) then call opacit (2,wave1(1)) do i=1,ntau - if (taulam(i) .ge. 0.5) go to 180 + if (taulam(i) >= 0.5) go to 180 enddo 180 jtau5 = i endif @@ -41,13 +41,13 @@ c if iterating "abfind" on a species affected by mol. eq., c do just that species (numpass=2); if doing a fake line, or c maybe one line for a special purpose, do just the first line c in the "list" (numpass=3) - if (numpass .eq. 1) then + if (numpass == 1) then j1 = 1 j2 = nlines + nstrong - elseif (numpass .eq. 2) then + elseif (numpass == 2) then j1 = lim1line j2 = lim2line - elseif (numpass .eq. 3) then + elseif (numpass == 3) then j1 = 1 j2 = 1 endif @@ -61,7 +61,7 @@ c*****now make the calculations: set up some parameters c*****compute the Doppler factors - if (numpass.eq.1 .or. numpass.eq.3) then + if (numpass==1 .or. numpass==3) then do i=1,ntau dopp(j,i) = dsqrt(1.6631d8*t(i)/amass(j)+vturb(i)**2) enddo @@ -76,7 +76,7 @@ c*****either: compute lower state number densities for atomic lines; c q21 is the ion/neutral ratio, etc., and q is the ratio of the total c to the species of interest; do the Saha equation first, then c the Boltzmann equation - if (iatom .lt. 100) then + if (iatom < 100) then do i=1,ntau q21 = 4.825d15*u(iatom,2,i)/(u(iatom,1,i)*ne(i))* . t(i)**1.5*dexp(-chi(j,1)/tkev(i)) @@ -84,15 +84,15 @@ c the Boltzmann equation . t(i)**1.5*dexp(-chi(j,2)/tkev(i)) q43 = 4.825d15*u(iatom,4,i)/(u(iatom,3,i)*ne(i))* . t(i)**1.5*dexp(-chi(j,3)/tkev(i)) - if (neq .ne. 0) then + if (neq /= 0) then do n=1,neq - if (iatom .eq. iorder(n)) then + if (iatom == iorder(n)) then xxnum = xamol(n,i) - if (ich .eq. 1) then + if (ich == 1) then q = 1.0 - elseif (ich .eq. 2) then + elseif (ich == 2) then q = 1.0/q21 - elseif (ich .eq. 3) then + elseif (ich == 3) then q = 1.0/(q21*q32) + 1.0/q32 + 1.0 + q43 endif xnum(i) = xxnum/q*dexp(-e(j,1)/tkev(i))/ @@ -100,14 +100,14 @@ c the Boltzmann equation endif enddo endif - if (ich .eq. 1) then + if (ich == 1) then q = 1.0 + q21 + q32*q21 + q43*q32*q21 - elseif (ich .eq. 2) then + elseif (ich == 2) then q = 1.0/q21 + 1.0 + q32 + q43*q32 - elseif (ich .eq. 3) then + elseif (ich == 3) then q = 1.0/(q21*q32) + 1.0/q32 + 1.0 + q43 endif - if (control .eq. 'abandy') then + if (control == 'abandy') then xxab = xabund(iatom)*10**deltaabund else xxab = xabund(iatom) @@ -118,11 +118,11 @@ c the Boltzmann equation c*****or: compute lower state number densities for molecular lines - elseif (iatom .lt. 10000) then + elseif (iatom < 10000) then call sunder(atom1(j),iaa,ibb) do n=1,neq - if(iorder(n) .eq. iaa) ia = n - if(iorder(n) .eq. ibb) ib = n + if(iorder(n) == iaa) ia = n + if(iorder(n) == ibb) ib = n enddo do i=1,ntau psipri = @@ -133,11 +133,11 @@ c*****or: compute lower state number densities for molecular lines . xamol(ia,i)*xamol(ib,i) enddo else - if (iatom .eq. 10108) then + if (iatom == 10108) then do i=1,ntau xnum(i) = xnh2o(i)/uh2o(i)*dexp(-e(j,1)/tkev(i)) enddo - elseif (iatom .eq. 60808) then + elseif (iatom == 60808) then do i=1,ntau xnum(i) = xnco2(i)/uco2(i)*dexp(-e(j,1)/tkev(i)) enddo @@ -149,9 +149,9 @@ c*****or: compute lower state number densities for molecular lines c*****finally, compute line opacities at line centers - if (atom1(j)-float(iatom) .ge. 0.0) then + if (atom1(j)-float(iatom) >= 0.0) then do n=1,numiso - if (atom1(j) .eq. isotope(n)) then + if (atom1(j) == isotope(n)) then factoriso = isoabund(n,isorun) endif enddo @@ -167,9 +167,9 @@ c*****finally, compute line opacities at line centers c*****output regular line information, and strong line information c if appropriate; exit routine normally - if (numpass.eq.1 .or. numpass.eq.3) then - if (linprintopt .ge. 0) then - if (dostrong .gt. 0) call lineinfo (2) + if (numpass==1 .or. numpass==3) then + if (linprintopt >= 0) then + if (dostrong > 0) call lineinfo (2) call lineinfo (1) endif endif @@ -11,14 +11,14 @@ c****************************************************************************** c*****if a carriage return has been hit, return with -9999. - if (nchar .le. 0) then + if (nchar <= 0) then xnum = -9999. return endif c*****set the conversion format - if (nchar .lt. 10) then + if (nchar < 10) then write(form,1001) nchar else write(form,1002) nchar @@ -13,66 +13,66 @@ c****************************************************************************** do j=1,36 k = 80*(j-1) - if (head(k+1:k+8) .eq. 'SIMPLE ') then - if (head(k+30:k+30) .ne. 'T') then + if (head(k+1:k+8) == 'SIMPLE ') then + if (head(k+30:k+30) /= 'T') then write(array,1029) head(k+1:k+58) istat = ivwrite (line+2,3,array,79) go to 1007 endif - elseif (head(k+1:k+8) .eq. 'BITPIX ') then + elseif (head(k+1:k+8) == 'BITPIX ') then read (head(k+1:k+80),1025) ibits - if (ibits .eq. 16) then + if (ibits == 16) then nblock = 1440 - elseif (ibits .eq. 32) then + elseif (ibits == 32) then nblock = 720 - elseif (ibits .eq. -32) then + elseif (ibits == -32) then nblock = 720 else write(array,1026) ibits istat = ivwrite (line+2,3,array,32) go to 1007 endif - elseif (head(k+1:k+8) .eq. 'NAXIS ') then + elseif (head(k+1:k+8) == 'NAXIS ') then read (head(k+1:k+80),1025) naxis - if (naxis .ne. 1) then + if (naxis /= 1) then write(array,1028) head(k+1:k+58) go to 1007 endif - elseif (head(k+1:k+8) .eq. 'NAXIS1 ') then + elseif (head(k+1:k+8) == 'NAXIS1 ') then read (head(k+1:k+80),1025) lount - elseif (head(k+1:k+8) .eq. 'OBJECT ') then + elseif (head(k+1:k+8) == 'OBJECT ') then write (obsitle,1027) head(k+12:k+80) - elseif (head(k+1:k+8) .eq. 'BZERO ') then + elseif (head(k+1:k+8) == 'BZERO ') then read (head(k+1:k+80),1024) bzero - elseif (head(k+1:k+8) .eq. 'BSCALE ') then + elseif (head(k+1:k+8) == 'BSCALE ') then read (head(k+1:k+80),1024) bscale - elseif ((head(k+1:k+8) .eq. 'W0 ') .or. - . (head(k+1:k+8) .eq. 'CRVAL1 ')) then + elseif ((head(k+1:k+8) == 'W0 ') .or. + . (head(k+1:k+8) == 'CRVAL1 ')) then read (head(k+1:k+80),1024) disp(1) - elseif ((head(k+1:k+8) .eq. 'WPC ') .or. - . (head(k+1:k+8) .eq. 'CDELT1 ')) then + elseif ((head(k+1:k+8) == 'WPC ') .or. + . (head(k+1:k+8) == 'CDELT1 ')) then read (head(k+1:k+80),1024) dval - if (dval .ne. 1.) disp(2) = dval - elseif (head(k+1:k+8) .eq. 'CD1_1 ') then + if (dval /= 1.) disp(2) = dval + elseif (head(k+1:k+8) == 'CD1_1 ') then read (head(k+1:k+80),1024) disp(2) - elseif (head(k+1:k+8) .eq. 'FILENAME') then + elseif (head(k+1:k+8) == 'FILENAME') then write (obsitle(39:80),1023) head(k+12:k+53) - elseif (head(k+1:k+8) .eq. 'HISTORY ') then - if (head(k+24:k+28) .eq. 'DISP=') then + elseif (head(k+1:k+8) == 'HISTORY ') then + if (head(k+24:k+28) == 'DISP=') then read (head(k+1:k+80),1022) (disp(i),i=1,4) - elseif (head(k+20:k+26) .eq. 'D1,2,3:') then + elseif (head(k+20:k+26) == 'D1,2,3:') then read (head(k+1:k+80),1042) (disp(i),i=1,3) - elseif (head(k+20:k+26) .eq. 'D4,5,6:') then + elseif (head(k+20:k+26) == 'D4,5,6:') then read (head(k+1:k+80),1042) (disp(i),i=4,6) - elseif (head(k+20:k+26) .eq. 'D7,8,9:') then + elseif (head(k+20:k+26) == 'D7,8,9:') then read (head(k+1:k+80),1042) (disp(i),i=7,9) - if (disp(7).ne.0.0 .and. disp(8).eq.0.0 .and. - . disp(9).eq.0.0) then + if (disp(7)/=0.0 .and. disp(8)==0.0 .and. + . disp(9)==0.0) then disp(8) = 1.0 disp(9) = lount endif endif - elseif (head(k+1:k+8) .eq. 'END ') then + elseif (head(k+1:k+8) == 'END ') then iend = 1 return endif @@ -19,9 +19,9 @@ c the line profile c*****get started; calculate an initial step size; wavestep only is c*****used in synpop - if (imode .eq. 0) gf1(ncurve) = gf(lim1) + if (imode == 0) gf1(ncurve) = gf(lim1) dellam(1) = 0. - if (wavestep .eq. 0.) then + if (wavestep == 0.) then st1 = wave1(lim1)*dopp(lim1,jtau5)/2.997929e10/5. st1 = dfloat(ifix(10000.*sngl(st1)))/10000. else @@ -32,16 +32,16 @@ c*****used in synpop c*****calculate continuous opacity and intensity/flux at line wavelength wave = wave1(lim1) - if (abs(wave-waveold) .gt. 30.) then + if (abs(wave-waveold) > 30.) then waveold = wave call opacit(2,wave) - if (imode.ne.2 .and. modprintopt.ge.2) + if (imode/=2 .and. modprintopt>=2) . write(nf1out,1002) wave,(kaplam(i),i=1,ntau) call cdcalc(1) first = 0.4343*cd(1) flux = rinteg(xref,cd,dummy1,ntau,first) - if (imode .ne. 2) then - if (iunits .eq. 1) then + if (imode /= 2) then + if (iunits == 1) then write (nf1out,1003) 1.d-4*wave,flux else write (nf1out,1004) wave,flux @@ -51,14 +51,14 @@ c*****calculate continuous opacity and intensity/flux at line wavelength c*****check the wavelength step size; expand/contract as necessary - if (wavestep .eq. 0.) then + if (wavestep == 0.) then wave = wave1(lim1) call taukap call cdcalc(2) first = 0.4343*cd(1) d(1) = rinteg(xref,cd,dummy1,ntau,first) do k=1,30 - if (k .eq. 30) then + if (k == 30) then write (*,1010) wave stop endif @@ -68,17 +68,17 @@ c*****check the wavelength step size; expand/contract as necessary first = 0.4343*cd(1) d(2) = rinteg(xref,cd,dummy1,ntau,first) d2d1 = d(2)/d(1) - if (d2d1 .le. 0.2) then + if (d2d1 <= 0.2) then st1 = st1/1.5 - elseif (d2d1 .le. 0.6) then + elseif (d2d1 <= 0.6) then st1 = st1/1.2 - elseif (d2d1.gt.0.60 .and. d2d1.lt.0.80) then - if (imode .ne. 2) write (nf1out,1001) lim1, + elseif (d2d1>0.60 .and. d2d1<0.80) then + if (imode /= 2) write (nf1out,1001) lim1, . 1000.*width(lim1), storig, st1, k exit - elseif (d2d1 .ge. 0.80) then + elseif (d2d1 >= 0.80) then st1 = st1*1.6 - elseif (d2d1 .ge. 0.90) then + elseif (d2d1 >= 0.90) then st1 = st1*2.1 endif st1 = dble(idnint(10000.*st1))/10000. @@ -96,7 +96,7 @@ c until the depth is very small in the line wing call cdcalc(2) first = 0.4343*cd(1) d(n) = rinteg(xref,cd,dummy1,ntau,first) - if (linprintopt.ge.3 .and. n.eq.1 .and. imode.eq.2) then + if (linprintopt>=3 .and. n==1 .and. imode==2) then do i=1,ntau dummy1(i) = xref(i)*cd(i) enddo @@ -104,21 +104,21 @@ c until the depth is very small in the line wing cdmean = rinteg(xref,dummy1,dummy2,ntau,first)/ . rinteg(xref,cd,dummy2,ntau,first) do i=1,ntau - if (cdmean .lt. cd(i)) exit + if (cdmean < cd(i)) exit enddo write (nf1out,1005) lim1, cdmean, i, xref(i) do i=1,ntau - if (taunu(i)+taulam(i) .ge. 1.) exit + if (taunu(i)+taulam(i) >= 1.) exit enddo write (nf1out,1006) lim1, i, dlog10(tauref(i)), . dlog10(taulam(i)), dlog10(taunu(i)) endif - if (d(n)/d(1) .lt. 0.0050) then + if (d(n)/d(1) < 0.0050) then ndepths = n exit endif - if (n .eq. maxsteps) then - if (d(n).gt.0.001 .and. imode.ne.2) then + if (n == maxsteps) then + if (d(n)>0.001 .and. imode/=2) then write (nf1out,1007) ndepths = maxsteps endif @@ -127,7 +127,7 @@ c until the depth is very small in the line wing c*****finish the calculation - if (imode .ne. 2) write (nf1out,1008) (d(j),j=1,ndepths) + if (imode /= 2) write (nf1out,1008) (d(j),j=1,ndepths) do n=2,ndepths d(ndepths+n-1) = d(n) enddo @@ -142,7 +142,7 @@ c*****finish the calculation enddo first = 2*dellam(ndep)*d(ndep) w(ncurve) = rinteg(dellam,d,dinteg,ndep,first) - if (imode .ne. 2) then + if (imode /= 2) then ew = 1000.0*w(ncurve) gflog = dlog10(gf1(ncurve)) rwlog = dlog10(w(ncurve)/wave1(lim1)) diff --git a/OpacHelium.f b/OpacHelium.f index b02fc91..294d8fa 100755 --- a/OpacHelium.f +++ b/OpacHelium.f @@ -19,7 +19,7 @@ c****************************************************************************** save data freq1 /0./ - if (freq .ne. freq1) then + if (freq /= freq1) then freq1 = freq a1 = 3.397d-46 + (-5.216d-31+7.039d-15/freq)/freq b1 = -4.116d-42 + ( 1.067d-26+8.135d-11/freq)/freq diff --git a/OpacHydrogen.f b/OpacHydrogen.f index a193fb9..6ba3b1d 100755 --- a/OpacHydrogen.f +++ b/OpacHydrogen.f @@ -22,7 +22,7 @@ c****************************************************************************** data modcount/0/ c set up some data upon first entrance with a new model atmosphere - if (modelnum .ne. modcount) then + if (modelnum /= modcount) then modcount = modelnum do i=1,ntau do n=1,8 @@ -46,7 +46,7 @@ c set up some data upon first entrance with a new model atmosphere c = 2.815d29/freq3 do i=1,ntau ex = boltex(i) - if (freq .lt. 4.05933d13) ex = exlim(i)/evhkt(i) + if (freq < 4.05933d13) ex = exlim(i)/evhkt(i) h = (cont(7)*bolt(i,7) + cont8*bolt(i,8) + . (ex - exlim(i))*c + . coulff(1,tlog(i),freq)*freet(i)*cfree)*(1.-evhkt(i)) @@ -134,7 +134,7 @@ C Bell and Berrington (1987, J.Phys.B, 20, 801-806) data modcount,istart/ 0,0/ c fill some arrays once and for all - if (istart .eq. 0) then + if (istart == 0) then istart = 1 c 91.134 number taken from Bell & Berrington do iwave=1,22 @@ -147,7 +147,7 @@ c 91.134 number taken from Bell & Berrington c initialize some quantities for each new model atmosphere c .754209 Hotop & Lineberger (1985, J. Phys. Chem. Ref. Data, 14,731-752) - if (modelnum .ne. modcount) then + if (modelnum /= modcount) then modcount = modelnum do i=1,ntau xhmin(i) = dexp(.754209/tkev(i))/(2.*2.4148d15*t(i)**1.5)* @@ -165,7 +165,7 @@ c main opacity computation yielding "aHminus" enddo hminbf = 0. nnnn = 85 - if (freq .gt. 1.82365d14) + if (freq > 1.82365d14) . maxwave = map1(wbf,bf,nnnn,wave,hminbf,1) do i=1,ntau nnnn = 11 @@ -193,7 +193,7 @@ c****************************************************************************** iold = 2 do inew=1,nnew -1 if (xnew(inew).lt.xold(iold) .or. iold.eq.nold) then +1 if (xnew(inew)<xold(iold) .or. iold==nold) then ynew(inew) = yold(iold-1) + (yold(iold)-yold(iold-1))/ . (xold(iold)-xold(iold-1))*(xnew(inew)-xold(iold-1)) return @@ -222,20 +222,20 @@ c****************************************************************************** l = 2 ll = 0 do 50 k=1,nnew -10 if(xnew(k) .lt. xold(l)) go to 20 +10 if(xnew(k) < xold(l)) go to 20 l = l + 1 - if (l .gt. nold) go to 30 + if (l > nold) go to 30 go to 10 -20 if (l .eq. ll) go to 50 - if (l .eq. 2) go to 30 - if (l .eq. 3) go to 30 +20 if (l == ll) go to 50 + if (l == 2) go to 30 + if (l == 3) go to 30 l1 = l - 1 - if (l.gt.ll+1 .or. l.eq.3) go to 21 - if (l.gt.ll+1 .or. l.eq.4) go to 21 + if (l>ll+1 .or. l==3) go to 21 + if (l>ll+1 .or. l==4) go to 21 cbac = cfor bbac = bfor abac = afor - if (l .eq. nold) go to 22 + if (l == nold) go to 22 go to 25 21 l2 = l - 2 d = (fold(l1)-fold(l2))/(xold(l1)-xold(l2)) @@ -244,7 +244,7 @@ c****************************************************************************** . (xold(l1)-xold(l2)) bbac = d - (xold(l1)+xold(l2))*cbac abac = fold(l2) - xold(l2)*d + xold(l1)*xold(l2)*cbac - if (l .lt. nold) go to 25 + if (l < nold) go to 25 22 c = cbac b = bbac a = abac @@ -257,13 +257,13 @@ c****************************************************************************** bfor = d - (xold(l)+xold(l1))*cfor afor = fold(l1) - xold(l1)*d + xold(l)*xold(l1)*cfor wt = 0. - if (dabs(cfor) .ne. 0.) wt = dabs(cfor)/(dabs(cfor)+dabs(cbac)) + if (dabs(cfor) /= 0.) wt = dabs(cfor)/(dabs(cfor)+dabs(cbac)) a = afor + wt*(abac-afor) b = bfor + wt*(bbac-bfor) c = cfor + wt*(cbac-cfor) ll = l go to 50 -30 if (l .eq. ll) go to 50 +30 if (l == ll) go to 50 l = min0(nold,l) c = 0. b = (fold(l)-fold(l-1))/(xold(l)-xold(l-1)) diff --git a/Opaccouls.f b/Opaccouls.f index 3c4f59d..2410af8 100755 --- a/Opaccouls.f +++ b/Opaccouls.f @@ -19,12 +19,12 @@ c******************************************************************************` data c/ -2.268d30, 4.077d28, 1.035d28, 4.593d27, 2.371d27, . 1.229d27/ - if (freq .lt. z*z*3.28805d15/dfloat(n*n)) then + if (freq < z*z*3.28805d15/dfloat(n*n)) then coulx = 0. else coulx = 2.815d29/freq**3/dfloat(n**5)*z**4 - if (n .gt. 6) return - if (n .eq. 1) then + if (n > 6) return + if (n == 1) then coulx = coulx*coulbf1s(freq,z) else coulx = coulx*(a(n)+(b(n)+c(n)*(z*z/freq))*(z*z/freq)) @@ -64,7 +64,7 @@ c******************************************************************************` . 0.2701,0.2648,0.2595,0.2544,0.2493,0.2443,0.2394,0.2345,0.2298, . 0.2251,0.2205,0.2160,0.2115,0.2072,0.2029,0.1987/ - if (freq/z**2 .lt. 3.28805d15) then + if (freq/z**2 < 3.28805d15) then coulbf1s = 0. else elog = log10(freq/z**2/3.28805d15) @@ -62,7 +62,7 @@ c*****sum up all the opacities c*****write out the opacities - if (modprintopt .gt. 1) then + if (modprintopt > 1) then write (nf1out,1001) waveop do i=1,ntau write (nf1out,1002) i, nint(t(i)), dlog10(kaplam(i)), @@ -82,10 +82,10 @@ c*****write out the opacities c*****add in an arbitrary amount of "extra" continuous opacity; c this is a pure fudge factor, so experiment with care - if (fudge .gt. 0.0) then + if (fudge > 0.0) then do i=1,ntau kaplam(i) = kaplam(i)*((fudge*10000)/t(i)) - if (i .eq. 1) then + if (i == 1) then write(nf1out,1005) write(nf1out,1006) fudge endif @@ -94,7 +94,7 @@ c this is a pure fudge factor, so experiment with care c*****compute an optical depth array at this wavelength, and exit - if (modeop .ne. 1) then + if (modeop /= 1) then do i=1,ntau dummy1(i) = tauref(i)*kaplam(i)/(0.4343*kapref(i)) enddo @@ -104,7 +104,7 @@ c*****compute an optical depth array at this wavelength, and exit do i=2,ntau taulam(i) = taulam(i-1) + taulam(i) enddo - if (modprintopt .gt. 1) write(nf1out,1004) (taulam(i),i=1,ntau) + if (modprintopt > 1) write(nf1out,1004) (taulam(i),i=1,ntau) return endif @@ -113,7 +113,7 @@ c*****here is the assignment of opacities kaplam to kapref; exit normally do i=1,ntau kapref(i) = kaplam(i) enddo - if(modprintopt .lt. 2) return + if(modprintopt < 2) return c*****here is an internal tauref check on the externally read in tau scale; diff --git a/Opacmetals.f b/Opacmetals.f index 41a24e1..9901a50 100755 --- a/Opacmetals.f +++ b/Opacmetals.f @@ -21,7 +21,7 @@ c****************************************************************************** data freq1, modcount/0. ,0/ c initialize some quantities for each new model atmosphere - if (modelnum .ne. modcount) then + if (modelnum /= modcount) then modcount = modelnum do i=1,ntau c1240(i) = 5.*dexp(-1.264/tkev(i)) @@ -32,7 +32,7 @@ c initialize some quantities for each new model atmosphere c initialize some quantities for each new model atmosphere or new frequency; c Luo, D. and Pradhan, A.K. 1989, J.Phys. B, 22, 3377-3395. c Burke, P.G. and Taylor, K.T. 1979, J. Phys. B, 12, 2971-2984. - if (modelnum.ne.modcount .or. freq.ne.freq1) then + if (modelnum/=modcount .or. freq/=freq1) then freq1 = freq ryd = 109732.298 waveno = freq/2.99792458d10 @@ -47,11 +47,11 @@ c Burke, P.G. and Taylor, K.T. 1979, J. Phys. B, 12, 2971-2984. c P2 3P 1 c I AM NOT SURE WHETHER THE CALL TO SEATON IN THE NEXT STATEMENT IS c CORRECT, BUT IT ONLY AFFECTS THINGS BELOW 1100A - if (freq .ge. 2.7254d15) x1100 = + if (freq >= 2.7254d15) x1100 = . 10.**(-16.80-(waveno-90777.000)/3./ryd)* . seaton (2.7254d15,1.219d-17,2.d0,3.317d0) c P2 1D 2 - if (freq .ge. 2.4196d15) then + if (freq >= 2.4196d15) then xd0 = 10.**(-16.80-(waveno-80627.760)/3./ryd) eeps = (waveno-93917.)*2./9230. aa = 22.d-18 @@ -64,7 +64,7 @@ c P2 1D 2 x1240 = xd0 + xd1 + xd2 endif c P2 1S 3 - if (freq .ge. 2.0761d15) then + if (freq >= 2.0761d15) then xs0 = 10.**(-16.80-(waveno-69172.400)/3./ryd) eeps = (waveno-97700.)*2./2743. aa = 68.d-18 @@ -75,7 +75,7 @@ c P2 1S 3 endif do i=1,ntau - if (freq .ge. 2.0761d15) then + if (freq >= 2.0761d15) then aC1(i) = (x1100*9. + x1240*c1240(i) + x1444*c1444(i))* . numdens(3,1,i)/u(6,1,i) endif @@ -130,7 +130,7 @@ c 4000 5000 6000 7000 8000 9000 10000 data freq1, modcount/0. ,0/ c initialize some quantities for each new model atmosphere - if (modelnum .ne. modcount) then + if (modelnum /= modcount) then modcount = modelnum do i=1,ntau n = max0(min0(6,nint(t(i)/1000.)-3),1) @@ -140,14 +140,14 @@ c initialize some quantities for each new model atmosphere endif c initialize some quantities for each new model atmosphere or new frequency; - if (modelnum.ne.modcount .or. freq.ne.freq1) then + if (modelnum/=modcount .or. freq/=freq1) then freq1 = freq do n=1,7 - if (freq .gt. freqMg(n)) go to 23 + if (freq > freqMg(n)) go to 23 enddo n = 8 23 dd = (freqlg-flog(n))/(flog(n+1)-flog(n)) - if (n .gt. 2) n = 2*n -2 + if (n > 2) n = 2*n -2 dd1 = 1.0 - dd do it=1,7 xx(it) = peach(it,n+1)*dd + peach(it,n)*dd1 @@ -155,7 +155,7 @@ c initialize some quantities for each new model atmosphere or new frequency; endif do i=1,ntau - if (freq .ge. 2.997925d+14) then + if (freq >= 2.997925d+14) then n = nt(i) aMg1(i) = dexp(xx(n)*(1.d0-dt(i))+xx(n+1)*dt(i))* . numdens(4,1,i)/u(12,1,i) @@ -184,7 +184,7 @@ c****************************************************************************** data freq1, modcount/0. ,0/ c initialize some quantities for each new model atmosphere - if (modelnum .ne. modcount) then + if (modelnum /= modcount) then modcount = modelnum do i=1,ntau c1169(i) = 6.*dexp(-4.43d+0/tkev(i)) @@ -193,14 +193,14 @@ c initialize some quantities for each new model atmosphere c initialize some quantities for each new model atmosphere or new frequency; c there are two edges, one at 824 A and the other at 1169 A - if (modelnum.ne.modcount .or. freq.ne.freq1) then + if (modelnum/=modcount .or. freq/=freq1) then freq1 = freq - if (freq .ge. 3.635492d15) then + if (freq >= 3.635492d15) then x824 = seaton (3.635492d15,1.40d-19,4.d0,6.7d0) else x824 = 1.d-99 endif - if (freq .ge. 2.564306d15) then + if (freq >= 2.564306d15) then x1169 = 5.11d-19*(2.564306d15/freq)**3 else x1169 = 1.d-99 @@ -208,7 +208,7 @@ c there are two edges, one at 824 A and the other at 1169 A endif do i=1,ntau - if (x1169 .ge. 1.d-90) then + if (x1169 >= 1.d-90) then aMg2(i) = (x824*2. + x1169*c1169(i))* . numdens(4,2,i)/u(12,2,i) endif @@ -234,7 +234,7 @@ c****************************************************************************** include 'Kappa.com' do i=1,ntau - if (freq .ge. 1.443d15) then + if (freq >= 1.443d15) then aAl1(i) = 6.5d-17*(1.443d15/freq)**5*6.* . numdens(5,1,i)/u(13,1,i) endif @@ -296,7 +296,7 @@ c 3P, 1D, 1S, 1D, 3D, 3F, 1D, 3P data freq1, modcount/0. ,0/ c initialize some quantities for each new model atmosphere - if (modelnum .ne. modcount) then + if (modelnum /= modcount) then modcount = modelnum do i=1,ntau n = max0(min0(8,nint(t(i)/1000.)-3),1) @@ -306,14 +306,14 @@ c initialize some quantities for each new model atmosphere endif c initialize some quantities for each new model atmosphere or new frequency - if (modelnum.ne.modcount .or. freq.ne.freq1) then + if (modelnum/=modcount .or. freq/=freq1) then freq1 = freq do n=1,9 - if (freq .gt. freqSi(n)) go to 23 + if (freq > freqSi(n)) go to 23 enddo n = 10 23 dd = (freqlg-flog(n))/(flog(n+1)-flog(n)) - if (n .gt. 2) n = 2*n - 2 + if (n > 2) n = 2*n - 2 dd1 = 1.0 - dd do it=1,9 xx(it) = peach(it,n+1)*dd + peach(it,n)*dd1 @@ -321,7 +321,7 @@ c initialize some quantities for each new model atmosphere or new frequency endif do i=1,ntau - if (freq .ge. 2.997925d+14) then + if (freq >= 2.997925d+14) then n = nt(i) aSi1(i) = (dexp(-(xx(n)*(1.-dt(i)) + xx(n+1)*dt(i)))*9.)* . numdens(6,1,i)/u(14,1,i) @@ -373,7 +373,7 @@ c 2P, 2D, 2P, 2D, 2P data freq1, modcount/0., 0/ c set up some data upon first entrance with a new model atmosphere - if (modelnum .ne. modcount) then + if (modelnum /= modcount) then modcount = modelnum do i=1,ntau n = max0(min0(5,nint(t(i)/2000.)-4),1) @@ -383,15 +383,15 @@ c set up some data upon first entrance with a new model atmosphere endif c initialize some quantities for each new model atmosphere or new frequency - if (modelnum.ne.modcount .or. freq.ne.freq1) then + if (modelnum/=modcount .or. freq/=freq1) then freq1 = freq do n=1,7 - if (freq .gt. freqSi(n)) go to 23 + if (freq > freqSi(n)) go to 23 enddo n = 8 23 dd = (freqlg-flog(n))/(flog(n+1)-flog(n)) - if (n .gt. 2) n = 2*n - 2 - if (n .eq. 14) n = 13 + if (n > 2) n = 2*n - 2 + if (n == 14) n = 13 dd1 = 1.0 - dd do it=1,6 xx(it) = peach(it,n+1)*dd + peach(it,n)*dd1 @@ -399,7 +399,7 @@ c initialize some quantities for each new model atmosphere or new frequency endif do i=1,ntau - if (freq .ge. 7.6869872d14) then + if (freq >= 7.6869872d14) then n = nt(i) aSi2(i) = (dexp(xx(n)*(1.-dt(i)) + xx(n+1)*dt(i))*6.)* . numdens(6,2,i)/u(14,2,i) @@ -443,7 +443,7 @@ c****************************************************************************** data freq1, modcount/0., 0/ c set up some data upon first entrance with a new model atmosphere - if (modelnum .ne. modcount) then + if (modelnum /= modcount) then modcount = modelnum do i=1,ntau hkt = 6.6256d-27/(1.38065d-16*t(i)) @@ -455,20 +455,20 @@ c set up some data upon first entrance with a new model atmosphere c initialize some quantities for each new model atmosphere or new frequency; c the absorption begins at 4762 A. - if (modelnum.ne.modcount .or. freq.ne.freq1) then + if (modelnum/=modcount .or. freq/=freq1) then freq1 = freq waveno = freq/2.99792458d10 - if (waveno .ge. 21000.) then + if (waveno >= 21000.) then do k=1,48 xsect(k) = 0. - if (wno(k) .lt. waveno) xsect(k)= 3.d-18/ + if (wno(k) < waveno) xsect(k)= 3.d-18/ . (1.+((wno(k)+3000.-waveno)/wno(k)/.1)**4) enddo endif endif do i=1,ntau - if (waveno .ge. 21000.) then + if (waveno >= 21000.) then do k=1,48 aFe1(i) = aFe1(i) + xsect(k)*bolt(k,i)* . numdens(7,1,i)/u(26,1,i) @@ -21,7 +21,7 @@ c****************************************************************************** data modcount/0/ c compute scattering, but only if there is a new model atmosphere. - if (modelnum .ne. modcount) then + if (modelnum /= modcount) then modcount = modelnum do i=1,ntau sigel(i) = 0.6653d-24*ne(i) @@ -19,7 +19,7 @@ c****************************************************************************** data newcount, linecount /0, 0/ - if (linecount .eq. 0) oldcount = 0 + if (linecount == 0) oldcount = 0 c IF DOING MULTIPLE RUNS: if this is not the first reading of the @@ -28,7 +28,7 @@ c using "linecount", and then skip the re-initialization of the c various variables rewind nfparam read (nfparam,1001,end=100) arrayz - if (linecount .ne. 0) then + if (linecount /= 0) then do i=1,linecount read (nfparam,1001,end=100) arrayz enddo @@ -210,9 +210,9 @@ c rest of arrayz into "array" c keyword 'RUN' signals that there are either multiple syntheses being c done or multiple comparisons with observed spectra - if (keyword .eq. 'RUN') then + if (keyword == 'RUN') then read (array,*) newcount - if (newcount .gt. oldcount+1) then + if (newcount > oldcount+1) then linecount = linecount - 1 oldcount = syncount go to 100 @@ -229,116 +229,116 @@ c the default value, then the old-style formatted input will be used; c If freeform = 1, unformatted read will be used, BUT the user must then c give values for all quantities (that is, explicit zeros will need to c be put instead of blank spaces. - if (keyword .eq. 'freeform') then + if (keyword == 'freeform') then read (array,*) linfileopt c keyword 'standard_out' controls the name of the verbose standard output - elseif (keyword .eq. 'standard_out') then + elseif (keyword == 'standard_out') then read (array,*) f1out c keyword 'summary_out' controls the name of either the EW summary or c the raw synthesis output - elseif (keyword .eq. 'summary_out') then + elseif (keyword == 'summary_out') then read (array,*) f2out c keyword 'hardpost_out' controls the name of a postscript plot output - elseif (keyword .eq. 'hardpost_out') then + elseif (keyword == 'hardpost_out') then read (array,*) f5out c keyword 'speccomp_out' controls the name of a text file containing the c comparisons (wavelength shifts, sigmas, etc.) between observed and c synthetic spectra - elseif (keyword .eq. 'speccomp_out') then + elseif (keyword == 'speccomp_out') then read (array,*) f6out c keyword 'bin_raw_out' controls the name of a file containing the c raw synthesis of a spectroscopic binary, with an appropriate velocity c difference and luminosity ratio dialed in - elseif (keyword .eq. 'bin_raw_out') then + elseif (keyword == 'bin_raw_out') then read (array,*) f9out c keyword 'bin_smo_out' controls the name of a file containing the c smoothed synthesis of a spectroscopic binary - elseif (keyword .eq. 'bin_smo_out') then + elseif (keyword == 'bin_smo_out') then read (array,*) f10out c keyword 'summary_in' controls the name of the raw synthesis file, c created previously, that will be read in for plotting purposes - elseif (keyword .eq. 'summary_in') then + elseif (keyword == 'summary_in') then read (array,*) f2out c keyword 'smoothed_out' controls the name of the smoothed synthesis output - elseif (keyword .eq. 'smoothed_out') then + elseif (keyword == 'smoothed_out') then read (array,*) f3out c keyword 'keeplines_out' controls the name of the list of kept lines c for future synthetic spectrum runs - elseif (keyword .eq. 'keeplines_out') then + elseif (keyword == 'keeplines_out') then read (array,*) f8out c keyword 'tosslines_out' controls the name of the list of discarded lines c that are too weak to keep in future synthetic spectrum runs - elseif (keyword .eq. 'tosslines_out') then + elseif (keyword == 'tosslines_out') then read (array,*) f9out c keyword 'iraf_out' controls the name of the optional IRAF output - elseif (keyword .eq. 'iraf_out') then + elseif (keyword == 'iraf_out') then read (array,*) f4out c keyword 'model_in' controls the name of input model atmosphere file - elseif (keyword .eq. 'model_in') then + elseif (keyword == 'model_in') then read (array,*) fmodel c keyword 'lines_in' controls the name of the input line list - elseif (keyword .eq. 'lines_in') then + elseif (keyword == 'lines_in') then read (array,*) flines c keyword 'stronglines_in' controls the name of the input strong line list - elseif (keyword .eq. 'stronglines_in') then + elseif (keyword == 'stronglines_in') then read (array,*) fslines c keyword 'observed_in' controls the name of the input observed spectrum - elseif (keyword .eq. 'observed_in') then + elseif (keyword == 'observed_in') then read (array,*) fobs c keyword 'table_in' controls the name of the extra input instruction file - elseif (keyword .eq. 'table_in ') then + elseif (keyword == 'table_in ') then read (array,*) ftable c keyword 'table_out' controls the name of the extra input instruction file - elseif (keyword .eq. 'table_out ') then + elseif (keyword == 'table_out ') then read (array,*) f7out c keyword 'popsyn_out' controls the name of the extra input instruction file - elseif (keyword .eq. 'popsyn_out ') then + elseif (keyword == 'popsyn_out ') then read (array,*) f9out c keyword 'rawbin_out ' controls the name of the input observed spectrum - elseif (keyword .eq. 'rawbin_out ') then + elseif (keyword == 'rawbin_out ') then read (array,*) f9out c keyword 'smoobin_out' controls the name of the input observed spectrum - elseif (keyword .eq. 'smoobin_out') then + elseif (keyword == 'smoobin_out') then read (array,*) f10out @@ -347,7 +347,7 @@ c 0 = do not print out the atmosphere c 1 = print out the standard things about an atmsophere c 2 = print standard things and additional stuff like continuous c opacities, etc. - elseif (keyword .eq. 'atmosphere') then + elseif (keyword == 'atmosphere') then read (array,*) modprintopt @@ -355,13 +355,13 @@ c keyword 'molecules ' controls the molecular equilibrium calculations c 0 = do not do molecular equilibrium c 1 = do molecular equilibrium but do not print results c 2 = do molecular equilibrium and print results - elseif (keyword .eq. 'molecules') then + elseif (keyword == 'molecules') then read (array,*) molopt - if (molopt .eq. 0) then + if (molopt == 0) then nchars = 64 write (array,1009) call getasci (nchars,l0) - if (chinfo(1:1) .eq. 'n') then + if (chinfo(1:1) == 'n') then stop else molopt = 1 @@ -373,7 +373,7 @@ c keyword 'molset' controls the choice of which set of molecules will be c used in molecular equilibrium calculations. c 1 = the small set involving H, C, N, O, Mg, Ti (DEFAULT) c 2 = the large set more useful for very cool stars - elseif (keyword .eq. 'molset') then + elseif (keyword == 'molset') then read (array,*) molset @@ -381,7 +381,7 @@ c keyword 'deviations' controls whether, for synthetic spectrum computations, c an 'obs-comp' plot will be made in addition to the normal spectrum plot c 0 = do not plot the obs-comp plot c 1 = plot the obs-comp plot - elseif (keyword .eq. 'deviations') then + elseif (keyword == 'deviations') then read (array,*) deviations @@ -389,7 +389,7 @@ c keyword 'lines ' controls the output of line data c 0 = print out nothing about the input lines c 1 = print out standard information about the input line list c 2 = gory line data print (usually for diagnostic purposes) - elseif (keyword .eq. 'lines') then + elseif (keyword == 'lines') then read (array,*) linprintopt linprintalt = linprintopt @@ -397,7 +397,7 @@ c 2 = gory line data print (usually for diagnostic purposes) c keyword 'gfstyle ' controls the output of line data c 0 = base-10 logarithms of the gf values (DEFAULT) c 1 = straight gf values - elseif (keyword .eq. 'gfstyle') then + elseif (keyword == 'gfstyle') then read (array,*) gfstyle @@ -405,7 +405,7 @@ c keyword 'contnorm ' allows multiplicative adjustment of the c continuum; useful probably only for batch syntheses c the numbers employed should be around 1.0; c default is 1.000000 - elseif (keyword .eq. 'contnorm') then + elseif (keyword == 'contnorm') then read (array,*) contnorm @@ -416,13 +416,13 @@ c 1 = given in following lines as follows c xlow xhi ylo yhi c vshift lamshift obsadd obsmult c smooth-type FWHM-Gauss vsini L.D.C. FWHM-Macro FWHM-Loren - elseif (keyword .eq. 'plotpars') then + elseif (keyword == 'plotpars') then read (array,*) iscale - if (iscale .ne. 0) then + if (iscale /= 0) then read (nfparam,*) xlo, xhi, ylo, yhi linecount = linecount + 1 read (nfparam,*) veladd, xadd, yadd, ymult - if (xadd .ne. 0.) then + if (xadd /= 0.) then veladd = 3.0d5*xadd/((xlo+xhi)/2.) xadd = 0. endif @@ -436,7 +436,7 @@ c smooth-type FWHM-Gauss vsini L.D.C. FWHM-Macro FWHM-Loren c keyword 'trudamp ' should moog use the detailed line damping for c those transitions that have information stored in c subroutine trudamp? (Default is *no*) - elseif (keyword .eq. 'trudamp') then + elseif (keyword == 'trudamp') then read (array,*) itru @@ -444,7 +444,7 @@ c keyword 'veladjust ' shoud moog try to do a cross-correlation between c observed and synthetic spectra and use that to c align the spectra better in wavelength c (Default is *no*) - elseif (keyword .eq. 'veladjust') then + elseif (keyword == 'veladjust') then read (array,*) maxshift @@ -453,9 +453,9 @@ c outputs the final spectrum c 0 = angs c 1 = microns c 2 = 1/cm - elseif (keyword .eq. 'units') then + elseif (keyword == 'units') then read (array,*) iunits - if (iunits .ne. 0) then + if (iunits /= 0) then write (*,1010) stop endif @@ -465,7 +465,7 @@ c keyword 'iraf ' allows the user to output a raw spectrum in c a form suitable for IRAF's rtext input command c 0 = don't do this, make output the normal way. c 1 = make an IRAF-compatible output - elseif (keyword .eq. 'iraf') then + elseif (keyword == 'iraf') then read (array,*) iraf @@ -473,14 +473,14 @@ c keyword 'scat 'allows the user to employ a source function c which has both scattering and absorption components c 0 = NO scattering c 1 = scattering - elseif (keyword .eq. 'scat') then + elseif (keyword == 'scat') then read (array,*) scatopt c keyword 'flux/int ' choses integrated flux or central intensity c 0 = integrated flux calculations c 1 = central intensity calculations - elseif (keyword .eq. 'flux/int') then + elseif (keyword == 'flux/int') then read (array,*) fluxintopt @@ -505,7 +505,7 @@ c dampingopt = 3 and dampnum > 10^(-10) ---> c c6 = (c6_NEXTGEN for H I, He I, H2)*dampnum c for molecular lines (lacking a better idea) ---> c c6 done as in dampingopt = 0 - elseif (keyword .eq. 'damping') then + elseif (keyword == 'damping') then read (array,*) dampingopt @@ -517,9 +517,9 @@ c 2 = (not implemented yet) c 3 = read a true Fits file with the FITSIO package c 4 = (not implemented yet) c 5 = read a special MONGO style (wavelength, flux pair) file - elseif (keyword .eq. 'obspectrum') then + elseif (keyword == 'obspectrum') then read (array,*) specfileopt - if (specfileopt .lt. 0) then + if (specfileopt < 0) then byteswap = 1 specfileopt = iabs(specfileopt) endif @@ -527,14 +527,14 @@ c 5 = read a special MONGO style (wavelength, flux pair) file c keyword 'histogram' makes histogram plots of observed spectra if c histoyes = 1 - elseif (keyword .eq. 'histogram') then + elseif (keyword == 'histogram') then read (array,*) histoyes c keyword 'terminal ' gives the sm plotting window type c smterm = a character string of the sm window type (see the c appropriate sm manual for a list) - elseif (keyword .eq. 'terminal') then + elseif (keyword == 'terminal') then read (array,*) smterm @@ -547,7 +547,7 @@ c For line analyses: # = the minimum number of lines of a c species necessary to trigger a plot c For curves-of-growth: 1 = make plots c For flux curves: 1 = make plots - elseif (keyword .eq. 'plot') then + elseif (keyword == 'plot') then read (array,*) plotopt @@ -556,7 +556,7 @@ c # = the number of different syntheses to run c (the next line gives the different abundance factors c to use) c minimum error check: numatomsyn must equal numisosyn or code will stop - elseif (keyword .eq. 'abundances') then + elseif (keyword == 'abundances') then neq = 0 numpecatom = 0 numatomsyn = 0 @@ -573,8 +573,8 @@ c minimum error check: numatomsyn must equal numisosyn or code will stop enddo enddo read (array,*) numpecatom,numatomsyn - if (numisosyn .ne. 0) then - if (numatomsyn .ne. numisosyn) then + if (numisosyn /= 0) then + if (numatomsyn /= numisosyn) then write (array,1002) numatomsyn, numisosyn call putasci (77,6) stop @@ -583,7 +583,7 @@ c minimum error check: numatomsyn must equal numisosyn or code will stop do l=1,numpecatom read (nfparam,*) jatom,(deltalogab(kk),kk=1,numatomsyn) linecount = linecount + 1 - if (jatom .eq. 99) then + if (jatom == 99) then do kk=1,numatomsyn abfactor (kk) = deltalogab(kk) enddo @@ -594,13 +594,13 @@ c minimum error check: numatomsyn must equal numisosyn or code will stop pec(jatom) = 1 endif enddo - if (numpecatom.eq.1 .and. jatom.eq.99) ninetynineflag = 1 + if (numpecatom==1 .and. jatom==99) ninetynineflag = 1 c keyword 'isotopes ' gives the isotopes used in the line list and their c abundance relative to the parent spiecies c minimum error check: numatomsyn must equal numisosyn or code will stop - elseif (keyword .eq. 'isotopes') then + elseif (keyword == 'isotopes') then numiso = 0 numisosyn = 0 newnumiso = 0 @@ -614,8 +614,8 @@ c minimum error check: numatomsyn must equal numisosyn or code will stop enddo enddo read (array,*) numiso,numisosyn - if (numatomsyn .ne. 0) then - if (numatomsyn .ne. numisosyn) then + if (numatomsyn /= 0) then + if (numatomsyn /= numisosyn) then write (array,1002) numatomsyn, numisosyn call putasci (77,6) stop @@ -630,13 +630,13 @@ c minimum error check: numatomsyn must equal numisosyn or code will stop c keyword 'lumratio' gives the ratio of the luminosity of two stars at a c specific wavelength in a binary star system (used c only with driver "binary") - elseif (keyword .eq. 'lumratio') then + elseif (keyword == 'lumratio') then read (array,*) lumratio c keyword 'deltaradvel' gives the velocity difference between the stars c binary star system (used only with driver "binary") - elseif (keyword .eq. 'deltaradvel') then + elseif (keyword == 'deltaradvel') then read (array,*) deltaradvel @@ -646,14 +646,14 @@ c wavelengths, step is the step size in the c syntheses, and delta is the wavelength range c to either side of a synthesis point to consider c for line opacity calculations - elseif (keyword .eq. 'synlimits') then + elseif (keyword == 'synlimits') then read (nfparam,*) start, sstop, step, delta oldstart = start oldstop = sstop oldstep = step olddelta = delta step1000 = 1000.*step - if (dble(idnint(step1000))-step1000 .ne. 0.) then + if (dble(idnint(step1000))-step1000 /= 0.) then write (*,1008) step stop endif @@ -664,7 +664,7 @@ c keyword 'fluxlimits' gives the wavelength parameters for flux curves; c start and sstop are beginning and ending c wavelengths, and step is the step size in the c flux curve - elseif (keyword .eq. 'fluxlimits') then + elseif (keyword == 'fluxlimits') then read (nfparam,*) start, sstop, step linecount = linecount + 1 @@ -677,7 +677,7 @@ c step is the wavelength step size in the c computations; cogatom is the name of the c element whose abundance should be varied c to achieve an EW match with observations. - elseif (keyword .eq. 'blenlimits') then + elseif (keyword == 'blenlimits') then read (nfparam,*) delwave, step, cogatom linecount = linecount + 1 @@ -692,7 +692,7 @@ c for spectrum synthesis curves-of-growth, c and wavestep is a forced (if desired) step size c in wavelength along the line (this applies to c single line computations only - elseif (keyword .eq. 'coglimits') then + elseif (keyword == 'coglimits') then read (nfparam,*) rwlow, rwhigh, rwstep, wavestep, cogatom linecount = linecount + 1 @@ -700,7 +700,7 @@ c single line computations only c keyword 'limits ' old limits format...tell the user to change the c keyword and quit. - elseif (keyword .eq. 'limits') then + elseif (keyword == 'limits') then write(*,*) 'Warning: keyword changed to *synlimits*, *coglimits*' write(*,*) 'for Syntesis and COG calculations.' write(*,*) 'Here are the proper formats:' @@ -713,7 +713,7 @@ c keyword and quit. c keyword of strong for lines which are to be considered for all of the c synthesis - elseif (keyword .eq. 'strong') then + elseif (keyword == 'strong') then read (array,*) dostrong @@ -721,7 +721,7 @@ c keyword word of opacit which takes the continuus opacity and scales it c with the form of kaplam(i)= kaplam(i)*((factor*10000)/t(i)) c in Opacit.f after it calulates the normal kaplam c if value is <= 0 then it does not do it - elseif (keyword .eq. 'opacit') then + elseif (keyword == 'opacit') then read (array,*) fudge @@ -739,33 +739,33 @@ c loop back to get another parameter c wrap things up with a few assignments -98 if (control.eq.'gridsyn' .or. control.eq.'gridplo' .or. - . control.eq.'binary ' .or. control.eq.'abandy ') then +98 if (control=='gridsyn' .or. control=='gridplo' .or. + . control=='binary ' .or. control=='abandy ') then control = 'gridend' endif c assign plotting window type; if no type has been given in the c parameter file, then ask for it -100 if (smterm .eq. ' ') then +100 if (smterm == ' ') then array = 'GIVE THE SM TERMINAL NAME : ' nchar = 28 call getasci (nchar,12) smterm = chinfo(1:nchar) ivstat = ivcleof(12,1) endif - if (smterm.eq.'x11' .or. smterm.eq.'X11') then - if (control .eq. 'synth ' .or. - . control .eq. 'synpop ' .or. - . control .eq. 'synplot' .or. - . control .eq. 'isoplot' .or. - . control .eq. 'gridsyn' .or. - . control .eq. 'gridplo' .or. - . control .eq. 'doflux ' .or. - . control .eq. 'cogsyn ' .or. - . control .eq. 'cog ' .or. - . control .eq. 'isotop ' .or. - . control .eq. 'binary ') then + if (smterm=='x11' .or. smterm=='X11') then + if (control == 'synth ' .or. + . control == 'synpop ' .or. + . control == 'synplot' .or. + . control == 'isoplot' .or. + . control == 'gridsyn' .or. + . control == 'gridplo' .or. + . control == 'doflux ' .or. + . control == 'cogsyn ' .or. + . control == 'cog ' .or. + . control == 'isotop ' .or. + . control == 'binary ') then smterm = smt1 else smterm = smt2 @@ -774,14 +774,14 @@ c parameter file, then ask for it c for syntheses, store the plotting parameters - if (control.eq.'synth ' .or. control.eq.'synplot' .or. - . control.eq.'gridsyn' .or. control.eq.'gridplo' .or. - . control.eq.'binary ' .or. control.eq.'synpop ') then - if (oldstart .eq. 0) then + if (control=='synth ' .or. control=='synplot' .or. + . control=='gridsyn' .or. control=='gridplo' .or. + . control=='binary ' .or. control=='synpop ') then + if (oldstart == 0) then write (*,1011) stop endif - if (iscale .eq. 0) call plotremember (0) + if (iscale == 0) call plotremember (0) call plotremember (1) endif @@ -22,7 +22,7 @@ c*****compute partition functions for 4 ionization states of an element. do k=1,4 iat = iat + 1 at = dfloat(iat)/10. - if (partflag(iatom,k) .gt. 0) then + if (partflag(iatom,k) > 0) then do i=1,ntau u(jmark,k,i) = partnew(at,k,i) enddo @@ -21,7 +21,7 @@ c****************************************************************************** iatom = nint(atom) iarray = partflag(iatom,k) - if (level .gt. 500) then + if (level > 500) then temp = dlog(dble(level)) else temp = tlog(level) @@ -21,7 +21,7 @@ c spectra (output), and (if desired) IRAF-style smoothed spectra (output) nchars = 19 call infile ('input ',nf2out,'formatted ',0,nchars, . f2out,lscreen) - if (plotopt .ne. 0) then + if (plotopt /= 0) then nf3out = 22 lscreen = lscreen + 2 array = 'SMOOTHED SYNTHESES OUTPUT' @@ -35,7 +35,7 @@ c spectra (output), and (if desired) IRAF-style smoothed spectra (output) call infile ('output ',nf5out,'formatted ',0,nchars, . f5out,lscreen) endif - if (plotopt .gt. 1) then + if (plotopt > 1) then nf6out = 27 lscreen = lscreen + 2 array = 'SPECTRUM COMPARISON OUTPUT' @@ -43,7 +43,7 @@ c spectra (output), and (if desired) IRAF-style smoothed spectra (output) call infile ('output ',nf6out,'formatted ',0,nchars, . f6out,lscreen) endif - if (iraf .ne. 0) then + if (iraf /= 0) then nf4out = 23 lscreen = lscreen + 2 array = 'IRAF ("rtext") OUTPUT' @@ -54,12 +54,12 @@ c spectra (output), and (if desired) IRAF-style smoothed spectra (output) c*****now plot the spectrum - if (plotopt.eq.2 .and. specfileopt.gt.0) then + if (plotopt==2 .and. specfileopt>0) then nfobs = 33 lscreen = lscreen + 2 array = 'THE OBSERVED SPECTRUM' nchars = 21 - if (specfileopt.eq.1 .or. specfileopt.eq.3) then + if (specfileopt==1 .or. specfileopt==3) then call infile ('input ',nfobs,'unformatted',2880,nchars, . fobs,lscreen) else @@ -67,7 +67,7 @@ c*****now plot the spectrum . fobs,lscreen) endif endif - if (plotopt .ne. 0) then + if (plotopt /= 0) then line = 10 ncall = 1 call pltspec (line,ncall) diff --git a/Plotremember.f b/Plotremember.f index 5f1619f..dc346a9 100755 --- a/Plotremember.f +++ b/Plotremember.f @@ -14,8 +14,8 @@ c****************************************************************************** c*****initialize the plot parameters; iscale=0 is simple default; c iscale=1 is when these parameters have been read from the c parameter file - if (option .eq. 0) then - if (iscale .eq. 0) then + if (option == 0) then + if (iscale == 0) then xlo = oldstart xhi = oldstop ylo = 0. @@ -35,7 +35,7 @@ c parameter file c*****store the original plot parameters - elseif (option .eq. 1) then + elseif (option == 1) then origxlo = xlo origxhi = xhi origylo = ylo @@ -54,7 +54,7 @@ c*****store the original plot parameters c*****re-set the plot parameters to their original values - elseif (option .eq. 2) then + elseif (option == 2) then xlo = origxlo xhi = origxhi ylo = origylo @@ -73,7 +73,7 @@ c*****re-set the plot parameters to their original values c*****store the plot parameters from the last entry into pltspec - elseif (option .eq. 3) then + elseif (option == 3) then oldxlo = xlo oldxhi = xhi oldylo = ylo @@ -92,7 +92,7 @@ c*****store the plot parameters from the last entry into pltspec c*****re-set plot parameters to values from the last entry into pltspec - elseif (option .eq. 4) then + elseif (option == 4) then xlo = oldxlo xhi = oldxhi ylo = oldylo @@ -14,8 +14,8 @@ c****************************************************************************** real*4 shortnum - if (kount .ge. plotopt) then - if (plotopt .eq. 0) return + if (kount >= plotopt) then + if (plotopt == 0) return 10 choice = 'y' plotroutine = 'term_port_abun' lscreen = maxline -2 @@ -25,11 +25,11 @@ c****************************************************************************** nchars = 35 call getasci (nchars,maxline) choice = chinfo(1:1) - if (choice.eq.'n' .or. nchars.le.0) then + if (choice=='n' .or. nchars<=0) then return - elseif (choice .eq. 'm') then + elseif (choice == 'm') then return - elseif (choice .eq. 'v') then + elseif (choice == 'v') then write (array,*) 'What is the new microturbulence (km/s)? ' nchars = 41 call getnum (nchars,lscreen,xnum,shortnum) @@ -42,12 +42,12 @@ c****************************************************************************** lim1obs = 0 lim2obs = 0 return - elseif (choice .eq. 'h') then + elseif (choice == 'h') then plotroutine = 'hard_port_abun' call makeplot (lscreen) - elseif (choice .eq. 'r') then + elseif (choice == 'r') then go to 10 - elseif (choice .eq. 'f') then + elseif (choice == 'f') then plotroutine = 'file_port_abun' call makeplot (lscreen) endif @@ -13,7 +13,7 @@ c****************************************************************************** c call up the curve-of-growth plot - if (plotopt .eq. 0) return + if (plotopt == 0) return 10 choice = 'y' plotroutine = 'term_land_cog ' lscreen = 12 @@ -28,13 +28,13 @@ c screen, try a new model atmosphere, or replot? nchars = 37 call getasci (nchars,lscreen) choice = chinfo(1:1) - if (choice.eq.'n' .or. nchars.le.0) then + if (choice=='n' .or. nchars<=0) then return - elseif (choice .eq. 'h') then + elseif (choice == 'h') then plotroutine = 'hard_land_cog ' call makeplot (lscreen) go to 10 - elseif (choice .eq. 'v') then + elseif (choice == 'v') then write (array,*) 'What is the new microturbulence (km/s)? ' nchars = 41 lscreen = lscreen + 2 @@ -48,17 +48,17 @@ c screen, try a new model atmosphere, or replot? rewind nfmodel rewind nflines return - elseif (choice .eq. 'm') then + elseif (choice == 'm') then return - elseif (choice .eq. 'r') then + elseif (choice == 'r') then go to 10 - elseif (choice .eq. 'p') then + elseif (choice == 'p') then array = 'MARK THE POSITION WITH THE CURSOR' istat=ivcleof(21,1) istat=ivwrite(13,3,array,34) call drawcurs go to 1 - elseif (choice .eq. 'f') then + elseif (choice == 'f') then plotroutine = 'file_land_cog ' call makeplot (lscreen) go to 10 @@ -12,7 +12,7 @@ c****************************************************************************** c call up the flux plot - if (plotopt .eq. 0) return + if (plotopt == 0) return 10 choice = 'y' plotroutine = 'term_land_flux' lscreen = 12 @@ -25,14 +25,14 @@ c make a hardcopy, write to a postscript file, or replot? nchars = 33 call getasci (nchars,lscreen) choice = chinfo(1:1) - if (choice.eq.'n' .or. nchars.le.0) then + if (choice=='n' .or. nchars<=0) then return - elseif (choice .eq. 'h') then + elseif (choice == 'h') then plotroutine = 'hard_land_flux' call makeplot (lscreen) - elseif (choice .eq. 'r') then + elseif (choice == 'r') then go to 10 - elseif (choice .eq. 'f') then + elseif (choice == 'f') then plotroutine = 'file_land_flux' call makeplot (lscreen) endif @@ -17,8 +17,8 @@ c****************************************************************************** c*****initialize some variables, or re-set them to old values -3 if (ncall .eq. 1) then - if (iscale .eq. 0) then +3 if (ncall == 1) then + if (iscale == 0) then choice = 's' else choice = '1' @@ -30,16 +30,16 @@ c*****initialize some variables, or re-set them to old values c*****make a special "choice" for grid syntheses, which go directly to c*****postscript files if desired - if (control.eq.'gridsyn' .or. control.eq.'gridend' .or. - . control.eq.'gridplo') then + if (control=='gridsyn' .or. control=='gridend' .or. + . control=='gridplo') then choice = 'g' endif c*****begin by reading in an observed spectrum - if (specfileopt.ge.1 .and. plotopt.eq.2 .and. ncall.eq.1) then + if (specfileopt>=1 .and. plotopt==2 .and. ncall==1) then call readobs (line) - if (lount .eq. -1) then + if (lount == -1) then array = 'OBSERVED SPECTRUM FILE PROBLEM! I QUIT.' istat = ivwrite (line+4,3,array,40) stop @@ -49,22 +49,22 @@ c*****begin by reading in an observed spectrum c*****go through the option list; the routine may be exited at this point istat = ivcleof(4,1) -1 if (choice .eq. 'q') return +1 if (choice == 'q') return c*****or a default plot may be made upon entering the routine - if (choice.eq.'1' .or. choice.eq.'g') then + if (choice=='1' .or. choice=='g') then c*****make a cross correlation to line up synthetic and observed spectra c*****in velocity (wavelength) space; the user can turn this off/on c*****not working for you call smooth (-1,ncall) c---------------------------------------------------------------------------- - if (plotopt.eq.2 .and. ncall.eq.1) then + if (plotopt==2 .and. ncall==1) then vfactor = 1.0 + veladd/2.99795e+5 do i=1,lount xobs(i) = vfactor*xobs(i) enddo - if (maxshift .gt. 0) then + if (maxshift > 0) then call correl (maxshift) vfactor = 1.0 + deltavel/2.99795e+5 do i=1,lount @@ -73,23 +73,23 @@ c---------------------------------------------------------------------------- endif endif c---------------------------------------------------------------------------- - if (ncall .eq. 1) then + if (ncall == 1) then wmiddle = (start + sstop)/2. - if (iunits .eq. 1) wmiddle = 1.d-4*wmiddle - if (ymult .eq. 0.0) ymult = 1.0 + if (iunits == 1) wmiddle = 1.d-4*wmiddle + if (ymult == 0.0) ymult = 1.0 do i=1,lount yobs(i) = ymult*yobs(i) yobs(i) = yadd+yobs(i) enddo - if (xlo .eq. xhi) then + if (xlo == xhi) then xlo = start xhi = sstop - if (iunits .eq. 1) then + if (iunits == 1) then xlo = 1.d-4*xlo xhi = 1.d-4*xhi endif endif - if (ylo .eq. yhi) then + if (ylo == yhi) then ylo = 0. yhi = 1.1 endif @@ -101,24 +101,24 @@ c---------------------------------------------------------------------------- c*****or the synthetic spectra may be resmoothed; if a problem develops in c a user-specified parameter (like a Gaussian FWHM that is too large), c then output a warning and let user decide what to do next - if (choice .eq. 's') then + if (choice == 's') then call smooth (line+2,1) -2 if (smtype .eq. 'e') then +2 if (smtype == 'e') then array = 'REDO THE SMOOTHING (y/n)? ' nchars = 26 call getasci (nchars,line+9) smtype = chinfo(1:1) - if (smtype .eq. 'n') then + if (smtype == 'n') then go to 100 else istat = ivcleof (10,1) go to 2 endif endif - if (xlo .eq. 0.0 .and. xhi .eq. 0.0) then + if (xlo == 0.0 .and. xhi == 0.0) then xlo = start xhi = sstop - if (iunits .eq. 1) then + if (iunits == 1) then xlo = 1.d-4*xlo xhi = 1.d-4*xhi endif @@ -127,7 +127,7 @@ c then output a warning and let user decide what to do next c*****or the observations may be rescaled - if (choice .eq. 'r') then + if (choice == 'r') then array = 'MULTIPLY THE OBSERVED POINTS BY WHAT FACTOR? ' nchars = 45 call getnum (nchars,13,xnum,yymult) @@ -139,7 +139,7 @@ c*****or the observations may be rescaled c*****or the observations may be shifted by an additive constant - if (choice .eq. 'a') then + if (choice == 'a') then array = 'ADD WHAT NUMBER TO THE OBSERVED POINTS? ' nchars = 40 call getnum (nchars,13,xnum,yyadd) @@ -151,7 +151,7 @@ c*****or the observations may be shifted by an additive constant c*****or the observations may be shifted by a constant wavelength - if (choice .eq. 'w') then + if (choice == 'w') then array = 'SHIFT THE OBSERVED POINTS BY WHAT WAVELENGTH? ' nchars = 46 call getnum (nchars,13,xnum,xxadd) @@ -165,7 +165,7 @@ c*****or the observations may be shifted by a constant wavelength c*****or the observations may be shifted by a constant velocity - if (choice .eq. 'v') then + if (choice == 'v') then array = 'SHIFT THE OBSERVED POINTS BY WHAT VELOCITY (KM/S)? ' nchars = 51 call getnum (nchars,13,xnum,vveladd) @@ -178,7 +178,7 @@ c*****or the observations may be shifted by a constant velocity c*****or the plot boundaries may be changed - if (choice .eq. 'c') then + if (choice == 'c') then write (array,1001) xlo nchars = 29 call getnum (nchars,15,xnum,xlo) @@ -196,7 +196,7 @@ c*****or the plot boundaries may be changed c*****or the cross hairs can be used to zoom in on a part of the plot - if (choice .eq. 'z') then + if (choice == 'z') then array = 'MARK THE LOWER LEFT HAND CORNER WITH THE CURSOR' istat = ivcleof(13,1) istat = ivwrite(13,3,array,47) @@ -210,7 +210,7 @@ c*****or the cross hairs can be used to zoom in on a part of the plot xhi = xplotpos yhi = yplotpos call boxit - if (iunits .eq. 1) then + if (iunits == 1) then xlo = 1.d-4*xlo xhi = 1.d-4*xhi endif @@ -220,7 +220,7 @@ c*****or the cross hairs can be used to zoom in on a part of the plot c*****or cursor position can be returned - if (choice .eq. 'p') then + if (choice == 'p') then array = 'MARK THE POSITION WITH THE CURSOR' istat=ivcleof(21,1) istat=ivwrite(13,3,array,34) @@ -230,7 +230,7 @@ c*****or cursor position can be returned c*****or the title of the model can be changed - if (choice .eq. 't') then + if (choice == 't') then array = 'ENTER THE NEW TITLE' istat = ivcleof(21,1) istat = ivwrite (13,3,array,19) @@ -240,7 +240,7 @@ c*****or the title of the model can be changed c*****or the spectra can be replotted, with a separate plot showing the c observed/synthtic spectrum differences - if (choice .eq. 'd') then + if (choice == 'd') then deviations = 1 whichwin = '2of2' endif @@ -248,10 +248,10 @@ c observed/synthtic spectrum differences c*****or the plot boundaries may be reset to the original values; c this is a basic starting over plot - if (choice .eq. 'o') then + if (choice == 'o') then xlo = start xhi = sstop - if (iunits .eq. 1) then + if (iunits == 1) then xlo = 1.d-4*xlo xhi = 1.d-4*xhi endif @@ -267,19 +267,19 @@ c this is a basic starting over plot c*****or the plot can simply be redone - if (choice .eq. 'm') then + if (choice == 'm') then go to 90 endif c*****now either here make a hardcopy plot - if (choice .eq. 'h') then - if (control .eq. 'binary ') then + if (choice == 'h') then + if (control == 'binary ') then plotroutine = 'hard_land_bin ' else plotroutine = 'hard_land_spec' endif - if (deviations .eq. 0) then + if (deviations == 0) then whichwin = '1of1' else whichwin = '2of2' @@ -291,20 +291,20 @@ c*****now either here make a hardcopy plot c*****or write the plot to a postscript file - if (choice.eq.'f' .or. choice.eq.'g') then - if (control .eq. 'binary ') then + if (choice=='f' .or. choice=='g') then + if (control == 'binary ') then plotroutine = 'file_land_bin ' else plotroutine = 'file_land_spec' endif - if (deviations .eq. 0) then + if (deviations == 0) then whichwin = '1of1' else whichwin = '2of2' endif lscreen = 12 call makeplot (lscreen) - if (choice .eq. 'g') then + if (choice == 'g') then return else go to 90 @@ -313,7 +313,7 @@ c*****or write the plot to a postscript file c*****or return to the calling routine in order to change abundances - if (choice .eq. 'n') then + if (choice == 'n') then call plotremember (3) return endif @@ -323,7 +323,7 @@ c*****or add an additional uniform amount of flux, expressed in terms of c the current continuum flux; this only is approximately physically c correct if spectrograph smoothing is negligible compared to other c smoothing - if (choice .eq. 'l') then + if (choice == 'l') then write (array,*) . 'WHAT IS THE ADDITIONAL FLUX IN TERMS OF CONTINUUM [0.0]? ' nchars = 57 @@ -334,7 +334,7 @@ c smoothing c*****or if total confusion has happened in the plotting, reset all parameters c to their original values, and replot - if (choice .eq. 'u') then + if (choice == 'u') then call plotremember (2) ncall = 1 go to 3 @@ -342,16 +342,16 @@ c to their original values, and replot c*****or plot on the terminal -90 if (control.eq.'gridsyn' .or. control.eq.'gridend' .or. - . control .eq. 'gridplo') return - if (control .eq. 'binary ') then +90 if (control=='gridsyn' .or. control=='gridend' .or. + . control == 'gridplo') return + if (control == 'binary ') then plotroutine = 'term_land_bin ' else plotroutine = 'term_land_spec' endif lscreen = 12 - if (choice.eq.'f' .or. choice.eq.'h') choice = 'm' - if (deviations .eq. 0) then + if (choice=='f' .or. choice=='h') choice = 'm' + if (deviations == 0) then whichwin = '1of1' else whichwin = '2of2' @@ -387,15 +387,15 @@ c*****finally, print the option table c*****reprint the option table if the choice is not understood c or take action on the choice - if (choice.eq.'s' .or. choice.eq.'r' .or. - . choice.eq.'a' .or. choice.eq.'h' .or. - . choice.eq.'c' .or. choice.eq.'q' .or. - . choice.eq.'m' .or. choice.eq.'o' .or. - . choice.eq.'v' .or. choice.eq.'w' .or. - . choice.eq.'z' .or. choice.eq.'p' .or. - . choice.eq.'t' .or. choice.eq.'f' .or. - . choice.eq.'n' .or. choice.eq.'d' .or. - . choice.eq.'l' .or. choice.eq.'u') then + if (choice=='s' .or. choice=='r' .or. + . choice=='a' .or. choice=='h' .or. + . choice=='c' .or. choice=='q' .or. + . choice=='m' .or. choice=='o' .or. + . choice=='v' .or. choice=='w' .or. + . choice=='z' .or. choice=='p' .or. + . choice=='t' .or. choice=='f' .or. + . choice=='n' .or. choice=='d' .or. + . choice=='l' .or. choice=='u') then go to 1 else go to 100 diff --git a/Pointcurs.f b/Pointcurs.f index d4bb692..fa46709 100755 --- a/Pointcurs.f +++ b/Pointcurs.f @@ -9,9 +9,9 @@ c****************************************************************************** call sm_graphics - if (whichwin .eq. '1of1') then + if (whichwin == '1of1') then call sm_window (1,1,1,1,1,1) - elseif (whichwin .eq. '2of2') then + elseif (whichwin == '2of2') then call sm_defvar ('y_gutter','0.0') call sm_window (1,2,1,1,1,1) endif @@ -11,25 +11,25 @@ c****************************************************************************** c*****do not print out information in real time if the code is in c batch mode - if (silent .eq. 'y') return + if (silent == 'y') return - if (line .eq. 1) then + if (line == 1) then istat = ivcleof(4,1) endif - if (line .eq. maxline-5) then + if (line == maxline-5) then errm1 = errmess array1 = array 10 array = 'WANT TO SEE MORE ([y]/n)? ' nchars = 26 call getasci (nchars,4+line) - if (chinfo(1:1).eq.'y' .or. nchars.le.0) then + if (chinfo(1:1)=='y' .or. nchars<=0) then istat = ivcleof(4,1) line = 1 array = array1 errmess = errm1 - elseif (chinfo(1:1) .eq. 'n') then + elseif (chinfo(1:1) == 'n') then errmess = 'stopinfo!' return else @@ -42,8 +42,8 @@ c first get the header records and search for key parameters 100 irec = 1 101 read (unit=nfobs,rec=irec,err=1002,iostat=ierr) head call obshead (head,iend,line) - if (lount .eq. -1) return - if (iend .eq. 0) then + if (lount == -1) return + if (iend == 0) then irec = irec + 1 go to 101 endif @@ -52,13 +52,13 @@ c first get the header records and search for key parameters c next read the flux array from the file nrec = lount/nblock ipt = 0 - if (mod(lount,nblock) .ne. 0) nrec = nrec + 1 + if (mod(lount,nblock) /= 0) nrec = nrec + 1 do j=1,nrec irec = irec + 1 jpt = min0(nblock,lount-ipt) - if (ibits .eq. 16) then + if (ibits == 16) then read (unit=nfobs,rec=irec,err=1006,iostat=ierr) int1 - if (byteswap .eq. 1) then + if (byteswap == 1) then do k=2,2880,2 onebyte = int1(k) int1(k) = int1(k-1) @@ -68,9 +68,9 @@ c next read the flux array from the file do k=1,jpt yobs(ipt+k) = bzero + bscale*real(int2(k)) enddo - elseif (ibits .eq. 32) then + elseif (ibits == 32) then read (unit=nfobs,rec=irec,err=1006,iostat=ierr) int1 - if (byteswap .eq. 1) then + if (byteswap == 1) then do k=4,2880,4 onebyte = int1(k) int1(k) = int1(k-3) @@ -83,9 +83,9 @@ c next read the flux array from the file do k=1,jpt yobs(ipt+k) = bzero + bscale*real(int4(k)) enddo - elseif(ibits .eq. -32) then + elseif(ibits == -32) then read (unit=nfobs,rec=irec,err=1006,iostat=ierr) int1 - if (byteswap .eq. 1) then + if (byteswap == 1) then do k=4,2880,4 onebyte = int1(k) int1(k) = int1(k-3) @@ -115,7 +115,7 @@ c now fill in the wavelength array enddo lim1obs = 1 lim2obs = lount - if(xobs(2) .lt. xobs(1)) then + if(xobs(2) < xobs(1)) then do j=1,lount/2 xtemp = xobs(j) ytemp = yobs(j) @@ -160,7 +160,7 @@ c*****here is a MONGO-style input array 525 lount = i - 1 lim1obs = 1 lim2obs = lount - if(xobs(2) .lt. xobs(1)) then + if(xobs(2) < xobs(1)) then do j=1,lount/2 xtemp = xobs(j) ytemp = yobs(j) @@ -35,7 +35,7 @@ c****************************************************************************** c(n)=0. b(n)=(f(n)-f(n1))/(x(n)-x(n1)) a(n)=f(n)-x(n)*b(n) - if(n.eq.2)return + if(n==2)return do 1 j=2,n1 j1=j-1 d=(f(j)-f(j1))/(x(j)-x(j1)) @@ -49,7 +49,7 @@ c****************************************************************************** c(3)=0. b(3)=(f(4)-f(3))/(x(4)-x(3)) a(3)=f(3)-x(3)*b(3) - if(c(j).eq.0.)go to 2 + if(c(j)==0.)go to 2 j1=j+1 wt=abs(c(j1))/(abs(c(j1))+abs(c(j))) a(j)=a(j1)+wt*(a(j)-a(j1)) @@ -24,7 +24,7 @@ c with the H_2 molecule c first do the neutrals do j=1,7 do k=1,neq - if (nel(j) .eq. iorder(k)) then + if (nel(j) == iorder(k)) then do i=1,ntau numdens(j,1,i) = xamol(k,i) enddo @@ -39,7 +39,7 @@ c*****then do the ions ispec10 = nint(dble(10*nel(j)+1)) do k=1,nmol kmol10 = nint(10.*amol(k)) - if (ispec10 .eq. kmol10) then + if (ispec10 == kmol10) then do i=1,ntau numdens(j,2,i) = xmol(k,i) enddo @@ -52,7 +52,7 @@ c*****then do the ions c*****finally add in H_2 do k=1,nmol ispec = 101 - if (ispec .eq. nint(amol(k))) then + if (ispec == nint(amol(k))) then do i=1,ntau numdens(8,1,i) = xmol(k,i) enddo @@ -65,7 +65,7 @@ c*****compute partitiion functions for H_2O and CO_2; do i=1,ntau h2olog = 0. co2log = 0. - if (t(i) .gt. 5000.) then + if (t(i) > 5000.) then uh2o(i) = 1.d8 uco2(i) = 1.d8 else @@ -83,8 +83,8 @@ c*****transfer H_2O and CO_2 number densities from the molecular c equilibrium output c note: HITRAN partition functions are given only for T < 5000K do j=1,nmol - if (nint(amol(j)) .eq. 10108) ih2o = j - if (nint(amol(j)) .eq. 60808) ico2 = j + if (nint(amol(j)) == 10108) ih2o = j + if (nint(amol(j)) == 60808) ico2 = j enddo do i=1,ntau xnh2o(i) = xmol(ih2o,i) @@ -19,18 +19,18 @@ c*****initialize parameters write (abitle(1:400),1081) write (isoitle(1:240),1082) nsyn = 1 - if (ncall .eq. 1) then + if (ncall == 1) then gaussflag = 'f' rotateflag = 'f' lorenflag = 'f' macroflag = 'f' - if (choice .ne. 'l') addflux = 0. + if (choice /= 'l') addflux = 0. endif c*****on entering, figure out what kind of smoothing is desired, unless c the default smoothing options have been set for first pass - if (line .gt. 0) then + if (line > 0) then 2 write (array,1007) istat = ivwrite (line,3,array,67) write (array,1004) @@ -39,17 +39,17 @@ c the default smoothing options have been set for first pass nchars = 21 call getasci (nchars,line+2) smtype = chinfo(1:1) - if (smtype.ne.'n' .and. smtype.ne.'g' .and. - . smtype.ne.'l' .and. smtype.ne.'v' .and. - . smtype.ne.'m' .and. smtype.ne.'c' .and. - . smtype.ne.'d' .and. smtype.ne.'r' .and. - . smtype.ne.'p') go to 2 + if (smtype/='n' .and. smtype/='g' .and. + . smtype/='l' .and. smtype/='v' .and. + . smtype/='m' .and. smtype/='c' .and. + . smtype/='d' .and. smtype/='r' .and. + . smtype/='p') go to 2 endif c if a user-specified variable Gaussian smoothing over the spectral range c is called for, option 'p', branch to a different routine - if (smtype .eq. 'p') then + if (smtype == 'p') then call vargauss (line+1) return endif @@ -61,12 +61,12 @@ c files, and get the synthesis range parameters from the 'dump' file rewind nf3out do i=1,20 read (nf2out,1002) array - if (array(1:7).eq.'element' .or. - . array(1:7).eq.'Changin' .or. - . array(1:7).eq.'ALL abu' .or. - . array(1:7).eq.'Isotopi') then + if (array(1:7)=='element' .or. + . array(1:7)=='Changin' .or. + . array(1:7)=='ALL abu' .or. + . array(1:7)=='Isotopi') then cycle - elseif (array(1:7).eq.'MODEL: ') then + elseif (array(1:7)=='MODEL: ') then moditle(1:73) = array(8:80) read (nf2out,*) start, sstop, step kount = int((sstop - start + (step/4.0) )/step) + 1 @@ -78,21 +78,21 @@ c files, and get the synthesis range parameters from the 'dump' file c*****branch to the desired smoothing function write (smitle,1010) smtype ism = 11 - if (smtype .eq. 'l') then + if (smtype == 'l') then lorenflag = 't' - elseif (smtype .eq. 'g') then + elseif (smtype == 'g') then gaussflag = 't' - elseif (smtype .eq. 'v') then + elseif (smtype == 'v') then rotateflag = 't' - elseif (smtype .eq. 'c') then + elseif (smtype == 'c') then rotateflag = 't' gaussflag = 't' - elseif (smtype .eq. 'm') then + elseif (smtype == 'm') then macroflag = 't' - elseif (smtype .eq. 'd') then + elseif (smtype == 'd') then macroflag = 't' gaussflag = 't' - elseif (smtype .eq. 'r') then + elseif (smtype == 'r') then macroflag = 't' rotateflag = 't' gaussflag = 't' @@ -101,27 +101,27 @@ c*****branch to the desired smoothing function c*****compute a stellar rotational broadening function; this follows c D. F. Gray, 1976, "The Obs. & Anal. of Stell. Phot", p394-9 - if (rotateflag .eq. 't') then + if (rotateflag == 't') then 32 array = 'GIVE THE STELLAR vsini [0.0]: ' nchars = 30 - if (line .gt. 0) then + if (line > 0) then call getnum (nchars,line+2,vsini,shortnum) - if (vsini .eq. -9999.) vsini = 0. + if (vsini == -9999.) vsini = 0. endif - if (vsini .lt. 0.0) go to 32 + if (vsini < 0.0) go to 32 write (smitle(ism+1:ism+13),1011) vsini ism = ism + 13 31 array = 'GIVE THE LIMB DARKENING COEFFICIENT [0.0]: ' nchars = 43 - if (line .gt. 0) then + if (line > 0) then call getnum (nchars,line+2,limbdark,shortnum) - if (limbdark .eq. -9999.) limbdark = 0. + if (limbdark == -9999.) limbdark = 0. endif - if (limbdark .lt. 0.0) go to 31 + if (limbdark < 0.0) go to 31 write (smitle(ism+1:ism+13),1012) limbdark ism = ism + 13 dlamlim = (start+sstop)/2.*vsini/3.0e5 - if (step .ge. dlamlim) then + if (step >= dlamlim) then rotateflag = 'f' else pi = 3.141527 @@ -131,11 +131,11 @@ c D. F. Gray, 1976, "The Obs. & Anal. of Stell. Phot", p394-9 prot0 = c1 + c2 powerrot = prot0 jdelrot = idint(dlamlim/step) - if (jdelrot .gt. 1000) then + if (jdelrot > 1000) then write (*,1026) smtype = 'e' return - elseif (jdelrot .gt. kount/4) then + elseif (jdelrot > kount/4) then write (*,1028) smtype = 'e' return @@ -151,17 +151,17 @@ c D. F. Gray, 1976, "The Obs. & Anal. of Stell. Phot", p394-9 c*****compute a macroturbulent smoothing function (uses subroutine vmacro) - if (macroflag .eq. 't') then + if (macroflag == 't') then 51 array = 'GIVE THE MACROTURBULENT VELOCITY [0.0]: ' nchars = 39 - if (line .gt. 0) then + if (line > 0) then call getnum (nchars,line+2,vmac,shortnum) - if (vmac .eq. -9999.) vmac = 0. + if (vmac == -9999.) vmac = 0. endif - if (vmac .lt. 0.0) go to 51 + if (vmac < 0.0) go to 51 write (smitle(ism+1:ism+13),1013) vmac ism = ism + 13 - if (vmac .eq. 0.) then + if (vmac == 0.) then macroflag = 'f' else wavemac = (start+sstop)/2.*vmac/3.0e5 @@ -170,16 +170,16 @@ c*****compute a macroturbulent smoothing function (uses subroutine vmacro) wavei = step*i/wavemac pmac(i) = vmacro(wavei) powermac = powermac + 2.0 *pmac(i) - if (pmac(i) .lt. 0.002) then + if (pmac(i) < 0.002) then jdelmac = i exit endif enddo - if (jdelmac .gt. 1000) then + if (jdelmac > 1000) then write (*,1025) wavemac smtype = 'e' return - elseif (jdelmac .gt. kount/4) then + elseif (jdelmac > kount/4) then write (*,1022) smtype = 'e' return @@ -189,17 +189,17 @@ c*****compute a macroturbulent smoothing function (uses subroutine vmacro) c*****compute a Gaussian smoothing function - if (gaussflag .eq. 't') then + if (gaussflag == 't') then 11 array = 'GIVE THE FWHM OF THE GAUSSIAN FUNCTION [0.0]: ' nchars = 46 - if (line .gt. 0) then + if (line > 0) then call getnum (nchars,line+2,fwhmgauss,shortnum) - if (fwhmgauss .eq. -9999.) fwhmgauss = 0. + if (fwhmgauss == -9999.) fwhmgauss = 0. endif - if (fwhmgauss .lt. 0.0) go to 11 + if (fwhmgauss < 0.0) go to 11 write (smitle(ism+1:ism+18),1014) fwhmgauss ism = ism + 18 - if (fwhmgauss .eq. 0.) then + if (fwhmgauss == 0.) then gaussflag = 'f' else sigma = fwhmgauss/2. @@ -208,16 +208,16 @@ c*****compute a Gaussian smoothing function do i=1,1000 p(i) = dexp(-aa*(step*i)**2 ) power = power + 2*p(i) - if (p(i) .lt. 0.02) then + if (p(i) < 0.02) then jdel = i exit endif enddo - if (jdel .gt. 1000) then + if (jdel > 1000) then write (*,1029) sigma smtype = 'e' return - elseif (jdel .gt. kount/4) then + elseif (jdel > kount/4) then write (*,1021) smtype = 'e' return @@ -227,17 +227,17 @@ c*****compute a Gaussian smoothing function c*****compute a Lorenzian smoothing function - if (lorenflag .eq. 't') then + if (lorenflag == 't') then 21 array = 'GIVE THE FWHM OF THE LORENTZIAN FUNCTION [0.0]: ' nchars = 48 - if (line .gt. 0) then + if (line > 0) then call getnum (nchars,line+2,fwhmloren,shortnum) - if (fwhmloren .eq. -9999.) fwhmloren = 0. + if (fwhmloren == -9999.) fwhmloren = 0. endif - if (fwhmloren .lt. 0.0) go to 21 + if (fwhmloren < 0.0) go to 21 write (smitle(ism+1:ism+20),1015) fwhmloren ism = ism + 20 - if (fwhmloren .eq. 0.) then + if (fwhmloren == 0.) then lorenflag = 'f' else sigma = fwhmloren/2. @@ -245,16 +245,16 @@ c*****compute a Lorenzian smoothing function do i=1,1000 p(i) = ((sigma**2)/((sigma**2)+((step*i)**2))) power = power + 2.0 *p(i) - if (p(i) .lt. 0.02) then + if (p(i) < 0.02) then jdel = i exit endif enddo - if (jdel .gt. 1000) then + if (jdel > 1000) then write (*,1030) sigma smtype = 'e' return - elseif (jdel .gt. kount/4) then + elseif (jdel > kount/4) then write (*,1031) smtype = 'e' return @@ -277,29 +277,29 @@ c*****here is the reading/grabbing of stuff preceding the depth array: nisos = 0 do i=1,20 read (nf2out,1002,end=2000) array - if (array(1:7).eq.'ALL abu') then + if (array(1:7)=='ALL abu') then cycle - elseif (array(1:7).eq.'Changin') then + elseif (array(1:7)=='Changin') then abitle (naboff+1:naboff+23) = '[M/H] FOR ALL ELEMENTS:' abitle (naboff+24:naboff+29) = array(32:37) cycle - elseif (array(1:7).eq.'element') then + elseif (array(1:7)=='element') then nabunds = nabunds + 1 - if (control .eq. 'binary ') then - if (nabunds .le. 5) then + if (control == 'binary ') then + if (nabunds <= 5) then ioff = naboff + 16*(nabunds-1) abitle(ioff+1:ioff+2) = array(9:10) abitle(ioff+3:ioff+14) = array(26:37) abitle(ioff+15:ioff+16) = ' ' endif else - if (nabunds .le. 8) then + if (nabunds <= 8) then ioff = naboff + 9*(nabunds-1) abitle(ioff+1:ioff+2) = array(9:10) read (array(26:32),*) abnum - if (abnum .gt. 0) then + if (abnum > 0) then write (abchars,1061) abnum - elseif (abnum .le. -10.) then + elseif (abnum <= -10.) then write (abchars,1062) abnum else write (abchars,1063) abnum @@ -309,20 +309,20 @@ c*****here is the reading/grabbing of stuff preceding the depth array: endif endif cycle - elseif (array(1:7).eq.'Isotopi') then + elseif (array(1:7)=='Isotopi') then nisos = nisos + 1 - if (nisos .le. 6) then + if (nisos <= 6) then read (array(37:46),1050) ratio - if (ratio .ge. 1000.) then + if (ratio >= 1000.) then write (isochars,1054) int(ratio) - elseif (ratio .ge. 100.) then + elseif (ratio >= 100.) then write (isochars,1051) ratio - elseif (ratio .ge. 10.) then + elseif (ratio >= 10.) then write (isochars,1052) ratio else write (isochars,1053) ratio endif - if (nsyn .eq. 1) then + if (nsyn == 1) then ioff = 40*(nisos-1) + 5*(nsyn-1) isoitle(ioff+1:ioff+10) = array(23:32) isoitle(ioff+11:ioff+12) = ': ' @@ -334,7 +334,7 @@ c*****here is the reading/grabbing of stuff preceding the depth array: isoitle(ioff+5:ioff+5) = '/' endif endif - elseif (array(1:7).eq.'MODEL: ') then + elseif (array(1:7)=='MODEL: ') then read (nf2out,1002) array exit endif @@ -349,7 +349,7 @@ c*****here is the actual reading of the depth array c*****here a veiling addition can be added in - if (addflux .gt. 0.0) then + if (addflux > 0.0) then do i=1,kount y(i) = (y(i) + addflux)/(1.0+addflux) enddo @@ -357,7 +357,7 @@ c*****here a veiling addition can be added in c*****apply the rotational broadening if desired - if (rotateflag .eq. 't') then + if (rotateflag == 't') then min = jdelrot + 1 max = kount - jdelrot do i=1,jdelrot @@ -378,7 +378,7 @@ c*****apply the rotational broadening if desired c*****apply the macroturbulent broadening if desired - if (macroflag .eq. 't') then + if (macroflag == 't') then min = jdelmac + 1 max = kount - jdelmac do i=1,jdelmac @@ -400,7 +400,7 @@ c*****apply the macroturbulent broadening if desired c*****apply the Gaussian or Lorenzian smoothing if desired (this c is an either/or situation; only one of these can apply. - if (gaussflag .eq. 't' .or. lorenflag .eq. 't') then + if (gaussflag == 't' .or. lorenflag == 't') then min = jdel + 1 max = kount - jdel do i=1,jdel @@ -429,7 +429,7 @@ c been done to the y-array) into the appropriate array c*****compute the wavelength array; must be done for each synthetic c spectrum because of the way the equivalences were set up - if (iunits .eq. 1) then + if (iunits == 1) then do i=1,kount xsyn(i) = 1.d-4*(start + (i-1)*step) enddo @@ -443,7 +443,7 @@ c spectrum because of the way the equivalences were set up c*****dump the smoothed spectrum in a MONGO-style set of c (wavelength,flux) point pairs write (nf3out,1005) kount,start,sstop,step - if (xsyn(1) .le. 100.0) then + if (xsyn(1) <= 100.0) then write (nf3out,1009) (xsyn(i),chunk(i,nsyn),i=1,kount) else write (nf3out,1008) (xsyn(i),chunk(i,nsyn),i=1,kount) @@ -18,7 +18,7 @@ c****************************************************************************** real*8 inx1,inx2,iny1,iny2 c*****for grid syntheses, dump out relevant information to a file - if (choice .eq. 'g') then + if (choice == 'g') then write (nf6out,3001) syncount write (nf6out,3002) obsitle, moditle, linitle, smitle endif @@ -44,24 +44,24 @@ c*****write smoothing information at the top of the plot call sm_label (smitle) - if (isoitle(1:10) .eq. ' ') then + if (isoitle(1:10) == ' ') then isoitle(1:16) = 'no isotopic data' endif call sm_relocate (-0.120,1.075) call sm_label (isoitle(1:120)) - if (numiso .gt. 3) then + if (numiso > 3) then call sm_relocate (-0.120,1.045) call sm_label (isoitle(121:240)) endif - if (control .eq. 'gridplo' .or. - . control .eq. 'gridsyn' .or. - . control .eq. 'gridend') then + if (control == 'gridplo' .or. + . control == 'gridsyn' .or. + . control == 'gridend') then write (nf6out,3002) isoitle(1:120) write (nf6out,3002) isoitle(121:240) endif c*****define the real plot limits - if (xlo .lt. xhi) then + if (xlo < xhi) then call sm_limits (xlo,xhi,ylo,yhi) iflip = 0 else @@ -75,7 +75,7 @@ c*****define the real plot limits c*****draw and label the box for the spectra call defcolor (1) - if (whichwin .eq. '1of1') then + if (whichwin == '1of1') then idev = 1 call sm_window (1,1,1,1,1,1) else @@ -89,7 +89,7 @@ c*****draw and label the box for the spectra call sm_lweight (2.0) call sm_expand (0.8) call sm_box (1,2,4,4) - if (iflip .eq. 1) then + if (iflip == 1) then array = 'Wavenumber' else array = 'Wavelength' @@ -103,15 +103,15 @@ c*****plot the synthetic spectra call sm_lweight (2.2) call sm_expand (0.7) do i=1,100 - if (pec(i) .ne. 0) go to 111 + if (pec(i) /= 0) go to 111 enddo 111 do j=1,nsyn - if (choice.eq.'h' .or. choice.eq.'f' .or. - . choice.eq.'g') then + if (choice=='h' .or. choice=='f' .or. + . choice=='g') then call defcolor (8) call sm_ltype (j-1) else - if (smterm(1:3) .eq. 'x11') then + if (smterm(1:3) == 'x11') then call defcolor (j+1) call sm_ltype (0) else @@ -120,7 +120,7 @@ c*****plot the synthetic spectra endif endif call sm_connect (xsyn,chunk(1,j),kount) - if (iflip .eq. 1) then + if (iflip == 1) then call sm_relocate (xhi+0.045*(xlo-xhi), . ylo+(0.12+0.06*j)*(yhi-ylo)) call sm_draw (xhi+0.005*(xlo-xhi), @@ -138,10 +138,10 @@ c*****plot the synthetic spectra noff = 80*(j-1) call sm_lweight (2.2) call sm_label (abitle(noff+1:noff+80)) - if ((control .eq. 'gridplo' .or. - . control .eq. 'gridsyn' .or. - . control .eq. 'gridend') .and. - . whichwin.eq.'1of1') then + if ((control == 'gridplo' .or. + . control == 'gridsyn' .or. + . control == 'gridend') .and. + . whichwin=='1of1') then write (nf6out,3002) abitle(noff+1:noff+80) endif enddo @@ -150,10 +150,10 @@ c*****plot the synthetic spectra c*****plot the observed spectrum - if (plotopt .eq. 2) then + if (plotopt == 2) then call defcolor (1) - if (choice.eq.'h' .or. choice.eq.'f' .or. - . choice.eq.'g') then + if (choice=='h' .or. choice=='f' .or. + . choice=='g') then call sm_lweight (4.0) else call sm_lweight (2.2) @@ -163,10 +163,10 @@ c*****plot the observed spectrum style(1) = 43.5 call sm_ptype (style,1) mount = lim2obs - lim1obs + 1 - if (mount .lt. 500) then + if (mount < 500) then call sm_points (xobs(lim1obs),yobs(lim1obs),mount) else - if (histoyes .eq. 1) then + if (histoyes == 1) then call sm_histogram (xobs(lim1obs),yobs(lim1obs),mount) else call sm_connect (xobs(lim1obs),yobs(lim1obs),mount) @@ -174,54 +174,54 @@ c*****plot the observed spectrum endif call sm_lweight (2.2) call sm_expand (0.7) - if (iflip .eq. 1) then + if (iflip == 1) then call sm_relocate (xhi+0.05*(xlo-xhi),ylo+0.12*(yhi-ylo)) else call sm_relocate (xlo+0.05*(xhi-xlo),ylo+0.12*(yhi-ylo)) endif call sm_label (obsitle) endif - if (iflip .eq. 1) then + if (iflip == 1) then call sm_relocate (xhi+0.05*(xlo-xhi),ylo+0.06*(yhi-ylo)) else call sm_relocate (xlo+0.05*(xhi-xlo),ylo+0.06*(yhi-ylo)) endif call sm_label (moditle) - if (whichwin.eq.'1of1' .or. plotopt.ne.2) then + if (whichwin=='1of1' .or. plotopt/=2) then return endif c*****this section of code is executed only if a deviations plot is desired; c find the starting and stopping points in the arrays for the deviations - if (xsyn(kount) .le. xobs(lim1obs)) go to 1000 - if (xsyn(1) .gt. xobs(lim2obs)) go to 1000 - if (xsyn(1) .gt. xobs(lim1obs)) go to 150 + if (xsyn(kount) <= xobs(lim1obs)) go to 1000 + if (xsyn(1) > xobs(lim2obs)) go to 1000 + if (xsyn(1) > xobs(lim1obs)) go to 150 lim3obs = lim1obs do k=2,kount - if (xsyn(k) .gt. xobs(lim3obs)) then + if (xsyn(k) > xobs(lim3obs)) then lim1syn = k - 1 go to 155 endif enddo 150 lim1syn = 1 do l=lim1obs,lim2obs - if (xsyn(lim1syn) .le. xobs(l)) then + if (xsyn(lim1syn) <= xobs(l)) then lim3obs = l go to 155 endif enddo -155 if (xsyn(kount) .lt. xobs(lim2obs)) go to 160 +155 if (xsyn(kount) < xobs(lim2obs)) go to 160 lim4obs = lim2obs do k=lim1syn,kount - if (xsyn(k) .gt. xobs(lim4obs)) then + if (xsyn(k) > xobs(lim4obs)) then lim2syn = k go to 165 endif enddo 160 lim2syn = kount do l=lim3obs,lim2obs - if (xsyn(lim2syn) .lt. xobs(l)) then + if (xsyn(lim2syn) < xobs(l)) then lim4obs = l - 1 go to 165 endif @@ -235,7 +235,7 @@ c of the synthetic spectra is considered sufficient lpoint = lim1syn devsigma = 0. do i=lim3obs,lim4obs -170 if (xsyn(lpoint+1) .lt. xobs(i)) then +170 if (xsyn(lpoint+1) < xobs(i)) then lpoint = lpoint + 1 go to 170 endif @@ -253,7 +253,7 @@ c of the synthetic spectra is considered sufficient c from first set of deviations, define the plot limits, draw and label box - if (j .eq. 1) then + if (j == 1) then yup = -1000. ydown = +1000. do i=lim3obs,lim4obs @@ -283,12 +283,12 @@ c from first set of deviations, define the plot limits, draw and label box c plot the array of deviations - if (choice.eq.'h' .or. choice.eq.'f' .or. - . choice.eq.'g') then + if (choice=='h' .or. choice=='f' .or. + . choice=='g') then call defcolor (8) call sm_ltype (j-1) else - if (smterm(1:3) .eq. 'x11') then + if (smterm(1:3) == 'x11') then call defcolor (j+1) call sm_ltype (0) else @@ -307,7 +307,7 @@ c plot the array of deviations call sm_draw(xhi-0.215*(xhi-xlo), . ydown+(0.10+0.06*j)*(yup-ydown)) call sm_label (array) - if (choice .eq. 'g') then + if (choice == 'g') then noff = 80*(j-1) write (nf6out,3002) abitle(noff+1:noff+80) write (nf6out,3003) devsigma, velsh @@ -317,7 +317,7 @@ c plot the array of deviations c reset the spectrum plot boundaries before exiting - if(xlo .lt. xhi) then + if(xlo < xhi) then call sm_limits (xlo,xhi,ylo,yhi) iflip = 0 else @@ -13,18 +13,18 @@ c*****compute the average and standard deviation average = 0. kount = 0 do l=lim1obs,lim2obs - if (abundout(l) .ne. 999.99) then + if (abundout(l) /= 999.99) then average = average + abundout(l) kount = kount + 1 endif enddo - if (kount .gt. 0) then + if (kount > 0) then average = average/kount endif deviate = 0. - if (kount .gt. 1) then + if (kount > 1) then do l=lim1obs,lim2obs - if (abundout(l) .ne. 999.99) then + if (abundout(l) /= 999.99) then deviate = deviate + (abundout(l)-average)**2 endif enddo @@ -34,7 +34,7 @@ c*****compute the average and standard deviation c*****correlate the abundances with excitation potential, equivalent width, c and wavelength - if (kount .gt. 2) then + if (kount > 2) then epmin = 999. epmax = -999. rwmin = 999. @@ -54,7 +54,7 @@ c and wavelength za = 0. do l=lim1obs,lim2obs - if (abundout(l) .ne. 999.99) then + if (abundout(l) /= 999.99) then c rw = dlog10(wid1comp(l)/wave1(l)) rw = dlog10(width(l)/wave1(l)) x1 = x1 + e(l,1) @@ -68,12 +68,12 @@ c rw = dlog10(wid1comp(l)/wave1(l)) xy = xy + e(l,1)*abundout(l) yz = yz + rw*abundout(l) za = za + wave1(l)*abundout(l) - if (e(l,1) .lt. epmin) epmin = e(l,1) - if (e(l,1) .gt. epmax) epmax = e(l,1) - if (rw .lt. rwmin) rwmin = rw - if (rw .gt. rwmax) rwmax = rw - if (wave1(l) .lt. wvmin) wvmin = wave1(l) - if (wave1(l) .gt. wvmax) wvmax = wave1(l) + if (e(l,1) < epmin) epmin = e(l,1) + if (e(l,1) > epmax) epmax = e(l,1) + if (rw < rwmin) rwmin = rw + if (rw > rwmax) rwmax = rw + if (wave1(l) < wvmin) wvmin = wave1(l) + if (wave1(l) > wvmax) wvmax = wave1(l) endif enddo @@ -13,7 +13,7 @@ c****************************************************************************** do i=1,5 i1 = im/itest(i) - if (i1 .eq. 0) cycle + if (i1 == 0) cycle i2 = im - i1*itest(i) exit enddo @@ -40,11 +40,11 @@ c*****FIRST PASS: For each model, compute a raw synthetic spectrum; c*****the starting do/if loops are for isotopic things only xhyd = 10.0**xsolar(1) do mmod=1,modtot - if (numiso .gt. 0) then - if (nisos .gt. 0) then + if (numiso > 0) then + if (nisos > 0) then do k=1,numiso do l=1,nisos - if (isotope(k) .eq. isospecial(l)) then + if (isotope(k) == isospecial(l)) then do m=1,numisosyn isoabund(k,m) = fracspecial(mmod,l) enddo @@ -59,7 +59,7 @@ c*****read in the model atmospheres and their summary output files line = synpre num = 80 call getcount (num,line) - if (mmod .lt. 10) then + if (mmod < 10) then write (line(num+1:num+1),1013) mmod else write (line(num+1:num+2),1014) mmod @@ -75,7 +75,7 @@ c*****read in the model atmospheres and their summary output files line = modpre num = 80 call getcount (num,line) - if (mmod .lt. 10) then + if (mmod < 10) then write (line(num+1:num+1),1013) mmod else write (line(num+1:num+2),1014) mmod @@ -106,7 +106,7 @@ c*****open the line list file and the strong line list file nchars = 13 call infile ('input ',nflines,'formatted ',0,nchars, . flines,lscreen) - if (dostrong .gt. 0) then + if (dostrong > 0) then nfslines = 32 lscreen = 18 array = 'THE STRONG LINE LIST' @@ -118,7 +118,7 @@ c*****open the line list file and the strong line list file c*****do the syntheses ncall = 1 - if (numpecatom .eq. 0 .or. numatomsyn .eq. 0) then + if (numpecatom == 0 .or. numatomsyn == 0) then isorun = 1 nlines = 0 mode = 3 @@ -188,11 +188,11 @@ c*****read back the syntheses, compute the weighted average lincount = 0 50 read (newunit,1001) line lincount = lincount + 1 - if (mmod .eq. 1) write (holdline(j,lincount),1001) line - if (line(1:5) .eq. 'MODEL') then + if (mmod == 1) write (holdline(j,lincount),1001) line + if (line(1:5) == 'MODEL') then read (newunit,1001) line lincount = lincount + 1 - if (mmod .eq. 1) write (holdline(j,lincount),1001) line + if (mmod == 1) write (holdline(j,lincount),1001) line read (line,*) wavemod1, wavemod2, wavestep nw1 = nint(1000.*wavemod1) nw2 = nint(1000.*wavemod2) @@ -223,12 +223,12 @@ c write the average spectrum back to disk. c*****now plot the spectrum - if (plotopt.eq.2 .and. specfileopt.gt.0) then + if (plotopt==2 .and. specfileopt>0) then nfobs = 33 lscreen = 12 array = 'THE OBSERVED SPECTRUM' nchars = 21 - if (specfileopt.eq.1 .or. specfileopt.eq.3) then + if (specfileopt==1 .or. specfileopt==3) then call infile ('input ',nfobs,'unformatted',2880,nchars, . fobs,lscreen) else @@ -236,7 +236,7 @@ c*****now plot the spectrum . fobs,lscreen) endif endif - if (plotopt .ne. 0) then + if (plotopt /= 0) then nf2out = nf9out f2out = f9out nf2out = 21 @@ -251,7 +251,7 @@ c*****now plot the spectrum nchars = 25 call infile ('output ',nf3out,'formatted ',0,nchars, . f3out,lscreen) - if (f5out .ne. 'optional_output_file') then + if (f5out /= 'optional_output_file') then nf5out = 26 lscreen = lscreen + 2 array = 'POSTSCRIPT PLOT OUTPUT' @@ -259,7 +259,7 @@ c*****now plot the spectrum call infile ('output ',nf5out,'formatted ',0,nchars, . f5out,lscreen) endif - if (plotopt .eq. 3) then + if (plotopt == 3) then call smooth (-1,ncall) choice = 'q' else @@ -268,9 +268,9 @@ c*****now plot the spectrum c*****if needed, loop back with abundance changes - if (choice .eq. 'n') then + if (choice == 'n') then call chabund - if (choice .eq. 'q') go to 60 + if (choice == 'q') go to 60 endif endif @@ -15,19 +15,19 @@ c****************************************************************************** c*****initialize the synthesis
write (nf1out,1101)
write (nf2out,1002) moditle(1:73)
- if (iunits .eq. 1) then
+ if (iunits == 1) then
write (nf2out,1103) start/1.d4, sstop/1.d4,
. step/1.d4, delta/1.d4
else
write (nf2out,1102) start,sstop,step,delta
endif
- if (iraf .eq. 1) then
+ if (iraf == 1) then
npoint = (sstop-start)/step
write (nf4out,1104) npoint,wave,wave,step,step
write (nf4out,1105)
write (nf4out,1106) moditle
do j=1,93
- if (pec(j) .gt. 0 ) then
+ if (pec(j) > 0 ) then
dummy1(j) = dlog10(xabund(j)) + 12.0
write (nf4out,1107) names(j),dummy1(j)
endif
@@ -35,7 +35,7 @@ c*****initialize the synthesis write (nf4out,1108) vturb(1)
write (nf4out,1109)
endif
- if (mode .ne. 4) then
+ if (mode /= 4) then
lim1line = 0
lim2line = 0
lim1obs = 0
@@ -56,15 +56,15 @@ c spectrum wavelength, if needed do n=1,kount
num = num + 1
wave = oldstart + (n-1)*step
- if (dabs(wave-wavl)/wave .ge. 0.001) then
+ if (dabs(wave-wavl)/wave >= 0.001) then
wavl = wave
call opacit (2,wave)
- if (modprintopt .ge. 2)
+ if (modprintopt >= 2)
. write (nf1out,1001) wave,(kaplam(i),i=1,ntau)
call cdcalc (1)
first = 0.4343*cd(1)
flux = rinteg(xref,cd,dummy1,ntau,first)
- if (iunits .eq. 1) then
+ if (iunits == 1) then
write (nf1out,1003) 1.d-4*wave,flux
else
write (nf1out,1004) wave,flux
@@ -75,9 +75,9 @@ c spectrum wavelength, if needed c*****find the appropriate set of lines for this wavelength, reading
c in a new set if this is the initial depth calculation or if
c needed because the line list end has been reached
- if (mode .eq. 3) then
+ if (mode == 3) then
20 call linlimit
- if (lim2line .lt. 0) then
+ if (lim2line < 0) then
call inlines (2)
call nearly (1)
go to 20
@@ -90,7 +90,7 @@ c needed because the line list end has been reached c*****compute a spectrum depth at this point; if there are no absorption
c lines in the interval then just set the depth to zero without
c extensive line calculations
- if (lineflag .lt. 0) then
+ if (lineflag < 0) then
d(num) = 0.
else
call taukap
@@ -98,45 +98,45 @@ c extensive line calculations first = 0.4343*cd(1)
d(num) = rinteg(xref,cd,dummy1,ntau,first)
endif
- if (mod(n,10) .eq. 0) then
- if (iraf .eq. 1) then
+ if (mod(n,10) == 0) then
+ if (iraf == 1) then
do j = 1,10
dd(num-10+j) = 1. - d(num-10+j)
enddo
write (nf4out,1110) (dd(num-10+j),j=1,10)
endif
- if (iunits .eq. 1) then
+ if (iunits == 1) then
wave3 = 1.d-4*(wave - 9.0*step)
write (nf1out,1112) wave3,(d(num-10+j),j=1,10)
else
wave3 = wave - 9.0*step
write (nf1out,1111) wave3,(d(num-10+j),j=1,10)
endif
- if (nf2out .gt. 0) write (nf2out,1110) (d(num-10+j),j=1,10)
+ if (nf2out > 0) write (nf2out,1110) (d(num-10+j),j=1,10)
endif
- if (num .ge. 5000) num = 0
+ if (num >= 5000) num = 0
enddo
c*****finish the synthesis
nn = mod(num,10)
- if (nn .ne. 0) then
- if (iraf .eq. 1) then
+ if (nn /= 0) then
+ if (iraf == 1) then
do j=1,nn
dd(num-nn+j) = 1. - d(num-nn+j)
enddo
write (nf4out,1110) (dd(num-nn+j),j=1,nn)
endif
- if (iunits .eq. 1) then
+ if (iunits == 1) then
wave3 = 1.d-4*(wave - 9.0*step)
write (nf1out,1112) wave3,(d(num-nn+j),j=1,nn)
else
wave3 = wave - 9.0*step
write (nf1out,1111) wave3,(d(num-nn+j),j=1,nn)
endif
- if (nf2out .gt. 0) write (nf2out,1110) (d(num-nn+j),j=1,nn)
+ if (nf2out > 0) write (nf2out,1110) (d(num-nn+j),j=1,nn)
endif
- if (iunits .eq. 1) then
+ if (iunits == 1) then
write (nf1out,1113) 1.d-4*wave
else
write (nf1out,1114) wave
@@ -32,14 +32,14 @@ c spectra, and (if desired) IRAF-style smoothed spectra nchars = 20 call infile ('output ',nf2out,'formatted ',0,nchars, . f2out,lscreen) - if (plotopt .ne. 0) then + if (plotopt /= 0) then nf3out = 22 lscreen = lscreen + 2 array = 'SMOOTHED SYNTHESES OUTPUT' nchars = 25 call infile ('output ',nf3out,'formatted ',0,nchars, . f3out,lscreen) - if (f5out .ne. 'optional_output_file') then + if (f5out /= 'optional_output_file') then nf5out = 26 lscreen = lscreen + 2 array = 'POSTSCRIPT PLOT OUTPUT' @@ -48,7 +48,7 @@ c spectra, and (if desired) IRAF-style smoothed spectra . f5out,lscreen) endif endif - if (iraf .ne. 0) then + if (iraf /= 0) then nf4out = 23 lscreen = lscreen + 2 array = 'IRAF ("rtext") OUTPUT' @@ -75,7 +75,7 @@ c*****open the line list file and the strong line list file nchars = 13 call infile ('input ',nflines,'formatted ',0,nchars, . flines,lscreen) - if (dostrong .gt. 0) then + if (dostrong > 0) then nfslines = 32 lscreen = lscreen + 2 array = 'THE STRONG LINE LIST' @@ -88,7 +88,7 @@ c*****open the line list file and the strong line list file c*****do the syntheses choice = '1' do i=1,100 - if (i .eq. 100) then + if (i == 100) then write (*,1002) stop endif @@ -100,21 +100,21 @@ c*****now either don't make a plot (plotopt = 0) c plot the synthetic spectrum, (plotopt = 1) c plot syntheses and observation (plotopt = 2) c or just smooth the syntheses (plotopt = 3) - if (choice .eq. 'n') then + if (choice == 'n') then ncall = 2 else ncall = 1 endif - if (plotopt .eq. 0) then + if (plotopt == 0) then choice = 'q' - elseif (plotopt .eq. 1) then + elseif (plotopt == 1) then call pltspec (lscreen,ncall) - elseif (plotopt .eq. 2) then + elseif (plotopt == 2) then nfobs = 33 lscreen = lscreen + 2 array = 'THE OBSERVED SPECTRUM' nchars = 21 - if (specfileopt .eq. 1) then + if (specfileopt == 1) then call infile ('input ',nfobs,'unformatted',2880,nchars, . fobs,lscreen) else @@ -122,14 +122,14 @@ c or just smooth the syntheses (plotopt = 3) . fobs,lscreen) endif call pltspec (lscreen,ncall) - elseif (plotopt .eq. 3) then + elseif (plotopt == 3) then call smooth (-1,ncall) choice = 'q' else write (*,1001) stop endif - if (choice .eq. 'q') then + if (choice == 'q') then call finish (0) exit endif @@ -30,27 +30,27 @@ c*****open the model table input file and the summary table output file c*****read the table input for integrated light EW matching - if (option .eq. 1) then + if (option == 1) then read (nftable,1001) line - if (line(1:5) .ne. 'abpop' ) then + if (line(1:5) /= 'abpop' ) then write(*,*) 'OOPS! WRONG TABLE FOR ABPOP; I QUIT!' stop endif do i=1,1000 read (nftable,1001) line - if (line(1:5) .eq. 'modpr') then + if (line(1:5) == 'modpr') then call blankstring (modpre) modpre(1:70) = line(11:80) - elseif (line(1:5) .eq. 'synpr') then + elseif (line(1:5) == 'synpr') then call blankstring (synpre) synpre(1:70) = line(11:80) - elseif (line(1:5) .eq. 'title') then + elseif (line(1:5) == 'title') then call blankstring (abitle) popitle(1:74) = line(7:80) write (nf7out,1002) popitle(1:73) - elseif (line(1:5) .eq. 'model') then + elseif (line(1:5) == 'model') then do mmod=1,100 - if (mmod .eq. 100) then + if (mmod == 100) then write(*,*) 'MORE THAN 99 MODELS; I QUIT!' stop endif @@ -66,60 +66,60 @@ c*****read the table input for integrated light spectrum syntheses nabs = 0 nisos = 0 read (nftable,1001) line(1:80) - if (line(1:6) .ne. 'synpop') then + if (line(1:6) /= 'synpop') then write (*,*) 'OOPS! WRONG TABLE FOR SYNPOP; I QUIT!' stop endif do k=1,1000 call blankstring (line) read (nftable,1001) line(1:80) - if (line(1:5) .eq. 'modpr') then + if (line(1:5) == 'modpr') then call blankstring (modpre) modpre(1:70) = line(11:80) - elseif (line(1:5) .eq. 'synpr') then + elseif (line(1:5) == 'synpr') then call blankstring (synpre) synpre(1:70) = line(11:80) - elseif (line(1:5) .eq. 'title') then + elseif (line(1:5) == 'title') then call blankstring (popitle) popitle(1:74) = line(7:80) write (nf7out,1003) popitle(1:74) - elseif (line(1:5) .eq. 'abund') then + elseif (line(1:5) == 'abund') then read (line(12:80),*) nabs - if (nabs .gt. 0) then + if (nabs > 0) then read (nftable,1001) line(1:80) read (line(1:80),*) (elspecial(i),i=1,nabs) write (nf7out,1009) (nint(elspecial(i)),i=1,nabs) else write (nf7out,1010) endif - elseif (line(1:5) .eq. 'isoto') then + elseif (line(1:5) == 'isoto') then read (line(10:80),*) nisos - if (nisos .gt. 0) then + if (nisos > 0) then read (nftable,1001) line(1:80) read (line(1:80),*) (isospecial(i),i=1,nisos) write (nf7out,1011) (isospecial(i),i=1,nisos) else write (nf7out,1012) endif - elseif (line(1:5) .eq. 'model') then + elseif (line(1:5) == 'model') then do mmod=1,100 - if (mmod .eq. 100) then + if (mmod == 100) then write (*,*) 'MORE THAN 99 MODELS; I QUIT!' stop endif - if (nabs.le.0 .and. nisos.le.0) then + if (nabs<=0 .and. nisos<=0) then read (nftable,*,end=10) j, radius(mmod), . relcount(mmod) write (nf7out,1006) j, radius(mmod), . relcount(mmod) - elseif (nabs.gt.0 .and. nisos.le.0) then + elseif (nabs>0 .and. nisos<=0) then read (nftable,*,end=10) j, radius(mmod), . relcount(mmod), . (abspecial(mmod,i),i=1,nabs) write (nf7out,1006) j, radius(mmod), . relcount(mmod), . (abspecial(mmod,i),i=1,nabs) - elseif (nabs.le.0 .and. nisos.gt.0) then + elseif (nabs<=0 .and. nisos>0) then read (nftable,*,end=10) j, radius(mmod), . relcount(mmod), . (fracspecial(mmod,i),i=1,nisos) @@ -23,7 +23,7 @@ c*****compute the total line opacity at each depth dummy1(i) = tauref(i)*kapnu(i)/(0.4343*kapref(i)) c*****do the same for the strong lines - if (dostrong .gt. 0) then + if (dostrong > 0) then do j=nlines+1,nlines+nstrong v = 2.997929d10*dabs(wave-wave1(j))/ . (wave1(j)*dopp(j,i)) @@ -15,7 +15,7 @@ c****************************************************************************** c*****compute the wavelength array ntotal = (sstop - start)/step + 1.3 - if (ntotal .gt. 5000) then + if (ntotal > 5000) then write (nf1out,1002) ntotal write (nf2out,1002) ntotal return @@ -34,7 +34,7 @@ c*****use the RINTEG routine to do an integration c*****Then recheck using Simpson's rule ntot = ntotal - if(ntotal/2*2 - ntotal .eq. 0) ntot = ntotal - 1 + if(ntotal/2*2 - ntotal == 0) ntot = ntotal - 1 answer = d(1) + 4.*d(2) + d(ntot) ntot = ntot - 2 do i=3,ntot,2 @@ -22,7 +22,7 @@ c*****begin with some calculations leading to a c6 value ("unsold") c*****Ca II "IR triplet" lines at the Ca II K line at 3934 A - if (iatom.eq.201 .and. iwave.eq.3933) then + if (iatom==201 .and. iwave==3933) then do i=1,ntau gammaa = 1.45d+8 gnature = gammaa + 0.5*1.5d-9*t(i)**(1/3)*numdens(1,1,i) @@ -37,10 +37,10 @@ c*****Ca II "IR triplet" lines at the Ca II K line at 3934 A c*****Ca II "IR triplet" lines at 8498, 8542, and 8662 A - elseif (iatom10.eq.201 .and. - . (iwave.eq.8498.or. - . iwave.eq.8542.or. - . iwave.eq.8662)) then + elseif (iatom10==201 .and. + . (iwave==8498.or. + . iwave==8542.or. + . iwave==8662)) then write (nf1out,1000) do i=1,ntau nhe = xabund(2)*nhtot(i) @@ -60,7 +60,7 @@ c*****Ca II "IR triplet" lines at 8498, 8542, and 8662 A c*****Ca I 6717 A - elseif (iatom10.eq.200 .and. iwave.eq.6717) then + elseif (iatom10==200 .and. iwave==6717) then write (nf1out,1002) iwave do i=1,ntau nhe = xabund(2)*nhtot(i) @@ -80,13 +80,13 @@ c*****Ca I 6717 A c*****Ca I 6318, 6343, 6361 A autoionization lines - elseif (iatom10.eq.200 .and. - . (iwave.eq.6318 .or. - . iwave.eq.6343 .or. - . iwave.eq.6361)) then + elseif (iatom10==200 .and. + . (iwave==6318 .or. + . iwave==6343 .or. + . iwave==6361)) then write (nf1out,1005) iwave do i=1,ntau - if (dampnum(j) .eq. 0) then + if (dampnum(j) == 0) then gnature = 1.5d12 else gnature = dampnum(j)*1.5d12 @@ -99,7 +99,7 @@ c*****Ca I 6318, 6343, 6361 A autoionization lines c*****Na I lines - elseif (iatom .eq. 110) then + elseif (iatom == 110) then write (nf1out,1003) iwave do i=1,ntau gnature = 2.21e+15/wave1(j)**2 @@ -116,11 +116,11 @@ c*****Na I lines c*****CH autoionization line at 3693 A - elseif (iatom10.eq.1060 .and. - . iwave.eq.3693) then + elseif (iatom10==1060 .and. + . iwave==3693) then write (nf1out,1006) iwave do i=1,ntau - if (dampnum(j) .eq. 0) then + if (dampnum(j) == 0) then gnature = 4.0d11 else gnature = dampnum(j)*4.0d11 @@ -11,7 +11,7 @@ c****************************************************************************** dimension scale(4) data scale/0.001,0.01,0.1,1.0/ - if (level .gt. 500) then + if (level > 500) then temp = level else temp = t(level) @@ -21,9 +21,9 @@ c****************************************************************************** ion = nint(10.*(atom - float(iatom))) j = 4*(iatom-1) + ion + 1 - if (ion .eq. 0) then + if (ion == 0) then chix = xchi1(iatom) - elseif (ion .eq. 1) then + elseif (ion == 1) then chix = xchi2(iatom) else chix = xchi3(iatom) @@ -38,7 +38,7 @@ c****************************************************************************** k2 = nudata(i,j) - k1*100000 k3 = k2/10 kscale = k2 - k3*10 - if (mod(it,2) .eq. 0) then + if (mod(it,2) == 0) then p1 = float(k3)*scale(kscale) k1 = nudata(i+1,j)/100000 kscale = mod(nudata(i+1,j),10) @@ -46,10 +46,10 @@ c****************************************************************************** else p1 = float(k1)*scale(kscale) p2 = float(k3)*scale(kscale) - if (dt .ge. 0.) go to 13 - if (kscale .gt. 1.) go to 13 + if (dt >= 0.) go to 13 + if (kscale > 1.) go to 13 kp1 = p1 - if (kp1 .ne. idint(p2+0.5)) go to 13 + if (kp1 /= idint(p2+0.5)) go to 13 pmin = kp1 endif 13 ucalc = dmax1(pmin,p1+(p2-p1)*dt) @@ -29,16 +29,16 @@ c files, and get the synthesis range parameters from the 'dump' file rewind nf2out rewind nf3out 55 read (nf2out,1002) moditle - if (moditle(1:15).eq.' element' .or. - . moditle(1:15).eq.' ALL abundances' .or. - . moditle(1:15).eq.'Isotope Ratio: ') go to 55 + if (moditle(1:15)==' element' .or. + . moditle(1:15)==' ALL abundances' .or. + . moditle(1:15)=='Isotope Ratio: ') go to 55 read (nf2out,*) start, sstop, step kount = nint((sstop - start + (step/4.0) )/step) + 1 rewind nf2out c*****the first time through, read in the Gaussian FWHM array - if (istart .eq. 0) then + if (istart == 0) then istart = 1 nfsmooth = 35 array = 'SMOOTHING FWHM DATA' @@ -48,7 +48,7 @@ c*****the first time through, read in the Gaussian FWHM array j = 1 39 read (nfsmooth,*,end=40) wavefwhm(j), fwhm(j) j = j + 1 - if (j .le. 100) then + if (j <= 100) then go to 39 else istat = ivcleof (line,1) @@ -66,26 +66,26 @@ c*****now read in the raw spectrum and flip to a depth scale abitle(noff+1:noff+12) = '[M/H] 0.00 ' nabunds = 0 41 read (nf2out,1002,end=2000) array - if (array(1:15).eq.' ALL abundances') then + if (array(1:15)==' ALL abundances') then abitle(noff+6:noff+10) = array(56:60) go to 41 - elseif (array(1:15).eq.' element') then + elseif (array(1:15)==' element') then nabunds = nabunds + 1 - if (nabunds .le. 7) then + if (nabunds <= 7) then ioff = noff + 12 + 9*(nabunds-1) abitle(ioff+1:ioff+2) = array(17:18) abitle(ioff+3:ioff+7) = array(34:38) abitle(ioff+8:ioff+9) = ' ' endif go to 41 - elseif (array(1:15).eq.'Isotope Ratio: ') then + elseif (array(1:15)=='Isotope Ratio: ') then nabunds = nabunds + 1 - if (nabunds .le. 7) then + if (nabunds <= 7) then ioff = noff + 12 + 9*(nabunds-1) abitle(ioff+1:ioff+4) = array(27:30) abitle(ioff+5:ioff+5) = ' ' do k=33,44 - if (array(k:k) .ne. ' ') then + if (array(k:k) /= ' ') then abitle(ioff+6:ioff+8) = array(k:k+2) go to 60 endif @@ -106,22 +106,22 @@ c Gaussian smoothing will need to be different at each step i = 0 oldhalf = 0. 25 i = i + 1 - if (i .gt. kount) go to 90 + if (i > kount) go to 90 synstep = start + (i-1)*step c*****interpolate linearly in the FWHM array to get the appropriate value c for the current wavelength step - if (synstep .le. wavefwhm(1)) then + if (synstep <= wavefwhm(1)) then half = fwhm(1) - elseif (synstep .ge. wavefwhm(jtotfwhm)) then + elseif (synstep >= wavefwhm(jtotfwhm)) then half = fwhm(jtotfwhm) else do j=2,jtotfwhm - if (synstep .le. wavefwhm(j)) then + if (synstep <= wavefwhm(j)) then half = fwhm(j-1) + (synstep-wavefwhm(j-1))* . (fwhm(j)-fwhm(j-1))/(wavefwhm(j)-wavefwhm(j-1)) - if (half .gt. 0.) then + if (half > 0.) then go to 10 else go to 15 @@ -132,7 +132,7 @@ c for the current wavelength step c*****compute the Gaussian smoothing function, if needed -10 if (dabs(half-oldhalf)/half .lt. 0.03) go to 50 +10 if (dabs(half-oldhalf)/half < 0.03) go to 50 oldhalf = half sigma = half/2 aa = 0.6932/sigma**2 @@ -140,7 +140,7 @@ c*****compute the Gaussian smoothing function, if needed do k=1,1000 p(k) = dexp(-aa*(step*k)**2 ) power = power + 2*p(k) - if (p(k) .lt. 0.05) then + if (p(k) < 0.05) then jdel = k min = jdel + 1 max = kount - jdel @@ -161,7 +161,7 @@ c*****if no smoothing, just equate the smoothed to the unsmoothed point c*****otherwise smooth the spectrum -50 if (i.lt.min .or. i.gt.max) then +50 if (i<min .or. i>max) then z(i) = y(i) else z(i) = p0*y(i) @@ -189,7 +189,7 @@ c spectrum because of the way the equivalences were set up c*****dump the smoothed spectrum in a MONGO-style set of c (wavelength,flux) point pairs write (nf3out,1005) kount, start, sstop, step - if (xsyn(1) .le. 100.0) then + if (xsyn(1) <= 100.0) then write (nf3out,1009) (xsyn(i),z(i),i=1,kount) else write (nf3out,1008) (xsyn(i),z(i),i=1,kount) @@ -85,7 +85,7 @@ c . 0.000190,0.000180,0.000170,0.000161,0.000152/ - if (xprofile.lt.xrt(1) .or. xprofile.ge.xrt(250)) then + if (xprofile<xrt(1) .or. xprofile>=xrt(250)) then vmacro = 0. return endif @@ -9,9 +9,9 @@ c****************************************************************************** implicit real*8 (a-h,o-z) a2 = a*a v2 = v*v - if(a .eq. 0.) go to 62 - if(a .le. 0.2) go to 50 - if(a.le.1.4 .and. a+v.le.3.2) go to 30 + if(a == 0.) go to 62 + if(a <= 0.2) go to 50 + if(a<=1.4 .and. a+v<=3.2) go to 30 c case 1, a .gt. 1.4 or (a .gt. 0.2 and a+v .gt. 3.2) @@ -22,8 +22,8 @@ c case 1, a .gt. 1.4 or (a .gt. 0.2 and a+v .gt. 3.2) return 30 u=0.979895023-0.962846325*a+0.532770573*a2-0.122727278*a*a2 31 h0 = dexp(-v2) - if(v .ge. 2.4) go to 35 - if(v .ge. 1.3) go to 33 + if(v >= 2.4) go to 35 + if(v >= 1.3) go to 33 c case 2, 0.2 .lt. a .le. 1.4 and a+v .le. 3.2 @@ -36,7 +36,7 @@ c case 2, 0.2 .lt. a .le. 1.4 and a+v .le. 3.2 35 h1=(0.554153432+0.278711796*v-0.188325687*v2+0.042991293*v*v2 1 -0.003278278*v2*v2)/(v2 - 1.5) 36 h2 = (1.0 - 2.0*v2)*h0 - if(a .le. 0.2) go to 52 + if(a <= 0.2) go to 52 h1 = h1 + 1.1283790*h0 h2p = h2 h2 = h2 - h0 + 1.1283790*h1 @@ -48,7 +48,7 @@ c case 2, 0.2 .lt. a .le. 1.4 and a+v .le. 3.2 c case 3, a .le. 0.2 and v .lt. 5.0 -50 if(v .ge. 5.0) go to 60 +50 if(v >= 5.0) go to 60 go to 31 52 voigt = h0 + h1*a + h2*a2 voigt = voigt/1.772454 @@ -8,12 +8,12 @@ c****************************************************************************** real*4 point real*8 c(9) - if (c(9) .eq. 1.) go to 20 - if (c(9) .eq. 2.) go to 40 - if (c(9) .eq. 3.) go to 30 + if (c(9) == 1.) go to 20 + if (c(9) == 2.) go to 40 + if (c(9) == 3.) go to 30 c*****no wavelength solution exists - if (c(1).eq.0. .or. c(2).eq.0.) then + if (c(1)==0. .or. c(2)==0.) then wavecalc = point c*****an ordinary polynomial solution @@ -80,22 +80,22 @@ c*****compute the line opacities c*****calculate continuum quantities at the line list wavelength middle wave = (wave1(1)+wave1(nlines))/2. call opacit (2,wave) - if (modprintopt .ge. 2) + if (modprintopt >= 2) . write(nf1out,1006) wave,(kaplam(i),i=1,ntau) c*****divide the lines into keepers and discards do j=1,nlines residual = 10.*atom1(j) - dble(nint(10.*(atom1(j)))) - if (strength(j)/kaplam(jtau5) .ge. xratio) then - if (atom1(j) .lt. 100.) then - if (residual .gt. 0. .and. dampnum(j) .gt. 0.) then + if (strength(j)/kaplam(jtau5) >= xratio) then + if (atom1(j) < 100.) then + if (residual > 0. .and. dampnum(j) > 0.) then write (nf8out,1007) wave1(j), atom1(j), e(j,1), . dlog10(gf(j)), dlog10(dampnum(j)), dlog10(strength(j)) - else if (residual .gt. 0. .and. dampnum(j) .le. 0.) then + else if (residual > 0. .and. dampnum(j) <= 0.) then write (nf8out,1007) wave1(j), atom1(j), e(j,1), . dlog10(gf(j)), 0.0, dlog10(strength(j)) - else if (residual .le. 0. .and. dampnum(j) .gt. 0.) then + else if (residual <= 0. .and. dampnum(j) > 0.) then write (nf8out,1004) wave1(j), atom1(j), e(j,1), . dlog10(gf(j)), dlog10(dampnum(j)), dlog10(strength(j)) else @@ -103,7 +103,7 @@ c*****divide the lines into keepers and discards . dlog10(gf(j)), 0.0, dlog10(strength(j)) endif else - if (residual .gt. 0.) then + if (residual > 0.) then write (nf8out, 1008) wave1(j), atom1(j), e(j,1), . dlog10(gf(j)), d0(j), dlog10(strength(j)) else @@ -112,14 +112,14 @@ c*****divide the lines into keepers and discards endif endif else - if (atom1(j) .lt. 100.) then - if (residual .gt. 0. .and. dampnum(j) .gt. 0.) then + if (atom1(j) < 100.) then + if (residual > 0. .and. dampnum(j) > 0.) then write (nf9out,1007) wave1(j), atom1(j), e(j,1), . dlog10(gf(j)), dlog10(dampnum(j)), dlog10(strength(j)) - else if (residual .gt. 0. .and. dampnum(j) .le. 0.) then + else if (residual > 0. .and. dampnum(j) <= 0.) then write (nf9out,1007) wave1(j), atom1(j), e(j,1), . dlog10(gf(j)), 0.0, dlog10(strength(j)) - else if (residual .le. 0. .and. dampnum(j) .gt. 0.) then + else if (residual <= 0. .and. dampnum(j) > 0.) then write (nf9out,1004) wave1(j), atom1(j), e(j,1), . dlog10(gf(j)), dlog10(dampnum(j)), dlog10(strength(j)) else @@ -127,7 +127,7 @@ c*****divide the lines into keepers and discards . dlog10(gf(j)), 0.0, dlog10(strength(j)) endif else - if (residual .gt. 0.) then + if (residual > 0.) then write (nf9out, 1008) wave1(j), atom1(j), e(j,1), . dlog10(gf(j)), d0(j), dlog10(strength(j)) else @@ -137,7 +137,7 @@ c*****divide the lines into keepers and discards endif endif enddo - if (nlines +nstrong .eq. 2500) then + if (nlines +nstrong == 2500) then call inlines (6) go to 1 endif diff --git a/Writenumber.f b/Writenumber.f index 66cdb6c..3d01926 100755 --- a/Writenumber.f +++ b/Writenumber.f @@ -11,9 +11,9 @@ c****************************************************************************** lognum = alog10(abs(xnum)) - if (lognum .ge. 6.) then + if (lognum >= 6.) then write (errmess,1002) - elseif (lognum .ge. 0.) then + elseif (lognum >= 0.) then numdec = 5 - nint(lognum) write (errmess,1001) numdec else |