      PROGRAM WIDTH9mf
c      Modifed 27.04.95
ChSt   Letzter Kommentar 26.März.2004

      IMPLICIT REAL*8 (A-H,O-Z)
ChSt  kw ist die maximal erlaubte Schichtanzahl der Modellatmosphären.
ChSt  sollten wir auf 300 setzen!  
ChSt  im folgenden ist zu sehen, daß die meisten physikalischen Größen
ChSt  wie zb. Temperatur (T,TKEV,TK...) oder Druck(PTOTAL,PTURB,...)
ChSt  für jede Schicht separat bestimmt und gespeichert werden müssen
ChSt  (T(kw),...)
      PARAMETER (kw=99)
	common/vmag/vmag
	common /potion/potion(594)
      COMMON /ABROSS/ABROSS(kw),TAUROS(kw)
      COMMON /ABTOT/ABTOT(kw),ALPHA(kw)
      COMMON /CONT/ABTOTC(kw),ALPHAC(kw),TAUNUC(kw),SNUC(kw),HNUC(kw),
     1             JNUC(kw),JMINSC(kw),RESIDC(kw)
      REAL*8 JNUC,JMINSC
      COMMON /CONV/DLTDLP(kw),HEATCP(kw),DLRDLT(kw),VELSND(kw),
     1             GRDADB(kw),HSCALE(kw),FLXCNV(kw),VCONV(kw),MIXLTH,
     2             OVERWT,FLXCNV0(kw),FLXCNV1(kw),IFCONV
      REAL*8 MIXLTH
      COMMON /DEPART/BHYD(kw,6),BMIN(kw),NLTEON
      COMMON /ELEM/ABUND(99),ATMASS(99),ELEM(99)
      COMMON /EDENS/EDENS(kw),IFEDNS
      COMMON /FILENAME/FILENAME
      CHARACTER*60 FILENAME
      COMMON /FLUX/ FLUX,FLXERR(kw),FLXDRV(kw),FLXRAD(kw)
      COMMON /FREQ/FREQ,FREQLG,EHVKT(kw),STIM(kw),BNU(kw)
      COMMON /FRESET/FRESET(1563),RCOSET(1563),NULO,NUHI,NUMNU
      COMMON /HEIGHT/HEIGHT(kw)
      COMMON /IF/IFCORR,IFPRES,IFSURF,IFSCAT,TAUSCAT,IFMOL
      COMMON /IFOP/IFOP(20)
      COMMON /IONS/XNFPH(kw,2),XNFPHE(kw,3),XNFH(kw,2),XNFHE(kw,3)
      COMMON /ITER/ ITER,IFPRNT(15),IFPNCH(15),NUMITS
      COMMON /JUNK/TITLE(74),FREQID(6),WLTE,XSCALE
      COMMON /MUS/ANGLE(20),SURFI(20),NMU
      COMMON /OPS/AHYD(kw),AH2P(kw),AHMIN(kw),SIGH(kw),AHE1(kw),
     1       AHE2(kw),AHEMIN(kw),SIGHE(kw),ACOOL(kw),ALUKE(kw),AHOT(kw),
     2            SIGEL(kw),SIGH2(kw),AHLINE(kw),ALINES(kw),SIGLIN(kw),
     3            AXLINE(kw),SIGXL(kw),AXCONT(kw),SIGX(kw),SHYD(kw),
     4            SHMIN(kw),SHLINE(kw),SXLINE(kw),SXCONT(kw)
      COMMON /OPTOT/ACONT(kw),SCONT(kw),ALINE(kw),SLINE(kw),SIGMAC(kw),
     1              SIGMAL(kw)
      COMMON /PTOTAL/PTOTAL(kw)
      COMMON /PUT/PUT,IPUT
      COMMON /PZERO/PZERO,PCON,PRADK0,PTURB0,KNU(kw),PRADK(kw),RADEN(kw)
      REAL*8 KNU
      COMMON /RAD/ ACCRAD(kw),PRAD(kw)
      COMMON /RHOX/RHOX(kw),NRHOX
      COMMON /STATE/P(kw),XNE(kw),XNATOM(kw),RHO(kw)
      COMMON /TAUSHJ/TAUNU(kw),SNU(kw),HNU(kw),JNU(kw),JMINS(kw)
      REAL*8 JNU,JMINS
      COMMON /STEPLG/STEPLG,TAU1LG,KRHOX
      COMMON /TEFF/TEFF,GRAV,GLOG
      COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),ITEMP
      COMMON /TURBPR/VTURB(kw),PTURB(kw),TRBFDG,TRBCON,TRBPOW,TRBSND,
     1               IFTURB
      COMMON /WAVEY/WBEGIN,DELTAW,IFWAVE
      COMMON /XABUND/XABUND(99),WTMOLE
      COMMON /XNF/XNFC(kw,6),XNFN(kw,6),XNFO(kw,6),XNFNE(kw,6),
     1            XNFMG(kw,6),XNFSI(kw,6),XNFS(kw,6),XNFFE(kw,5)
      COMMON /XNFP/XNFPC(kw,4),XNFPN(kw,5),XNFPO(kw,6),XNFPNE(kw,6),
     1             XNFPAL(kw,1),XNFPMG(kw,2),XNFPSI(kw,2),XNFPCA(kw,2),
     2             XNFPFE(kw,1),XNFPCH(kw),XNFPOH(kw)
      DIMENSION PART(kw,6)
