457 lines
14 KiB
Fortran
457 lines
14 KiB
Fortran
SUBROUTINE INIBL0
|
|
C
|
|
C AUXILIARY INITIALIZATION PROCEDURE
|
|
C
|
|
INCLUDE 'PARAMS.FOR'
|
|
INCLUDE 'MODELP.FOR'
|
|
INCLUDE 'LINDAT.FOR'
|
|
INCLUDE 'SYNTHP.FOR'
|
|
INCLUDE 'WINCOM.FOR'
|
|
parameter (un=1.)
|
|
character*2 iu
|
|
character*6 ilab
|
|
DIMENSION CROSS(MCROSS,MFRQ),
|
|
* ABSO(MFREQ),EMIS(MFREQ),SCAT(MFREQ),
|
|
* ABSOC(MFREQC),EMISC(MFREQC),SCATC(MFREQC)
|
|
COMMON/LIMPAR/ALAM0,ALAM1,FRMIN,FRLAST,FRLI0,FRLIM
|
|
COMMON/BLAPAR/RELOP,SPACE0,CUTOF0,TSTD,DSTD,ALAMC
|
|
common/lasers/lasdel
|
|
common/linrej/ilne(mdepth),ilvi(mdepth)
|
|
common/velaux/velmax,iemoff,nltoff,itrad
|
|
common/alsave/ALAM0s,ALASTs,CUTOF0s,CUTOFSs,RELOPs,SPACEs
|
|
C
|
|
C --------------------------------------------------------------
|
|
C Parameters controlling an evaluation of the synthetic spectrum
|
|
C
|
|
C --------------------------------------------------------------
|
|
C
|
|
C ALAM0, ALAM1 - synthetic spectrum is evaluated between wavelengths
|
|
C ALAM0 (initial) and ALAM1 (final), given in Anstroms
|
|
C CUTOF0 - cutoff parameter for normal lines (given in Angstroms)
|
|
C ie the maximum distance from the line center, in
|
|
C which the opacity in the line is allowd to contribute
|
|
C to the total opacity (recommended 5 - 10)
|
|
C CUTOFS = SPACON
|
|
C SPACON - spacing of the continuum wavelength points
|
|
C (at the midpoint of teh total interval; actual spacing
|
|
C is equidistant in log(lambda)
|
|
C RELOP - the minimum value of the ratio (opacity in the line
|
|
C center)/(opacity in continuum), for which is the line
|
|
C taken into account (usually 1d-4 to 1d-3)
|
|
C SPACE - the maximum distance of two neighbouring frequency
|
|
C points for evaluating the spectrum; in Angstroms
|
|
C
|
|
C INLTE = 0 - pure LTE (no line in NLTE)
|
|
C ne.0 - NLTE option, ie one or more lines treated
|
|
C in the exact or approximate NLTE approach
|
|
C IFHE2 gt.0 - He II line opacity in the first four series
|
|
C (Lyman, Balmer, Paschen, Brackett)
|
|
C for lines with lambda < 3900 A
|
|
C is taken into account even if line list
|
|
C does not contain any He II lines (i.e.
|
|
C He II lines are treated as the hydrogen lines)
|
|
C
|
|
C IHYDPR = 0 - means that hydrogen lines Stark profiles
|
|
C are calculated by approximate formulae
|
|
C > 0 - hydrogen lines Stark profiles are calculated
|
|
C in detail, using the Schoening & Butler tables;
|
|
C (for 1-2 to 1-5; 2-3 to 2-10).
|
|
C the tables are stored in file FOR0xx.dat,
|
|
C where xx=IHYDPR;
|
|
C higher Balmer lines are calculated as before
|
|
C
|
|
C the meaning of other parameters is quite analogous, for the
|
|
C following lines
|
|
C
|
|
C IHE1PR - He I lines at 4471, 4026, 4387, and 4922 Angstroms
|
|
C (tables calculated by Barnard, Cooper, and Shamey)
|
|
C IHE2PR - for the He II lines calculated by Schoening and Butler,
|
|
C
|
|
if(ifeos.le.0) then
|
|
READ(55,*) IFREQ,INLTE,ICONTL,INLIST,IFHE2
|
|
IF(LTE) INLTE=0
|
|
READ(55,*) IHYDPR,IHE1PR,IHE2PR
|
|
READ(55,*) ALAM0,ALAST,CUTOF0,CUTOFS,RELOP,SPACE
|
|
end if
|
|
C
|
|
IF(IDSTD.EQ.0) THEN
|
|
ID1=5
|
|
NDSTEP=(ND-2*ID1)/2
|
|
IDSTD=2*ND/3
|
|
ELSE IF(IDSTD.LT.0) THEN
|
|
ID1=1
|
|
NDSTEP=-IDSTD
|
|
IDSTD=2*ND/3
|
|
END IF
|
|
if(imode.le.-3) ndstep=1
|
|
c
|
|
alam0s=alam0
|
|
alasts=alast
|
|
cutof0s=cutof0
|
|
cutofss=cutofs
|
|
relops=relop
|
|
spaces=space
|
|
C
|
|
C if ALAST.lt.0 - set up vacuum wavelengths everywhere
|
|
C
|
|
vaclim=2000.
|
|
if(alast.lt.0.) then
|
|
alast=abs(alast)
|
|
alasts=alast
|
|
vaclim=1.e18
|
|
end if
|
|
c
|
|
if(inlte.lt.10) then
|
|
lasdel=.true.
|
|
else if(inlte.le.20) then
|
|
inlte=inlte-10
|
|
lasdel=.false.
|
|
else if(inlte.le.30) then
|
|
inlte=inlte-20
|
|
ifreq=11
|
|
lasdel=.true.
|
|
else if(inlte.le.40) then
|
|
inlte=inlte-30
|
|
ifreq=11
|
|
lasdel=.false.
|
|
end if
|
|
C
|
|
ibin(0)=mod(inlist,10)
|
|
do ilist=1,mmlist
|
|
tmlim(ilist)=tmolim
|
|
ibin(ilist)=mod(inlist,10)
|
|
ivdwli(ilist)=0
|
|
iun=19+ilist
|
|
write(iu,622) iun
|
|
622 format(i2)
|
|
amlist(ilist)='fort.' // iu
|
|
end do
|
|
c
|
|
if(imode.ge.-3.and.imode.le.1) then
|
|
nmlist=0
|
|
numlis=0
|
|
read(55,*,err=5,end=5) nmlist,(iunitm(ilist),ilist=1,nmlist)
|
|
do ilist=1,nmlist
|
|
write(iu,622) iunitm(ilist)
|
|
amlist(ilist) ='fort.' // iu
|
|
end do
|
|
5 continue
|
|
c
|
|
ilist=0
|
|
amlist(0)='fort.19'
|
|
read(3,*,err=20,end=20) amlist(0),ibin(0)
|
|
c
|
|
ilist=0
|
|
10 continue
|
|
ilist=ilist+1
|
|
read(3,*,end=20) amlist(ilist),ibin(ilist),tmlim(ilist)
|
|
numlis=numlis+1
|
|
go to 10
|
|
20 continue
|
|
if(numlis.gt.0) nmlist=numlis
|
|
if(nmlist.gt.0.and.ifmol.eq.0) then
|
|
write(*,*) 'NEEDS TO SET IFMOL > 0 with NMLIST>0'
|
|
stop
|
|
end if
|
|
c
|
|
ilist=0
|
|
ilab='ATOMIC'
|
|
write(6,623) ilist,ilab,trim(amlist(ilist)),ibin(ilist)
|
|
ilab='MOLEC '
|
|
do ilist=1,nmlist
|
|
write(6,624) ilist,ilab,trim(amlist(ilist)),ibin(ilist),
|
|
* tmlim(ilist)
|
|
end do
|
|
623 format(/'************************'/
|
|
* ' LINE LISTS:'/
|
|
* /' ILIST',8x,'FILENAME IBIN TMLIM'/
|
|
* i4,2x,a6,2x,a,2x,i4,f11.1)
|
|
624 format( i4,2x,a6,2x,a,2x,i4,f11.1)
|
|
end if
|
|
c
|
|
C
|
|
c VTB - turbulent velocity (in km/s). In non-negative, this
|
|
C value overwrites the value given by the standard input
|
|
C
|
|
read(55,*,err=30,end=30) VTB
|
|
if(ifwin.le.0) then
|
|
if(vtb.ge.0.) then
|
|
WRITE(6,608) VTB
|
|
608 FORMAT(//' TURBULENT VELOCITY - CHANGED TO VTURB =',
|
|
* 1PE10.3,' KM/S'/' ------------------'/)
|
|
do id=1,nd
|
|
vturb(id)=vtb*vtb*1.e10
|
|
end do
|
|
end if
|
|
end if
|
|
C
|
|
TSTD=TEMP(IDSTD)
|
|
VTS=VTURB(IDSTD)
|
|
DSTD=SQRT(1.4E7*TSTD+VTS)
|
|
30 continue
|
|
C
|
|
C angle points (in case the specific intensities are evaluated
|
|
C
|
|
C NMU0 - number of angles:
|
|
C >0 - and if also ANG0>0, angles (mu's) equidistant
|
|
C between 1 and ANG0
|
|
C >0 - and if also ANG0<0, angles (mu's) equidistant
|
|
C between 0.7 and ANG0, and sinuses equidistatnt for
|
|
C others
|
|
C <0 - angles read in the next record
|
|
C ANG0 - minimum mu (see above)
|
|
C IFLUX - mode for evaluating angle-dependent intensities and
|
|
C the corresponding flux:
|
|
C =0 - no specifiec intensities are evaluated; only usual
|
|
C flux is stored (unit 7 and 17)
|
|
C =1 - specific intensities are evaluated;
|
|
C and stored on unit 18
|
|
C =2 - (interesting only for the case of macroscopic
|
|
C velocity field); specific intensities evaluated by
|
|
C a simple formal solution (RESOLV)
|
|
C
|
|
NMU0=1
|
|
ANG0=1.
|
|
ANGL(1)=1.
|
|
WANGL(1)=0.
|
|
IFLUX=0
|
|
velmax=3.e5
|
|
nltoff=0
|
|
iemoff=0
|
|
itrad=0
|
|
do id=1,nd
|
|
wdil(id)=un
|
|
end do
|
|
if(ifwin.le.0) then
|
|
READ(55,*,end=100,err=100) NMU0,ANG0,IFLUX
|
|
C
|
|
C determinantion of the angle points and weights
|
|
C
|
|
IF(NMU0.LT.0) THEN
|
|
NMU0=IABS(NMU0)
|
|
READ(55,*) (ANGL(IMU),IMU=1,NMU0)
|
|
DO IMU=2,NMU0-1
|
|
WANGL(IMU)=0.5*(ANGL(IMU-1)+ANGL(IMU+1))
|
|
END DO
|
|
WANGL(1)=0.5*(ANGL(1)-ANGL(2))
|
|
WANGL(NMU0)=0.5*(ANGL(NMU0-1)-ANGL(NMU0))
|
|
ELSE
|
|
IF(ANG0.GT.0.) THEN
|
|
IF(NMU0.GT.1) THEN
|
|
DMU=(1.-ANG0)/(NMU0-1)
|
|
DO IMU=1,NMU0
|
|
ANGL(IMU)=1.-(IMU-1)*DMU
|
|
WANGL(IMU)=DMU
|
|
END DO
|
|
WANGL(1)=0.5*DMU
|
|
WANGL(NMU0-1)=0.5*DMU
|
|
WANGL(NMU0)=2.*DMU
|
|
END IF
|
|
ELSE
|
|
ANGH=0.70710678
|
|
DMU=ANGH/(NMU0-1)
|
|
DO IMU=1,NMU0
|
|
ANGL(IMU)=(IMU-1)*DMU
|
|
ANGL(IMU)=SQRT(1.-ANGL(IMU)**2)
|
|
IF(IMU.GT.1.AND.IMU.LT.NMU0)
|
|
* WANGL(IMU)=0.5*(ANGL(IMU-1)+ANGL(IMU+1))
|
|
END DO
|
|
WANGL(1)=0.5*(ANGL(1)-ANGL(2))
|
|
WANGL(NMU0)=0.5*(ANGL(NMU0-1)-ANGL(NMU0))
|
|
IF(ANG0.LT.0.) DMU=(ANGH+ANG0)/(NMU0-1)
|
|
DO IMU=1,NMU0-2
|
|
ANGL(IMU+NMU0)=ANGH-IMU*DMU
|
|
WANGL(IMU+NMU0)=DMU
|
|
END DO
|
|
WANGL(NMU0)=WANGL(NMU0)+0.5*DMU
|
|
WANGL(2*NMU0-3)=0.5*DMU
|
|
WANGL(2*NMU0-2)=2.*DMU
|
|
NMU0=2*NMU0-2
|
|
END IF
|
|
END IF
|
|
IF(NMU0.LE.0) GO TO 100
|
|
WRITE(6,609) NMU0,(ANGL(I),I=1,NMU0)
|
|
609 FORMAT(//' SPECIFIC INTENSITIES COMPUTED FOR',I3,
|
|
* ' ANGLES mu=cos(theta) ='/
|
|
* ' ---------------------------------',
|
|
* '------------------------'//
|
|
* (10F7.2))
|
|
100 CONTINUE
|
|
else
|
|
itrad=1
|
|
read(55,*,end=110,err=110) velmax,ITRAD,nltoff,iemoff
|
|
110 write(6,602) velmax,itrad,nltoff,iemoff
|
|
if(velmax.lt.0.) then
|
|
velmax=3.e5
|
|
go to 120
|
|
end if
|
|
602 format(//' velmax (velocity for line rejection)',
|
|
* ' itrad,nltoff,iemoff',f10.1,2i3)
|
|
C
|
|
C Set up rays and weights
|
|
C
|
|
call velset
|
|
call radtem
|
|
CALL SETRAY
|
|
CALL WGTJH1
|
|
C
|
|
end if
|
|
C
|
|
120 CONTINUE
|
|
velmax=velmax*1.e5
|
|
do id=1,nd
|
|
ilvi(id)=0
|
|
ilne(id)=0
|
|
if(vel(id).gt.velmax.and.iemoff.eq.0) ilvi(id)=1
|
|
if(vel(id).gt.velmax.and.nltoff.gt.0.and.iemoff.gt.0)
|
|
* ilne(id)=1
|
|
end do
|
|
C
|
|
IF(IMODE.EQ.-1) THEN
|
|
INLTE=0
|
|
CUTOF0=0.
|
|
END IF
|
|
C
|
|
C continuum frequencies
|
|
C
|
|
if(ifwin.le.0) then
|
|
alam0=alam0s
|
|
if(alam0s.eq.0.) alam0=5.e7/temp(1)/10.
|
|
if(alam0s.lt.0.) alam0=-5.e7/temp(1)/alam0s
|
|
alast=alasts
|
|
if(alasts.eq.0.) alast=5.e7/temp(1)*20.
|
|
if(alasts.lt.0.) alast=-5.e7/temp(1)*alasts
|
|
c if(alast.gt.1.e5) alast=1.e5
|
|
ALAMC=(ALAM0+ALAST)*0.5
|
|
if(space.eq.0.) space=4.3e-8*sqrt(temp(idstd))*alamc
|
|
if(space.lt.0.) space=-5.72e-8*sqrt(temp(idstd))*alamc*space
|
|
SPACF=2.997925E18/ALAMC/ALAMC*SPACE
|
|
WRITE(6,601) ALAM0,ALAST,CUTOF0,RELOP,SPACF,SPACE
|
|
CUTOF0=0.1*CUTOF0
|
|
SPACE0=SPACE*0.1
|
|
ALAM0=1.D-1*ALAM0
|
|
ALAST=1.D-1*ALAST
|
|
ALAMC=ALAMC*0.1
|
|
ALST00=ALAST
|
|
FRLAST=2.997925D17/ALAST
|
|
NFREQ=2
|
|
FREQ(1)=2.997925D17/ALAM0
|
|
FREQ(2)=FRLAST
|
|
C
|
|
else
|
|
C
|
|
spacon=cutofs
|
|
IF(SPACON.EQ.0) SPACON=3.
|
|
XFR=(ALAST-ALAM0)/SPACON
|
|
NFREQC=int(XFR)+1
|
|
NFREQC=MIN(NFREQC,MFREQC)
|
|
NFREQC=MAX(NFREQC,2)
|
|
DLAMLO=LOG10(ALAST/ALAM0)/(NFREQC-1)
|
|
AL0L=LOG10(ALAM0)
|
|
alambe=alam0
|
|
DO IJ=1,NFREQC
|
|
AL=AL0L+(IJ-1)*DLAMLO
|
|
ALAM=EXP(2.3025851*AL)
|
|
WLAMC(IJ)=ALAM
|
|
FREQC(IJ)=2.997925E18/ALAM
|
|
END DO
|
|
ALAMC=(ALAM0+ALAST)*0.5
|
|
SPACF=2.997925E18/ALAMC/ALAMC*SPACE
|
|
WRITE(6,601) ALAM0,ALAST,CUTOF0,RELOP,SPACF,SPACE
|
|
CUTOF0=0.1*CUTOF0
|
|
SPACE0=SPACE*0.1
|
|
ALAM0=1.D-1*ALAM0
|
|
ALAST=1.D-1*ALAST
|
|
ALAMC=ALAMC*0.1
|
|
ALST00=ALAST
|
|
FRLAST=2.997925D17/ALAST
|
|
NFREQ=2
|
|
FREQ(1)=2.997925D17/ALAM0
|
|
FREQ(2)=FRLAST
|
|
c
|
|
end if
|
|
c
|
|
CALL SIGAVS
|
|
IF(IHYDPR.NE.0) THEN
|
|
CALL HYDINI
|
|
CALL XENINI
|
|
END IF
|
|
IF(IHE1PR.GT.0) CALL HE1INI
|
|
IF(IHE2PR.GT.0) CALL HE2INI
|
|
C
|
|
C auxiliary quantities for dissolved fractions
|
|
C
|
|
DO ID=1,ND
|
|
CALL DWNFR0(ID)
|
|
CALL WNSTOR(ID)
|
|
END DO
|
|
C
|
|
c pretabulate expansion coefficients for the Voigt function
|
|
c
|
|
CALL PRETAB
|
|
c
|
|
c calculate the characteristic standard opacity
|
|
c
|
|
IF(IMODE.LE.2) THEN
|
|
if(ifwin.le.0.and.ndstep.eq.0) then
|
|
c
|
|
c old procedure
|
|
c
|
|
CALL CROSET(CROSS)
|
|
DO ID=1,ND
|
|
CALL OPAC(ID,CROSS,ABSO,EMIS,SCAT)
|
|
ABSTD(ID)=MIN(ABSO(1),ABSO(2))
|
|
END DO
|
|
else
|
|
c
|
|
c new procedure
|
|
c
|
|
if(ifwin.le.0) then
|
|
nfreqc=ifix(real(cutofs,4))
|
|
if(nfreqc.eq.0) nfreqc=mfreq
|
|
all0=log(alam0)
|
|
all1=log(alast)
|
|
dlc=(all1-all0)/(nfreqc-1)
|
|
do ijc=1,nfreqc
|
|
wlamc(ijc)=exp(all0+(ijc-1)*dlc)
|
|
freqc(ijc)=2.997925e17/wlamc(ijc)
|
|
end do
|
|
CALL CROSEW(CROSS)
|
|
do id=1,nd
|
|
CALL OPACON(ID,CROSS,ABSOC,EMISC,SCATC)
|
|
do ijc=1,nfreqc
|
|
abstdw(ijc,id)=absoc(ijc)
|
|
end do
|
|
end do
|
|
|
|
c write(*,*) 'abstdw(1,ij)',(abstdw(ij,1),ij=1,nfreqc)
|
|
c write(*,*) 'abstdw(50,ij)',(abstdw(ij,50),ij=1,nfreqc)
|
|
c
|
|
else
|
|
CALL CROSEW(CROSS)
|
|
DO ID=1,ND
|
|
CALL OPACW(ID,CROSS,ABSO,EMIS,ABSOC,EMISC,SCATC,0)
|
|
DO IJ=1,NFREQC
|
|
ABSTDW(IJ,ID)=ABSOC(IJ)/DENSCON(ID)
|
|
END DO
|
|
END DO
|
|
end if
|
|
end if
|
|
END IF
|
|
C
|
|
601 FORMAT(//'----------------------------------------------'/
|
|
* ' BASIC INPUT PARAMETERS FOR SYNTHETIC SPECTRA'/
|
|
* ' ---------------------------------------------'/
|
|
* ' INITIAL LAMBDA',28X,1H=,F10.3,' ANGSTROMS'/
|
|
* ' FINAL LAMBDA',28X,1H=,F10.3,' ANGSTROMS'/
|
|
* ' CUTOFF PARAMETER',26X,1H=,F10.3,' ANGSTROMS'/
|
|
* ' MINIMUM VALUE OF (LINE OPAC.)/(CONT.OPAC) =',1PE10.1/
|
|
* ' MAXIMUM FREQUENCY SPACING',17X,1H=,1PE10.3,' I.E.',
|
|
* 0PF6.3,' ANGSTROMS'/
|
|
* ' ---------------------------------------------'/)
|
|
c
|
|
write(6,612) idstd,ndstep
|
|
612 format(/'IDSTD, NDSTEP = ',2i5/)
|
|
RETURN
|
|
END
|