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

347 lines
9.0 KiB
Fortran

SUBROUTINE INMOLI(ILIST)
C ========================
C
C read in the input molecular line list,
C selection of lines that may contribute,
C set up auxiliary fields containing line parameters,
C
C Input of line data - unit 20:
C
C For each line, one (or two) records, containing:
C
C ALAM - wavelength (in nm)
C ANUM - code of the modelcule (as in Kurucz)
C (eg. 101.00 = H2; 607.00 = CN)
C GF - log gf
C EXCL - excitation potential of the lower level (in cm*-1)
C GR - gamma(rad)
C GS - gamma(Stark)
C GW - gamma(VdW)
C
INCLUDE 'PARAMS.FOR'
INCLUDE 'MODELP.FOR'
INCLUDE 'SYNTHP.FOR'
INCLUDE 'LINDAT.FOR'
COMMON/LIMPAR/ALAM0,ALAM1,FRMIN,FRLAST,FRLI0,FRLIM
COMMON/BLAPAR/RELOP,SPACE0,CUTOF0,TSTD,DSTD,ALAMC
COMMON/NXTINM/ALMM00,ALSM00
common/alendm/alend(mmlist)
common/brdstd/gsstd,gwstd
character*80 dum
dimension x(9)
PARAMETER (PI4=7.95774715E-2)
PARAMETER (C1 = 2.3025851,
* C2 = 4.2014672,
* C3 = 1.4387886,
* CNM = 2.997925D17,
* EXT0 = 3.17,
* UN = 1.0,
* TEN = 10.,
* HUND = 1.D2,
* TENM4 = 1.D-4,
* TENM8 = 1.D-8,
* OP4 = 0.4,
* AGR0=2.4734E-22,
* XEH=13.595, XET=8067.6, XNF=25.,
* R02=2.5, R12=45., VW0=4.5E-9)
C
c DATA INLSET /0/
C
if(imode.ne.-3.and.temp(idstd).gt.tmolim) return
IUNIT=IUNITM(ILIST)
if(ibin(ilist).eq.0) then
open(unit=iunit,file=amlist(ilist),status='old')
else
open(unit=iunit,file=amlist(ilist),form='unformatted',
* status='old')
end if
C
c define a conversion table between Kurucz notation and Tsuji table
c through array MOLIND
C
do i=1,11000
molind(i)=0
end do
molind(101)=2
molind(106)=5
molind(107)=12
molind(108)=4
molind(111)=122
molind(112)=32
molind(114)=17
molind(116)=16
molind(120)=34
molind(124)=198
molind(126)=214
molind(606)=8
molind(607)=7
molind(608)=6
molind(614)=21
molind(616)=20
molind(707)=9
molind(708)=11
molind(714)=24
molind(716)=23
molind(808)=10
molind(812)=126
molind(813)=134
molind(814)=25
molind(816)=26
molind(820)=179
molind(822)=29
molind(823)=30
molind(10108)=3
c
c iunit=19+ilist
C
c ================================
c detect the type of the line list
c
ivdwli(ilist)=0
ibroad=1
c
c text list
c
if(ibin(ilist).eq.0) then
read(iunit,'(a80)') dum
read(dum,*,iostat=kst1) (x(i),i=1,9)
np=9
if(kst1.ne.0) then
read(dum,*,iostat=kst2) (x(i),i=1,7)
np=7
if(kst2.ne.0) then
read(dum,*,iostat=kst3) (x(i),i=1,4)
ibroad=0
np=4
if(kst3.ne.0) then
write(*,*) 'no applicable format of line list',ilist
end if
end if
end if
if(np.eq.9) ivdwli(ilist)=1
else
c
c binary list
c
read(iunit,err=110) (x(i),i=1,9)
np=9
go to 150
110 continue
read(iunit,err=120) (x(i),i=1,7)
np=7
go to 150
120 continue
read(iunit,err=130) (x(i),i=1,4)
ibroad=0
np=4
go to 150
130 continue
150 continue
if(np.eq.9) ivdwli(ilist)=1
if(np.eq.9) ivdwli(ilist)=1
end if
c =========================
c
ALAST=CNM/FRLAST
ALASTM(ILIST)=ALAST
IL=0
IF(NXTSEM(ILIST).EQ.1) THEN
ALAM0=ALM00
ALASTM(ILIST)=ALST00
FRLASM(ILIST)=CNM/ALASTM(ILIST)
NXTSEM(ILIST)=0
REWIND IUNIT
END IF
ALMM00=ALAM0
c ALASTM(ILIST)=CNM/FRLAST
c FRLASM(ILIST)=CNM/ALASTM(ILIST)
DOPSTD=1.E7/ALAM0*DSTD
DOPLAM=ALAM0*ALAM0/CNM*DOPSTD
AVAB=ABSTD(IDSTD)*RELOP
ASTD=1.0
c IF(GRAV.GT.6.) ASTD=0.1
CUTOFF=CUTOF0
ALAST=CNM/FRLAST
C
C first part of reading line list - read only lambda, and
C skip all lines with wavelength below ALAM0-CUTOFF
C
REWIND IUNIT
ALAM=0.
IJC=2
c
7 if(ibin(ilist).eq.0) then
READ(IUNIT,510) ALAM
else
read(iunit) alam
end if
510 FORMAT(F10.4)
IF(ALAM.LT.ALAM0-CUTOFF) GO TO 7
BACKSPACE(IUNIT)
GO TO 10
c
c read the line list
c
ill=0
8 continue
10 continue
ill=ill+1
c ivdwli(ilist)=1
if(ibin(ilist).eq.0) then
c if(ivdwli(ilist).ne.0) then
if(np.eq.9) then
read(iunit,*,end=100) alam,anum,gf,excl,gr,gh2,xnh2,ghe,xnhe
else if(np.eq.7) then
READ(IUNIT,*,END=100,err=8) ALAM,ANUM,GF,EXCL,GR,GS,GW
else
read(iunit,*,end=100,err=8) alam,anum,gf,excl
gr=2.4e13/alam**2
gs=gsstd
gw=gwstd
end if
else
c if(ivdwli(ilist).ne.0) then
if(np.eq.9) then
read(iunit,end=100) alam,anum,gf,excl,gr,gh2,xnh2,ghe,xnhe
else if(np.eq.7) then
READ(IUNIT,END=100) ALAM,ANUM,GF,EXCL,GR,GS,GW
else
read(iunit,end=100) alam,anum,gf,excl
gr=2.4e13/alam**2
gs=gsstd
gw=gwstd
end if
end if
C
c change wavelength to vacuum for lambda > 2000
c
if(alam.gt.200..and.vaclim.gt.2000.) then
wl0=alam*10.
ALM=1.E8/(WL0*WL0)
XN1=64.328+29498.1/(146.-ALM)+255.4/(41.-ALM)
WL0=WL0*(XN1*1.D-6+UN)
alam=wl0*0.1
END IF
C
C first selection : for a given interval
C
IF(ALAM.GT.ALASTM(ILIST)+CUTOFF) GO TO 100
C
C second selection : for line strengths
C
FR0=CNM/ALAM
icod=int(anum+tenm4)
c IF(ICOD.EQ.823) go to 10
imol=molind(icod)
if(imol.le.0.or.imol.gt.nmolec) go to 10
EXCL=ABS(EXCL)
GFP=C1*GF-C2
EPP=C3*EXCL
gx=gfp-epp/tstd
ab0=0.
c
if(ndstep.eq.0.and.ifwin.eq.0) then
c
c old procedure for line rejection
c
if(gx.gt.-30)
* AB0=EXP(GFP-EPP/TSTD)*RRMOL(IMOL,IDSTD)/DOPSTD/AVAB
IF(AB0.LT.UN) GO TO 10
else
c
c new procedure for line rejection
c
do ijcn=ijc,nfreqc
if(fr0.ge.freqc(ijcn)) go to 12
end do
12 continue
ijc=ijcn
if(ijc.gt.nfreqc) ijc=nfreqc
c
tkm=1.65e8/ammol(imol)
DP0=3.33564E-11*FR0
do id=1,nd,ndstep
td=temp(id)
gx=gfp-epp/td
ab0=0.
if(gx.gt.-30) then
dops=dp0*sqrt(tkm*td+vturb(id))
AB0=EXP(gx)*RRMOL(IMOL,ID)/(DOPS*abstdw(ijc,id)*relop)
end if
if(ab0.ge.un) go to 15
end do
GO TO 10
end if
c
C truncate line list if there are more lines than maximum allowable
C (given by MLIN0 - see include file LINDAT.FOR)
C
15 CONTINUE
IL=IL+1
IF(IL.GT.MLINM0) THEN
WRITE(6,601) ALAM
IL=MLINM0
ALASTM(ILIST)=CNM/FREQM(IL,ILIST)-CUTOFF
FRLASM(ILIST)=CNM/ALASTM(ILIST)
NXTSEM(ILIST)=1
GO TO 100
END IF
C
C =============================================
C line is selected, set up necessary parameters
C =============================================
C
C evaluation of EXTIN0 - the distance (in delta frequency) where
C the line is supposed to contribute to the total opacity
C
EX0=AB0*ASTD*10.
EXT=EXT0
IF(EX0.GT.TEN) EXT=SQRT(EX0)
EXTIN0=EXT*DOPSTD
C
C store parameters for selected lines
C
FREQM(IL,ILIST)=FR0
EXCLM(IL,ILIST)=real(EPP)
GFM(IL,ILIST)=real(GFP)
EXTINM(IL,ILIST)=real(EXTIN0)
INDATM(IL,ILIST)=imol
C
C ****** line broadening parameters *****
C assuming for Stark 1.e-8*effnsq**5/2, with effnsq=25
C
GRM(IL,ILIST)=real(GR*PI4)
GSM(IL,ILIST)=real(GS*PI4*3.125e-5)
GWM(IL,ILIST)=real(GW*PI4)
c IF(imol.eq.30) gwm(il,ilist)=0.
if(ivdwli(ilist).ne.0) then
gvdwh2(il,ilist)=real(gh2)
gexph2(il,ilist)=real(xnh2)
gvdwhe(il,ilist)=real(ghe)
gexphe(il,ilist)=real(xnhe)
gsm(il,ilist)=0.
gwm(il,ilist)=0.
end if
C
GO TO 10
100 NLINM0(ILIST)=IL
nlinmt(ilist)=nlinmt(ilist)+nlinm0(ilist)
alend(ilist)=cnm/fr0
C
xln=float(il)*1.e-6
WRITE(6,611) IUNIT,trim(amlist(ilist)),XLN
611 FORMAT(/' --------------------------------------------'/
*' MOLECULAR LINES - FROM UNIT ',i3,
*', FILE ',a,':',f8.3,' M'/
*' --------------------------------------------'/)
601 FORMAT('0 **** MORE LINES THAN MLINM0, LINE LIST TRUNCATED '/
*' AT LAMBDA',F15.4,' NM'/)
RETURN
END