C
      COMMON /CONTIN/CONTIN
      COMMON /CURVE/MINLOG,DABLOG,NABLOG
      REAL*8 MINLOG
  
      COMMON /GAM/GLOGR,GLOGS,GLOGW
      COMMON /IFPROF/IFPROF
      COMMON/LINDAT/WL,GFLOG,XJ,E,XJP,EP,CODE,GAMMAR,GAMMAS,GAMMAW,
     *ELO,GF,X1,X2,REF,OTHER1(3),OTHER2(3),LABEL(3),LABELP(3),NELION,       
     1               NBLO,NBUP,ISO1,ISO2           
      REAL*8 WL,E,EP
      COMMON /LINEY/LINE0(10,kw),DOPWAV(10,kw),ADAMP(10,kw)
      REAL*8 LINE0
  
      COMMON /OBS/EW,WLOBS,STAR(19)
      COMMON /RESIDU/RESIDU(15),TAUONE(15),DELWAV(15)
      COMMON /VTS/VTS(10),NVT

ChSt  diese strings (card,fff) werden beim Einlesen verwendet
ChSt  (SUBROUTINE READIN in readin.f).
	character*80 card,fff
C	character*80 card,fff*25
	character*60 file_m,id_ln*4,dumb*29
	integer xl(8)
	real*8 elmi(99)
      DATA elmi/   2HH , 2HHe,
     1 2HLi, 2HBe, 2HB , 2HC , 2HN , 2HO , 2HF , 2HNe, 2HNa, 2HMg,
     2 2HAl, 2HSi, 2HP , 2HS , 2HCl, 2HAr, 2HK , 2HCa, 2HSc, 2HTi,
     3 2HV , 2HCr, 2HMn, 2HFe, 2HCo, 2HNi, 2HCu, 2HZn, 2HGa, 2HGe,
     4 2HAs, 2HSe, 2HBr, 2HKr, 2HRb, 2HSr, 2HY , 2HZr, 2HNb, 2HMo,
     5 2HTc, 2HRu, 2HRh, 2HPd, 2HAg, 2HCd, 2HIn, 2HSn, 2HSb, 2HTe,
     6 2HI , 2HXe, 2HCs, 2HBa, 2HLa, 2HCe, 2HPr, 2HNd, 2HPm, 2HSm,
     7 2HEu, 2HGd, 2HTb, 2HDy, 2HHo, 2HEr, 2HTm, 2HYb, 2HLu, 2HHf,
     8 2HTa, 2HW , 2HRe, 2HOs, 2HIr, 2HPt, 2HAu, 2HHg, 2HTl, 2HPb,
     9 2HBi, 2HPo, 2HAt, 2HRn, 2HFr, 2HRa, 2HAc, 2HTh, 2HPa, 2HU ,
     T 2HNp, 2HPu, 2HAm, 2HCm, 2HBk, 2HCf, 2HEs/

