161 lines
5.3 KiB
Fortran
161 lines
5.3 KiB
Fortran
SUBROUTINE INPMOD
|
|
C =================
|
|
C
|
|
C Read an initial model atmosphere from unit 8
|
|
C File 8 contains:
|
|
C 1. NDPTH - number of depth points in which the initial model is
|
|
C given (if not equal to ND, routine interpolates
|
|
C automatically to the set DM by linear interpolation
|
|
C in log(DM)
|
|
C NUMPAR - number of input model parameters in each depth
|
|
C = 3 for LTE model - ie. N, T, N(electron);
|
|
C > 3 for NLTE model)
|
|
C 2. DEPTH(ID),ID=1,NDPTH - mass-depth points for the input model
|
|
C 3. for each depth:
|
|
C T - temperature
|
|
C ANE - electron density
|
|
C RHO - mass density
|
|
C level populations - only for NLTE input model
|
|
C Number of input level populations need not be
|
|
C equal to NLEVEL; in that case the procedure
|
|
C CHANGE is called from START to calculate the
|
|
C remaining level populations
|
|
C
|
|
C Note: The output file 7, which is created by this program
|
|
C (procedure OUTPUT) has the same structure as file 8
|
|
C and may thus be used as input to another run of the
|
|
C program
|
|
C INTRPL - switch indicating whether (and, if so, how) interpolate
|
|
C the initial model if the depth scales for the input model
|
|
C and the present depth scale are different
|
|
C = 0 - no interpolation, i.e. scale DEPTH coincides with DM
|
|
C > 0 - polynomial interpolation of the (INTRPL-1)th order
|
|
C
|
|
INCLUDE 'PARAMS.FOR'
|
|
INCLUDE 'MODELP.FOR'
|
|
PARAMETER (MINPUT=MLEVEL+4)
|
|
DIMENSION ESEMAT(MLEVEL,MLEVEL),BESE(MLEVEL),POPLTE(MLEVEL),
|
|
* TOTN(MDEPTH),PLTE(MLEVEL,MDEPTH)
|
|
COMMON ESEMAT,BESE,POPLTE,POPUL0(MLEVEL,MDEPTH),X(MINPUT),
|
|
* TEMP0(MDEPTH),ELEC0(MDEPTH),DENS0(MDEPTH),PPL0(MDEPTH),
|
|
* PPL(MDEPTH),DEPTH(MDEPTH),DM0(MDEPTH),DP(MDEPTH)
|
|
COMMON/NLTPOP/PNLT(MATOM,MION,MDEPTH)
|
|
common/quasex/iexpl(mlevel),iltot(mlevel)
|
|
C
|
|
NUMLT=3
|
|
IF(INMOD.EQ.2) NUMLT=4
|
|
READ(8,*) NDPTH,NUMPAR
|
|
READ(8,*) (DEPTH(I),I=1,NDPTH)
|
|
ND=NDPTH
|
|
NUMP=ABS(NUMPAR)
|
|
DO 30 ID=1,NDPTH
|
|
READ(8,*) (X(I),I=1,NUMP)
|
|
TEMP(ID)=X(1)
|
|
ELEC(ID)=X(2)
|
|
DENS(ID)=X(3)
|
|
TOTN(ID)=DENS(ID)/WMM(ID)+ELEC(ID)
|
|
CALL WNSTOR(ID)
|
|
CALL SABOLF(ID)
|
|
IP=NUMLT
|
|
IF(NUMPAR.LT.0) THEN
|
|
IP=IP+1
|
|
TOTN(ID)=X(IP)
|
|
END IF
|
|
IF(INMOD.EQ.2) IP=IP+1
|
|
c
|
|
c first compute LTE level populations for all levels,
|
|
c i.e. explicit, semi-explisit, and quasi-explicit
|
|
c
|
|
NLEV0=NLEVEL
|
|
TEMP(ID)=X(1)
|
|
ELEC(ID)=X(2)
|
|
DENS(ID)=X(3)
|
|
t=temp(id)
|
|
if(ifmol.gt.0.and.t.lt.tmolim) then
|
|
ipri=1
|
|
aein=elec(id)
|
|
an=totn(id)
|
|
call moleq(id,t,an,aein,ane,ipri)
|
|
else
|
|
if(imode.gt.-2) then
|
|
DO IAT=1,NATOM
|
|
ATTOT(IAT,ID)=DENS(ID)/WMM(ID)/YTOT(ID)*ABUND(IAT,ID)
|
|
END DO
|
|
else
|
|
DO IAT=1,NATOM
|
|
ATTOT(IAT,ID)=DENS(ID)/WMM(1)/YTOT(1)*ABUND(IAT,1)
|
|
END DO
|
|
end if
|
|
end if
|
|
CALL WNSTOR(ID)
|
|
CALL SABOLF(ID)
|
|
CALL RATMAT(ID,ESEMAT,BESE)
|
|
CALL LEVSOL(ESEMAT,BESE,POPLTE,NLEV0)
|
|
DO I=1,NLEV0
|
|
POPUL(I,ID)=POPLTE(I)
|
|
PLTE(I,ID)=POPLTE(I)
|
|
c if(id.eq.1) write(6,651) i,ip,popul(i,id),plte(i,id)
|
|
END DO
|
|
c
|
|
c if the input file fort.8 contains also NLTE level populations
|
|
c of b-factors, replace the LTE populations by those
|
|
c
|
|
IF(NUMP.GT.IP) THEN
|
|
NLEV0=NUMP-IP
|
|
DO I=1,NLEV0
|
|
j=iltot(i)
|
|
POPUL(J,ID)=X(IP+I)*RELAB(IATM(I),ID)
|
|
c if(id.eq.1) write(6,651) i,j,x(ip+i),popul(i,id)
|
|
c 651 format('in',2i4,1p2e12.4)
|
|
END DO
|
|
c DO I=1,NLEV0
|
|
c j=iltot(i)
|
|
c if(popul(j,id).le.0.) then
|
|
c IE=IEL(I)
|
|
c N0I=NFIRST(IE)
|
|
c NKI=NNEXT(IE)
|
|
c POPUL(J,ID)=ELEC(ID)*POPUL(iltot(NKI),ID)*SBF(I)
|
|
c end if
|
|
c END DO
|
|
c
|
|
c in the case the input "NLTE populations are in fact b-factors,
|
|
c compute the real populations
|
|
c
|
|
if(ibfac.eq.1) then
|
|
do i=1,nlev0
|
|
j=iltot(i)
|
|
popul(j,id)=popul(j,id)*plte(j,id)
|
|
end do
|
|
end if
|
|
END IF
|
|
30 CONTINUE
|
|
C
|
|
close(8)
|
|
c
|
|
write(6,600)
|
|
600 format(/' INPUT TLUSTY MODEL'/
|
|
* ' ------------------'/
|
|
* 1H ,8X,'MASS',9X,'T',9X,'NE',9X,'DENS'//)
|
|
nd=ndpth
|
|
DO 40 ID=1,ND
|
|
DM(ID)=DEPTH(ID)
|
|
write(6,601) id,dm(id),temp(id),elec(id),dens(id),
|
|
* popul(1,id)
|
|
601 format(i6,1pe10.3,0pf10.1,1p4e12.3)
|
|
40 CONTINUE
|
|
C
|
|
DO 100 ID=1,ND
|
|
BCON=ELEC(ID)/TEMP(ID)/SQRT(TEMP(ID))*2.0706E-16
|
|
DO 100 IONE=1,NION
|
|
ION=IZ(IONE)
|
|
IAT=NUMAT(IATM(NFIRST(IONE)))
|
|
NKI=NNEXT(IONE)
|
|
IF(ION.GT.0) PNLT(IAT,ION,ID)=POPUL(NKI,ID)/G(NKI)*BCON
|
|
100 CONTINUE
|
|
c
|
|
c check abundances
|
|
c
|
|
c CALL CHCKAB
|
|
RETURN
|
|
END
|