192 lines
5.2 KiB
Fortran
192 lines
5.2 KiB
Fortran
SUBROUTINE HYDINI
|
|
C
|
|
C Initializes necessary arrays for evaluating hydrogen line profiles
|
|
C from the Lemke, Tremblay-Bergeron, or Schoening-Butler tables
|
|
C
|
|
INCLUDE 'PARAMS.FOR'
|
|
INCLUDE 'MODELP.FOR'
|
|
c DIMENSION WLINE(4,22)
|
|
DIMENSION IILW(100),IIUP(100)
|
|
CHARACTER*1 CHAR
|
|
DATA INIT /0/
|
|
C
|
|
IF(INIT.EQ.0) THEN
|
|
DO I=1,4
|
|
DO J=I+1,22
|
|
CALL STARK0(I,J,IZZ,XK,WL0,FIJ,FIJ0)
|
|
WLINE(I,J)=WL0
|
|
c OSCH(I,J)=FIJ+FIJ0
|
|
END DO
|
|
END DO
|
|
INIT=1
|
|
END IF
|
|
DO I=1,4
|
|
DO J=1,22
|
|
ILIN0(I,J)=0
|
|
END DO
|
|
END DO
|
|
C
|
|
C --------------------------------------------
|
|
C Schoening-Butler tables - for IHYDPR < 0
|
|
C --------------------------------------------
|
|
C
|
|
IF(IHYDPR.LT.0) THEN
|
|
IHYDPR=67
|
|
ILEMKE=0
|
|
NLINE=12
|
|
c
|
|
OPEN(UNIT=IHYDPR,FILE='./data/hydprf.dat',STATUS='OLD')
|
|
write(6,*) ' reading Schoening-Butler tables'
|
|
C
|
|
DO I=1,12
|
|
READ(IHYDPR,500)
|
|
END DO
|
|
DO 100 ILINE=1,NLINE
|
|
C
|
|
C read the tables, which have to be stored in file
|
|
C unit IHYDPR (which is the input parameter in the progarm)
|
|
C
|
|
READ(IHYDPR,501) I,J
|
|
IF(ILINE.EQ.12) J=10
|
|
WL0=WLINE(I,J)
|
|
ILIN0(I,J)=ILINE
|
|
READ(IHYDPR,*) CHAR,NWL,(WL(I,ILINE),I=1,NWL)
|
|
READ(IHYDPR,*) CHAR,NT,(XT(I,ILINE),I=1,NT)
|
|
READ(IHYDPR,*) CHAR,NE,(XNE(I,ILINE),I=1,NE)
|
|
READ(IHYDPR,500)
|
|
NWLH(ILINE)=NWL
|
|
NWLHYD(ILINE)=NWL
|
|
NTH(ILINE)=NT
|
|
NEH(ILINE)=NE
|
|
C
|
|
DO I=1,NWL
|
|
IF(WL(I,ILINE).LT.1.E-4) WL(I,ILINE)=1.E-4
|
|
WLHYD(ILINE,I)=LOG10(WL(I,ILINE))
|
|
END DO
|
|
C
|
|
DO IE=1,NE
|
|
DO IT=1,NT
|
|
READ(IHYDPR,500)
|
|
READ(IHYDPR,*) (PRF(IWL,IT,IE,ILINE),IWL=1,NWL)
|
|
END DO
|
|
END DO
|
|
C
|
|
C coefficient for the asymptotic profile is determined from
|
|
C the input data
|
|
C
|
|
XCLOG=PRF(NWL,1,1,ILINE)+2.5*LOG10(WL(NWL,ILINE))+31.5304-
|
|
* XNE(1,ILINE)-2.*LOG10(WL0)
|
|
XKLOG=0.6666667*(XCLOG-0.176)
|
|
XK=EXP(XKLOG*2.3025851)
|
|
C
|
|
DO ID=1,ND
|
|
C
|
|
C temperature is modified in order to account for the
|
|
C effect of turbulent velocity on the Doppler width
|
|
C
|
|
T=TEMP(ID)+6.06E-9*VTURB(ID)
|
|
ANE=ELEC(ID)
|
|
TL=LOG10(T)
|
|
ANEL=LOG10(ANE)
|
|
F00=1.25E-9*ANE**0.666666667
|
|
FXK=F00*XK
|
|
DOP=1.E8/WL0*SQRT(1.65E8*T)
|
|
DBETA=WL0*WL0/2.997925E18/FXK
|
|
BETAD=DBETA*DOP
|
|
C
|
|
C interpolation to the actual values of temperature and electron
|
|
C density. The result is stored at array PRFHYD, having indices
|
|
C ILINE (line number: 1 for L-alpha,..., 4 for H-delta, etc.);
|
|
C 5 for H-alpha,..., 8 for H-delta, etc.)
|
|
C ID - depth index
|
|
C IWL - wavelength index
|
|
C
|
|
DO IWL=1,NWL
|
|
CALL INTHYD(PROF,TL,ANEL,IWL,ILINE)
|
|
PRFHYD(ILINE,ID,IWL)=PROF
|
|
END DO
|
|
END DO
|
|
100 CONTINUE
|
|
CLOSE(IHYDPR)
|
|
C
|
|
500 FORMAT(1X)
|
|
501 FORMAT(12X,I1,9X,I1)
|
|
C
|
|
IHYDPR=-IHYDPR
|
|
RETURN
|
|
END IF
|
|
C
|
|
C ---------------------------------
|
|
C read Lemke or Tremblay tables
|
|
C ---------------------------------
|
|
C
|
|
if(ihydpr.lt.20) ihydpr=ihydpr+20
|
|
if(ihydpr.eq.21) then
|
|
open(unit=ihydpr,file='./data/lemke.dat',status='old')
|
|
write(6,641) ihydpr
|
|
else if(ihydpr.eq.22) then
|
|
open(unit=ihydpr,file='./data/tremblay.dat',status='old')
|
|
write(6,642) ihydpr
|
|
end if
|
|
641 format(' -----------'/
|
|
* ' reading Lemke tables; ihydpr =',i3,/
|
|
* ' -----------')
|
|
642 format(' -----------'/
|
|
* ' reading Tremblay tables; ihydpr =',i3,/
|
|
* ' -----------')
|
|
C
|
|
ILEMKE=1
|
|
READ(IHYDPR,*) NTAB
|
|
write(6,611) ntab
|
|
611 format(' ntab',i4)
|
|
DO ITAB=1,NTAB
|
|
ILINEB=ILINE
|
|
READ(IHYDPR,*) NLLY
|
|
DO ILI=1,NLLY
|
|
ILINE=ILINE+1
|
|
READ(IHYDPR,*) I,J,ALMIN,ANEMIN,TMIN,DLA,DLE,DLT,
|
|
* NWL,NE,NT
|
|
WL0=WLINE(I,J)
|
|
ILIN0(I,J)=ILINE
|
|
NWLH(ILINE)=NWL
|
|
NWLHYD(ILINE)=NWL
|
|
NTH(ILINE)=NT
|
|
NEH(ILINE)=NE
|
|
iilw(iline)=i
|
|
iiup(iline)=j
|
|
DO IWL=1,NWL
|
|
WL(IWL,ILINE)=ALMIN+(IWL-1)*DLA
|
|
WLHYD(ILINE,IWL)=WL(IWL,ILINE)
|
|
WL(IWL,ILINE)=EXP(2.3025851*WL(IWL,ILINE))
|
|
END DO
|
|
DO INE=1,NE
|
|
XNE(INE,ILINE)=ANEMIN+(INE-1)*DLE
|
|
END DO
|
|
DO IT=1,NT
|
|
XT(IT,ILINE)=TMIN+(IT-1)*DLT
|
|
END DO
|
|
END DO
|
|
c
|
|
DO ILI=1,NLLY
|
|
ILNE=ILINEB+ILI
|
|
NWL=NWLH(ILNE)
|
|
READ(IHYDPR,500)
|
|
DO INE=1,NEH(ILNE)
|
|
DO IT=1,NTH(ILNE)
|
|
READ(IHYDPR,*) QLT,(PRF(IWL,IT,INE,ILNE),IWL=1,NWL)
|
|
END DO
|
|
END DO
|
|
C
|
|
i=iilw(ilne)
|
|
j=iiup(ilne)
|
|
DO ID=1,ND
|
|
CALL HYDTAB(I,J,ID)
|
|
END DO
|
|
END DO
|
|
END DO
|
|
NLIHYD=ILNE
|
|
CLOSE(IHYDPR)
|
|
C
|
|
RETURN
|
|
END
|