!      EXP10(X)=EXP(X*2.30258509299405E0)

ChSt  igetarg und getarg sind C-Routinen.
ChSt  funktioniert deßwegen, da width9 mit: > width9mf inputfile
ChSt  gestartet wird.
	mm=iargc()
	call getarg(1,fff)
	open(1,file=fff,status='unknown')
	k=index(fff,'.')
	k=k
	file_m='width.out'
	open(6,file=file_m,status='unknown')
	file_m=fff(1:k)//'abn'
	open(3,file=file_m,status='unknown')

      IFPROF=0
      ITEMP=0
      NABLOG=0
      IFSURF=1
      MODEAV=1

	  do i=1,3
	   read(1,'((a))')card
	  enddo
	
	nlin=0
1	read(1,'((a))',end=2)card
	if(card(1:1).ne."'")then
	 go to 2
	else
ChSt    über nlin wird die Anzahl der Linien im inputfile ermittelt.
	nlin=nlin+1
	k=index(card,"',")
	if(k.ne.0)then
        file_m=' '
	file_m(1:k-2)=card(2:k-1)
	endif
	go to 1
	endif
	
2	nlin=nlin-1
	rewind 1
	 nnnn=nlin+4
	  do i=1,nnnn
	   read(1,'((a))')card
	  enddo
	 read(1,*)nvt,(vts(i),i=1,nvt)
	  read(1,*)Bmean
	  Bmean=Bmean*1000.
	
	open(5,file=file_m,status='unknown')

ChSt  Modellatmosphäre wird eingelesen
    3 CALL READIN(20)
	dumb(1:29)=' '
      write(3,"(10x,'Teff=',f7.0,' lg g=',f6.2,a29)")teff,glog,dumb
	  write(3,'(63a1)')(title(i),i=1,63)
	close (5)
	rewind 1
	  do i=1,3
	   read(1,'((a))')card
	  enddo

ChSt  Bestimmung des physikalischen Zustandes der Sternatmosphäre

      ITEMP=ITEMP+1
      IFOP(15)=0
      IF(IFPRES.NE.0)CALL POPS(0.D0,1,XNE)
			
      CALL POPS(1.01D0,11,XNFPH)
      CALL POPS(2.02D0,11,XNFPHE)
      CALL POPS(6.03D0,11,XNFPC)
      CALL POPS(7.04D0,11,XNFPN)
      CALL POPS(8.05D0,11,XNFPO)
      CALL POPS(10.05D0,11,XNFPNE)
      CALL POPS(12.01D0,11,XNFPMG)
      CALL POPS(13.00D0,11,XNFPAL)
      CALL POPS(14.01D0,11,XNFPSI)
      CALL POPS(20.01D0,11,XNFPCA)
      CALL POPS(26.00D0,11,XNFPFE)
      IF(IFMOL.EQ.0)THEN
      CALL POPS(1.01D0,12,XNFH)
      CALL POPS(2.02D0,12,XNFHE)
      CALL POPS(6.05D0,12,XNFC)
      CALL POPS(7.05D0,12,XNFN)
      CALL POPS(8.05D0,12,XNFO)
      CALL POPS(10.05D0,12,XNFNE)
      CALL POPS(12.05D0,12,XNFMG)
      CALL POPS(14.05D0,12,XNFSI)
      CALL POPS(16.05D0,12,XNFS)
      CALL POPS(26.04D0,12,XNFFE)
      ENDIF
      IF(IFMOL.EQ.1)THEN
      CALL POPS(106.00D0,11,XNFPCH)
      CALL POPS(108.00D0,11,XNFPOH)
