342 lines
10 KiB
Fortran
342 lines
10 KiB
Fortran
SUBROUTINE INILAM
|
|
C =================
|
|
C
|
|
C Auxiliary procedure for RESOLV
|
|
C initialization of model parameters for further use in RESOLV
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ITERAT.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
DIMENSION COL(MTRANS),DENS0(MDEPTH),SBW(MLEVEL),XE(MDEPTH),
|
|
* CLOC(MTRANS),ANTC(MDEPTH)
|
|
DIMENSION GRD(MDEPTH),pra(mdepth)
|
|
C
|
|
IF(INIT.EQ.1) THEN
|
|
C
|
|
C Before first iteration of complete linearization, ie. just after
|
|
C initialization of basic parameters and reading the initial model
|
|
C atmosphere - i.e. if INIT=1:
|
|
C only evaluation of collisional rates for all transitions, at all
|
|
C depths
|
|
C
|
|
anerel=0.5
|
|
if(teff.lt.8000.) anerel=0.01
|
|
CALL TDPINI
|
|
IF(IDISK.EQ.1) THEN
|
|
AMUV0=DMVISC**(ZETA0+UN)
|
|
AMUV1=UN-AMUV0
|
|
DMTOT=DM(ND)
|
|
EDISC=SIG4P*TEFF**4/DMTOT
|
|
END IF
|
|
DO ID=1,ND
|
|
CALL WNSTOR(ID)
|
|
FCOOL(ID)=0.
|
|
FPRD(ID)=0.
|
|
call sabolf(id)
|
|
c
|
|
c intialize b-factors
|
|
c
|
|
DO I=1,NLEVEL
|
|
BFAC(I,ID)=UN
|
|
SBW(I)=ELEC(ID)*SBF(I)*WOP(I,ID)
|
|
END DO
|
|
IF(.NOT.LTE.AND.IPSLTE.EQ.0 .and. id.lt.idlte) THEN
|
|
DO ION=1,NION
|
|
DO I=NFIRST(ION),NLAST(ION)
|
|
IF(POPUL(NNEXT(ION),ID).GT.0..and.iltlev(i).eq.0)
|
|
* BFAC(I,ID)=POPUL(I,ID)/
|
|
* (POPUL(NNEXT(ION),ID)*SBW(I))
|
|
END DO
|
|
END DO
|
|
END IF
|
|
c
|
|
c fully fixed part of the total charge
|
|
c
|
|
qfix(id)=0.
|
|
do i=1,nlevel
|
|
if(imodl(i).lt.0.or.iifix(iatm(i)).gt.0) then
|
|
ch=iz(iel(i))-1
|
|
il=ilk(i)
|
|
if(il.gt.0) ch=iz(il)+(iz(il)-1)*usum(il)*elec(id)
|
|
qfix(id)=qfix(id)+ch*popul(i,id)
|
|
end if
|
|
end do
|
|
c
|
|
c collisional rates
|
|
c
|
|
IF(.NOT.LTE) THEN
|
|
CALL COLIS(ID,TEMP(ID),COL,CLOC)
|
|
DO I=1,NTRANS
|
|
COLRAT(I,ID)=COL(I)
|
|
COLTAR(I,ID)=CLOC(I)
|
|
END DO
|
|
END IF
|
|
END DO
|
|
CALL ODFMER
|
|
C
|
|
c initialize viscosity parameters for disks
|
|
c
|
|
IF(IDISK.EQ.1) CALL VISINI
|
|
c
|
|
c dielectronic recombination
|
|
c
|
|
if(ifdiel.gt.0) call dietot
|
|
C
|
|
C in case of ISPLIN>4 - (i.e. the DFE solver for transfer equation)
|
|
C evaluate Planck function for the first estimate of RAD for
|
|
C the electron scattering source function
|
|
C
|
|
IF(ISPLIN.GE.5) THEN
|
|
DO ID=1,ND
|
|
DO IJ=1,NFREQ
|
|
RAD(IJ,ID)=BNUE(IJ)/(EXP(HKT1(ID)*FREQ(IJ))-UN)
|
|
END DO
|
|
END DO
|
|
END IF
|
|
RETURN
|
|
END IF
|
|
C
|
|
C ==================================================================
|
|
C Immediately after a completed iteration of complete linearization:
|
|
C ==================================================================
|
|
C
|
|
PRAD=0.
|
|
DO ID=1,ND
|
|
C
|
|
C save some old quantities
|
|
C
|
|
AOLD=DENS(ID)/WMM(ID)+ELEC(ID)
|
|
XE(ID)=UN-ELEC(ID)/AOLD
|
|
C
|
|
C Depth by depth:
|
|
C set up individual arrays of model parameters, namely:
|
|
C temperature - TEMP (if temperature is one of variables,
|
|
C ie. if radiative equilibrium equation is solved)
|
|
C
|
|
IF(INRE.NE.0) TEMP(ID)=PSY0(NFREQE+INRE,ID)
|
|
C
|
|
C electron density - ELEC (if ELEC is one of variables, ie. if
|
|
C particle conservation equation is solved)
|
|
C
|
|
IF(INPC.NE.0) ELEC(ID)=PSY0(NFREQE+INPC,ID)
|
|
C
|
|
C density - DENS, calculated by means of total particle number
|
|
C density and electron density (if total particle
|
|
C number density is one of variables, ie. if hydrostatic
|
|
C equilibrium equation is solved - INHE.GT.0)
|
|
C
|
|
IF(INHE.NE.0) TOTN(ID)=PSY0(NFREQE+INHE,ID)
|
|
C
|
|
IF(IFIXDE.EQ.0) THEN
|
|
DENS0(ID)=DENS(ID)
|
|
IF(INHE.NE.0) THEN
|
|
DENS(ID)=WMM(ID)*(PSY0(NFREQE+INHE,ID)-ELEC(ID))
|
|
DPLP=DENS0(ID)*DPSILN
|
|
DPLM=DENS0(ID)/DPSILN
|
|
IF(DENS(ID).GT.DPLP) DENS(ID)=DPLP
|
|
IF(DENS(ID).LT.DPLM) DENS(ID)=DPLM
|
|
ELSE
|
|
if(ioptab.lt.-1) then
|
|
PGS(ID)=PTOTAL(ID)
|
|
DENS(ID)=RHOEOS(TEMP(ID),PGS(ID))
|
|
end if
|
|
END IF
|
|
END IF
|
|
C
|
|
C or again density, but calculated by means of the
|
|
C massive particle density,
|
|
C if this is one of variables - INMP.GT.0)
|
|
C
|
|
IF(INMP.NE.0) DENS(ID)=WMM(ID)*PSY0(NFREQE+INMP,ID)
|
|
C
|
|
C geometrical distance
|
|
C
|
|
IF(INZD.GT.0) ZD(ID)=PSY0(NFREQE+INZD,ID)
|
|
C
|
|
C the temperature logarithmic gradient delta
|
|
C
|
|
IF(INDL.NE.0) DELTA(ID)=PSY0(NFREQE+INDL,ID)
|
|
C
|
|
ANMA(ID)=DENS(ID)/WMM(ID)
|
|
ANTO(ID)=ANMA(ID)+ELEC(ID)
|
|
END DO
|
|
C
|
|
c avoid oscillations in temperature
|
|
c
|
|
IF(IOSCOR.NE.0) CALL OSCCOR
|
|
C
|
|
C if needed, formal solution of the hydrostatic equilibrium to get
|
|
C new total particle density, and possibly electron density
|
|
C
|
|
IF(IHECOR.GE.2) THEN
|
|
ID=1
|
|
PTUR=HALF*VTURB(ID)*VTURB(ID)*DENS(ID)
|
|
ANTC(ID)=(DM(ID)*GRAV-PRD0-PTUR)/BOLK/TEMP(ID)
|
|
IF(ANTC(ID).LE.0) ANTC(ID)=DENS(ID)/WMM(ID)+ELEC(ID)
|
|
DO ID=2,ND
|
|
PTUR=HALF*VTURB(ID)*VTURB(ID)*DENS(ID)
|
|
PTURM=HALF*VTURB(ID-1)*VTURB(ID-1)*DENS(ID-1)
|
|
ANTC(ID)=(GRAV*(DM(ID)-DM(ID-1))+
|
|
* BOLK*TEMP(ID-1)*ANTC(ID-1)-
|
|
* PRADT(ID)+PRADT(ID-1)-PTUR+PTURM)/
|
|
* BOLK/TEMP(ID)
|
|
END DO
|
|
DO ID=1,ND
|
|
ELEC(ID)=(UN-XE(ID))*ANTC(ID)
|
|
DENS(ID)=WMM(ID)*(ANTC(ID)-ELEC(ID))
|
|
ANMA(ID)=DENS(ID)/WMM(ID)
|
|
ANTO(ID)=ANMA(ID)+ELEC(ID)
|
|
END DO
|
|
END IF
|
|
C
|
|
C other depth-dependent quantities
|
|
C
|
|
do id=1,nd
|
|
PGS(ID)=(DENS(ID)/WMM(ID)+ELEC(ID))*BOLK*TEMP(ID)
|
|
end do
|
|
c
|
|
IF(IOPTAB.GE.0) THEN
|
|
DO ID=1,ND
|
|
C
|
|
C occupation probabilities
|
|
C
|
|
CALL WNSTOR(ID)
|
|
C
|
|
C Evaluation of collision rates for the new temperature
|
|
C
|
|
IF(.NOT.LTE) THEN
|
|
call sabolf(id)
|
|
CALL COLIS(ID,TEMP(ID),COL,CLOC)
|
|
DO I=1,NTRANS
|
|
COLRAT(I,ID)=COL(I)
|
|
COLTAR(I,ID)=CLOC(I)
|
|
END DO
|
|
END IF
|
|
END DO
|
|
END IF
|
|
C
|
|
CALL CONCOR
|
|
CALL TDPINI
|
|
C------------------------
|
|
C
|
|
C new populations
|
|
C
|
|
C 1. first possibility (useful for pure CL) --
|
|
C new populations, ie. those calculated directly
|
|
C by the complete linearization (as old population + corresponding
|
|
C corrections) are not used.
|
|
C Instead, new populations are calculated from the new radiation
|
|
C field basically by lambda iteration, ie.
|
|
C a) evaluation of new radiative rates for explicit transitions
|
|
C (radiative rates in fixed-option transitions are not updated
|
|
C at this step)
|
|
C
|
|
if(ifpopr.le.0.or.lte.OR.IFRYB.GT.0) then
|
|
IMOR=0
|
|
IF(IFRYB.GT.2) IMOR=1
|
|
IF(.NOT.LTE) CALL RATES1(IMOR)
|
|
C
|
|
C b) depth by depth evaluation of new populations by solving
|
|
C statistical equilibrium equation with previously evaluated
|
|
C rates, together with constraints - precisely as in RESOLV
|
|
C
|
|
DO ID=1,ND
|
|
CALL STEQEQ(ID,POP,1)
|
|
IF(.NOT.LCHC.and.iter.lt.ielcor) CALL ELCOR(ID)
|
|
END DO
|
|
if(iprind.eq.2) call output
|
|
c
|
|
C 2. second possibility -- populations from linearization corrections
|
|
C
|
|
ELSE
|
|
DO ID=1,ND
|
|
DO I=1,NLEVEL
|
|
if(iifix(iatm(i)).ne.1) then
|
|
II=IIEXP(I)
|
|
IF(II.GT.0) THEN
|
|
III=IINONZ(II)
|
|
IF(III.GT.0) THEN
|
|
POPUL(I,ID)=PSY0(NFREQE+INSE+III-1,ID)
|
|
ELSE
|
|
POPUL(I,ID)=0.
|
|
END IF
|
|
ELSE IF(II.LT.0) THEN
|
|
III=IINONZ(-II)
|
|
IF(III.GT.0) THEN
|
|
POPUL(I,ID)=PSY0(NFREQE+INSE+III-1,ID)*
|
|
* SBPSI(I,ID)
|
|
ELSE
|
|
POPUL(I,ID)=0.
|
|
END IF
|
|
ELSE
|
|
III=IINONZ(IIEXP(ILTREF(I,ID)))
|
|
IF(III.GT.0) THEN
|
|
POPUL(I,ID)=PSY0(NFREQE+INSE+III-1,ID)*
|
|
* SBPSI(I,ID)
|
|
ELSE
|
|
POPUL(I,ID)=0.
|
|
END IF
|
|
END IF
|
|
END IF
|
|
END DO
|
|
END DO
|
|
END IF
|
|
CALL ODFMER
|
|
c
|
|
c dielectronic recombination
|
|
c
|
|
if(ifdiel.gt.0) call dietot
|
|
C
|
|
if(ifryb.gt.0) call rybheq
|
|
C
|
|
C solution of the transfer equation with Compton scattering
|
|
C
|
|
if(icompt.gt.0) then
|
|
call opaini(1)
|
|
call comset
|
|
call rtecom
|
|
end if
|
|
c
|
|
IF(IDISK.EQ.1) call visini
|
|
c
|
|
c if(ifryb.eq.0) then
|
|
CALL OPAINI(1)
|
|
DO ID=1,ND
|
|
GRD(ID)=0.
|
|
pra(id)=0.
|
|
pradt(id)=0.
|
|
END DO
|
|
prd0=0.
|
|
DO IJ=1,NFREQ
|
|
CALL OPACF1(IJ)
|
|
CALL RTEFR1(IJ)
|
|
GRD(1)=GRD(1)+W(IJ)*ABSO1(1)*FH(IJ)*RAD1(1)
|
|
DO ID=2,ND
|
|
GRD(ID)=GRD(ID)+(RAD1(ID)*FAK1(ID)-RAD1(ID-1)*FAK1(ID-1))*
|
|
* W(IJ)
|
|
pra(id)=pra(id)+RAD1(ID)*FAK1(ID)*W(IJ)
|
|
END DO
|
|
pra(1)=pra(1)+RAD1(1)*FAK1(1)*W(IJ)-ABSO1(1)*W(IJ)*
|
|
* (RAD1(1)*FH(IJ)-HEXTRD(IJ))
|
|
END DO
|
|
GRD(1)=PCK*GRD(1)/DENS(1)
|
|
c
|
|
do id=1,nd
|
|
pra(id)=pra(id)*pck
|
|
pradt(id)=pradt(id)*pck
|
|
c pradfc(id)=pra(id)/(2.5213e-15*temp(id)**4)
|
|
end do
|
|
C
|
|
if(idisk.eq.0) then
|
|
PGS(1)=DM(1)*(GRAV-GRD(1))
|
|
DO ID=2,ND
|
|
PGS(ID)=PGS(ID-1)-PCK*GRD(ID)+GRAV*(DM(ID)-DM(ID-1))
|
|
END DO
|
|
end if
|
|
c
|
|
RETURN
|
|
END
|