SpectraRust/synspec/extracted/molini.f
2026-03-19 14:05:33 +08:00

79 lines
1.9 KiB
Fortran

subroutine molini
c =================
c
c Initialization of the molecular equilibrium
c
INCLUDE 'PARAMS.FOR'
INCLUDE 'MODELP.FOR'
common/moltst/pfmol(600,mdepth),anmol(600,mdepth),
* pfato(100,mdepth),anato(100,mdepth),
* pfion(100,mdepth),anion(100,mdepth)
dimension hpo(mdepth)
c
aeinit=1.0
c
do 10 id=1,nd
t=temp(id)
tln=log(t)*1.5
thl=11605./t
t32=exp(tln)
do i=1,MMOLEC
rrmol(i,id)=0.
end do
hpo(id)=DENS(ID)/WMM(ID)/YTOT(ID)
if(t.gt.tmolim) go to 10
HPOP=DENS(ID)/WMM(ID)/YTOT(ID)
an=dens(id)/wmm(id)+elec(id)
aeinit=0.1*an
if(t.lt.4000.) aeinit=0.01*an
call moleq(id,t,an,aeinit,ane,0)
c next initial guess will be the last ane determined for
c previous depth point
aeinit=ane
c
if (id.eq.idstd) then
write(6,600)
nmol=nmolec
if(id.eq.1) nmol=32
do i=1,nmol
write(6,601) i, cmol(i), rrmol(i,id), rrmol(i,id)/hpop
end do
end if
600 format(/ 'Molecular number densities at the standard depth'/)
601 format(i4,1x,A8,1x,1pe12.2,1x,e12.2)
10 continue
c update atomic populations once molecular densities are calculated
if(imode.lt.-4) then
do i=1,nlevel
iat=numat(iatm(i))
ion=iz(iel(i))
ii=nfirst(iel(i))
ener=(enion(ii)-enion(i))/bolk
if((enion(i).eq.0).and.(ilk(i).gt.0)) then
ener=0.
ion=ion+1
end if
if(ifwop(i).ge.0) then
do id=1,nd
popul(i,id)=rrr(id,ion,iat)*g(i)
* *exp(-ener/temp(id))
if(iat.eq.1.and.ion.eq.0) popul(i,id)=anhm(id)
end do
endif
end do
end if
c
return
end