!      CALL W(6HXNFPCH,XNFPCH,NRHOX)
!      CALL W(6HXNFPOH,XNFPOH,NRHOX)
C     THE POPS WILL NOT RETURN NUMBER DENSITIES WHEN MOLECULES ARE ON
C     SO WE COMPUTE NUMBER DENSITIES/PART FUNCTIONS  AND PART FUNCTIONS
      CALL POPS(6.05D0,11,XNFC)
      CALL POPS(7.05D0,11,XNFN)
      CALL POPS(8.05D0,11,XNFO)
      CALL POPS(10.05D0,11,XNFNE)
      CALL POPS(12.05D0,11,XNFMG)
      CALL POPS(14.05D0,11,XNFSI)
      CALL POPS(16.05D0,11,XNFS)
      CALL POPS(26.04D0,11,XNFFE)
      DO 444 J=1,NRHOX
      CALL PFSAHA(J,1,1,3,PART)
      XNFH(J,1)=XNFPH(J,1)*PART(J,1)
      XNFH(J,2)=XNFPH(J,2)
      CALL PFSAHA(J,2,2,13,PART)
      XNFHE(J,1)=XNFPHE(J,1)*PART(J,1)
      XNFHE(J,2)=XNFPHE(J,2)*PART(J,2)
      XNFHE(J,3)=XNFPHE(J,3)
      CALL PFSAHA(J,6,6,13,PART)
      XNFC(J,1)=XNFC(J,1)*PART(J,1)
      XNFC(J,2)=XNFC(J,2)*PART(J,2)
      XNFC(J,3)=XNFC(J,3)*PART(J,3)
      XNFC(J,4)=XNFC(J,4)*PART(J,4)
      XNFC(J,5)=XNFC(J,5)*PART(J,5)
      XNFC(J,6)=XNFC(J,6)*PART(J,6)
      CALL PFSAHA(J,7,6,13,PART)
      XNFN(J,1)=XNFN(J,1)*PART(J,1)
      XNFN(J,2)=XNFN(J,2)*PART(J,2)
      XNFN(J,3)=XNFN(J,3)*PART(J,3)
      XNFN(J,4)=XNFN(J,4)*PART(J,4)
      XNFN(J,5)=XNFN(J,5)*PART(J,5)
      XNFN(J,6)=XNFN(J,6)*PART(J,6)
      CALL PFSAHA(J,8,6,13,PART)
      XNFO(J,1)=XNFO(J,1)*PART(J,1)
      XNFO(J,2)=XNFO(J,2)*PART(J,2)
      XNFO(J,3)=XNFO(J,3)*PART(J,3)
      XNFO(J,4)=XNFO(J,4)*PART(J,4)
      XNFO(J,5)=XNFO(J,5)*PART(J,5)
      XNFO(J,6)=XNFO(J,6)*PART(J,6)
      CALL PFSAHA(J,10,6,13,PART)
      XNFNE(J,1)=XNFNE(J,1)*PART(J,1)
      XNFNE(J,2)=XNFNE(J,2)*PART(J,2)
      XNFNE(J,3)=XNFNE(J,3)*PART(J,3)
      XNFNE(J,4)=XNFNE(J,4)*PART(J,4)
      XNFNE(J,5)=XNFNE(J,5)*PART(J,5)
      XNFNE(J,6)=XNFNE(J,6)*PART(J,6)
      CALL PFSAHA(J,12,6,13,PART)
      XNFMG(J,1)=XNFMG(J,1)*PART(J,1)
      XNFMG(J,2)=XNFMG(J,2)*PART(J,2)
      XNFMG(J,3)=XNFMG(J,3)*PART(J,3)
      XNFMG(J,4)=XNFMG(J,4)*PART(J,4)
      XNFMG(J,5)=XNFMG(J,5)*PART(J,5)
      XNFMG(J,6)=XNFMG(J,6)*PART(J,6)
      CALL PFSAHA(J,14,6,13,PART)
      XNFSI(J,1)=XNFSI(J,1)*PART(J,1)
      XNFSI(J,2)=XNFSI(J,2)*PART(J,2)
      XNFSI(J,3)=XNFSI(J,3)*PART(J,3)
      XNFSI(J,4)=XNFSI(J,4)*PART(J,4)
      XNFSI(J,5)=XNFSI(J,5)*PART(J,5)
      XNFSI(J,6)=XNFSI(J,6)*PART(J,6)
      CALL PFSAHA(J,16,6,13,PART)
      XNFS(J,1)=XNFS(J,1)*PART(J,1)
      XNFS(J,2)=XNFS(J,2)*PART(J,2)
      XNFS(J,3)=XNFS(J,3)*PART(J,3)
      XNFS(J,4)=XNFS(J,4)*PART(J,4)
      XNFS(J,5)=XNFS(J,5)*PART(J,5)
      XNFS(J,6)=XNFS(J,6)*PART(J,6)
      CALL PFSAHA(J,26,5,13,PART)
      XNFFE(J,1)=XNFFE(J,1)*PART(J,1)
      XNFFE(J,2)=XNFFE(J,2)*PART(J,2)
      XNFFE(J,3)=XNFFE(J,3)*PART(J,3)
      XNFFE(J,4)=XNFFE(J,4)*PART(J,4)
      XNFFE(J,5)=XNFFE(J,5)*PART(J,5)
  444 CONTINUE
      ENDIF

ChSt  Nun beginnt die Häufigkeitsbestimmung

	call integ(rhox,abross,tauros,nrhox,abross(1)*rhox(1))

ChSt    Start der Hauptschleife mit dem Einlesen der Linieninformation
ChSt    Achtung alle Subroutinen der anderen Teile von width9 können
ChSt    auf diese Variablen zugreifen. Zwischen Groß- und Kleinschreibung
ChSt    wird nicht unterschieden. Beispiel:
ChSt     OBEWLG= LOG10(EW) ... in width9sb.f Zeile 38

	do ilin=1,nlin
	  read(1,*)id_ln,wla,expot,vtdummy,gflog,gammar,gammas,gammaw
     *   ,glande,ewma
ChSt    glande=99.0 -> nicht bekannt für diese Linie
	   if(glande.gt.90.)glande=0.
	  vmag=4.66e-13*2.997925e5*wla*Bmean*glande
	ew=ewma/10.
	  wl=wla/10.
	   elo=expot/1.23981e-4
	  read(id_ln,'(a2,i2)')el_id,ionn
	    do i=1,99
!	     if(el_id.eq.elem(i))go to 1111
	     if(el_id.eq.elem(i).or.el_id.eq.elmi(i))go to 1111
	    enddo
1111	nnn=i
	code=float(nnn)+float(ionn-1)/100.
	nelion=(nnn-1)*6+ionn

140	format(/,f10.3,f8.3,f14.3,f8.2,i4,3f8.2)
      GFLOG=GFLOG+X1+X2
ChSt  2.30258509299405 ist ln(10); 10^X gibts in F77 nicht.
ChSt  GF=EXP(GFLOG*2.30258509299405E0) entspricht GF=10^GFLOG
      GF=EXP(GFLOG*2.30258509299405E0)
                                                      
      GAMMAR=EXP(GAMMAR*2.30258509299405E0)
      GAMMAS=EXP(GAMMAS*2.30258509299405E0)
      GAMMAW=EXP(GAMMAW*2.30258509299405E0)

C     CLASSICAL DAMPING CONSTANT
      IF(GAMMAR.EQ.1.)GAMMAR=2.223E13/WL**2

      IF(GAMMAS.NE.1..AND.GAMMAW.NE.1.)GO TO 141
      IF(GAMMAS.NE.1.)GO TO 138
      IF(CODE.GE.100.)GO TO 137

	epup=1./(wl*1.e-7)+elo
      ZEFF=(CODE-FLOAT( INT(CODE)))*100.+1.
      EFFNSQ=25.
      DEPUP=POTION(NELION)-EPUP*1.23981e-4
      IF(DEPUP.GT.0.)EFFNSQ=13.595*ZEFF**2/DEPUP
      GAMMAS=1.0E-8*EFFNSQ**2*SQRT(EFFNSQ)
      GO TO 138
  137 GAMMAS=1.0E-5

  138 IF(GAMMAW.NE.1.)GO TO 141
      IF(CODE.GE.100.)GO TO 139
	epup=1./(wl*1.e-7)+elo
      ZEFF=(CODE-FLOAT( INT(CODE)))*100.+1.
      EFFNSQ=25.
      DEPUP=POTION(NELION)-EPUP*1.23981e-4
      DELO=POTION(NELION)-ELO*1.23981e-4
      IF(DEPUP.GT.0.)EFFNSQ=(13.595*ZEFF**2/DEPUP)**2-
     * (13.595*ZEFF**2/DELO)**2
      RSQ=2.5*EFFNSQ/ZEFF**2
      NSEQ=CODE-ZEFF+1.
c (Bikmaev)
c      IF(NSEQ.GT.20.AND.NSEQ.LT.29)RSQ=(45.-FLOAT(NSEQ))/ZEFF
      GAMMAW=4.5E-9*RSQ**.4
c      write(3,*) GAMMAW
      GO TO 141
  139 GAMMAW=1.E-7/ZEFF
                                                     
  141 GLOGR= LOG10(GAMMAR)
      GLOGS= LOG10(GAMMAS)
      GLOGW= LOG10(GAMMAW)

	write(6,140)wl,gflog,elo,code,nelion,glogr,glogs,glogw

      FREQ=2.997925E17/WL
      FREQLG= LOG(FREQ)
      DO 61 J=1,NRHOX
      EHVKT(J)=EXP(-FREQ*HKT(J))
      STIM(J)=1.-EHVKT(J)
      BNU(J)=1.47439E-2*(FREQ/1.E15)**3*EHVKT(J)/STIM(J)
   61 CONTINUE
C     CALCULATES CONTINUOUS OPACITIES
      CALL KAPP(1,NSTEPS,STEPWT)
C     CALCULATES THE SOURCE FUNCTION AND SURFACE FLUX OR SURFACE INTENSITY
      CALL JOSH(IFSCAT,IFSURF)
      IF(IFSURF.LT.2)CONTIN=HNU(1)
      IF(IFSURF.EQ.2)CONTIN=SURFI(1)
C     TREAT HYDROGEN LINES AS CONTINUUM
      IF(IFSCAT.EQ.0.)GO TO 162
      DO 63 J=1,NRHOX
      SCONT(J)=(ACONT(J)*SCONT(J)+ALINE(J)*SLINE(J))/(ACONT(J)+ALINE(J))
      ACONT(J)=ACONT(J)+ALINE(J)
      SLINE(J)=BNU(J)
   63 CONTINUE
      GO TO 164
C     SOURCE FUNCTION WILL BE WEIGHTED MEAN OF CONTINUUM SOURCE
C        FUNCTION (INCLUDING JNU) AND BNU
  162 DO 163 J=1,NRHOX
      SCONT(J)=SNU(J)
      SLINE(J)=BNU(J)
  163 ACONT(J)=ACONT(J)+ALINE(J)
  164 CONTINUE
C     CALCULATES THE LINE CENTER MASS ABSORPTION COEFFICIENT FOR UNIT ABUNDANCE

      CALL LINCEN(ABUND1,VTS,NVT)
      DO 65 IVT=1,NVT
C     CALCULATES THE CURVE OF GROWTH AND THE FINAL ABUNDANCE
C     FIRST GUESS FOR ABLG
      ABLG= LOG10(ABUND1)
      IF(IFPROF.EQ.1)WRITE(6,64)
   64 FORMAT(///23H PROFILES FOR NEXT LINE)
      CALL COG(IVT,ABLG,TAULG,xl,yy,n1)
C     CALCULATES THE AVERAGE ABUNDANCE
	write(3,1400)id_ln,wla,expot,vts(ivt),gflog,glande,
     *  ewma,taulg,ablg
1400    format(1x,a4,f11.4,f9.4,f5.2,f8.3,4f7.2)
65	continue
	enddo
	close(6,status='delete')
	stop 'Best regards !'
      END
