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

287 lines
10 KiB
Fortran

SUBROUTINE CHANGE
C =================
C
C This procedure controls an evaluation of initial level
C populations in case where the system of explicit levels
C (ie. the choice of explicit level, their numbering, or their
C total number) is not consistent with that for the input level
C populations read by procedure INPMOD.
C Obviously, this procedure need be used only for NLTE input models.
C
C ICHANG < 0 - general change of populations as described below
C > 0 - a simplified change; original data for the input
C model are required to assign the input NLTE populations
C to the levels in the new models; all additional
C levels are assumed having LTE populations.
C ICHANG is the unit number for the data file of old model.
C
C Case ICHANG < 0:
C
C Input from unit 5:
C For each explicit level, II=1,NLEVEL, the following parameters:
C IOLD - NE.0 - means that population of this level is
C contained in the set of input populations;
C IOLD is then its index in the "old" (i.e. input)
C numbering.
C All the subsequent parameters have no meaning
C in this case.
C - EQ.0 - means that this level has no equivalent in the
C set of "old" levels. Population of this level
C has thus to be evaluated.
C MODE - indicates how the population is evaluated:
C = 0 - population is equal to the population of the "old"
C level with index ISIOLD, multiplied by REL;
C = 1 - population assumed to be LTE, with respect to the
C first state of the next ionization degree whose
C population must be contained in the set of "old"
C (ie. input) populations, with index NXTOLD in the
C "old" numbering.
C The population determined of this way may further
C be multiplied by REL.
C = 2 - population determined assuming that the b-factor
C (defined as the ratio between the NLTE and
C LTE population) is the same as the b-factor of
C the level ISINEW (in the present numbering). The
C level ISINEW must have the equivalent in the "old"
C set; its index in the "old" set is ISIOLD, and the
C index of the first state of the next ionization
C degree, in the "old" numbering, is NXTSIO.
C The population determined of this way may further
C be multiplied by REL.
C = 3 - level corresponds to an ion or atom which was not
C explicit in the old system; population is assumed
C to be LTE.
C NXTOLD - see above
C ISINEW - see above
C ISIOLD - see above
C NXTSIO - see above
C REL - population multiplier - see above
C if REL=0, the program sets up REL=1
C
INCLUDE 'IMPLIC.FOR'
INCLUDE 'BASICS.FOR'
INCLUDE 'ATOMIC.FOR'
INCLUDE 'MODELQ.FOR'
character*20 fnstd
dimension n0old(30,30),n1old(30,30)
dimension katold(2,30),vtbold(mdepth)
COMMON POPUL0(MLEVEL,MDEPTH),POPULL(MLEVEL,MDEPTH),
* ESEMAT(MLEVEL,MLEVEL),BESE(MLEVEL),POPL(MLEVEL)
C
PARAMETER (S = 2.0706D-16)
IF(ICHANG.LT.0) THEN
IFESE=0
DO 100 II=1,NLEVEL
READ(IBUFF,*) IOLD,MODE,NXTOLD,ISINEW,ISIOLD,NXTSIO,REL
IF(REL.EQ.0.) REL=1.
IF(MODE.GE.3) IFESE=IFESE+1
DO 90 ID=1,ND
IF(IOLD.EQ.0) GO TO 10
POPUL0(II,ID)=POPUL(IOLD,ID)
GO TO 90
10 IF(MODE.NE.0) GO TO 20
POPUL0(II,ID)=POPUL(ISIOLD,ID)*REL
GO TO 90
20 T=TEMP(ID)
ANE=ELEC(ID)
IF(MODE.GE.3) GO TO 40
NXTNEW=NNEXT(IEL(II))
SB=S/T/SQRT(T)*G(II)/G(NXTNEW)*EXP(ENION(II)/T/BOLK)
IF(MODE.GT.1) GO TO 30
POPUL0(II,ID)=SB*ANE*POPUL(NXTOLD,ID)*REL
GO TO 90
30 KK=ISINEW
KNEXT=NNEXT(IEL(KK))
SBK=S/T/SQRT(T)*G(KK)/G(KNEXT)*EXP(ENION(KK)/T/BOLK)
POPUL0(II,ID)=SB/SBK*POPUL(NXTOLD,ID)/POPUL(NXTSIO,ID)*
* POPUL(ISIOLD,ID)*REL
GO TO 90
40 IF(IFESE.EQ.1) THEN
LTE0=LTE
LTE=.TRUE.
do iii=1,nlevel
if(wop(iii,id).eq.0.) wop(iii,id)=1.
end do
CALL STEQEQ(ID,POPL,0)
DO III=1,NLEVEL
POPULL(III,ID)=POPL(III)
END DO
LTE=LTE0
END IF
POPUL0(II,ID)=POPULL(II,ID)
90 CONTINUE
100 CONTINUE
DO I=1,NLEVEL
DO ID=1,ND
POPUL(I,ID)=POPUL0(I,ID)
END DO
END DO
C
C simplified change - no additional input (the case ICHANG > 0)
C
ELSE
LTE0=LTE
LTE=.TRUE.
DO ID=1,ND
do ii=1,nlevel
if(wop(ii,id).eq.0.) wop(ii,id)=1.
end do
CALL STEQEQ(ID,POPL,0)
DO II=1,NLEVEL
POPUL0(II,ID)=POPL(II)
END DO
END DO
C
IF(ICHANG.EQ.1) THEN
DO II=NLEV0+1,NLEVEL
DO ID=1,ND
POPUL(II,ID)=POPUL0(II,ID)
END DO
END DO
C
ELSE
modr=0
rewind 1
read(1,*,err=200,end=200) modr
200 continue
call readbf(ichang)
if(modr.eq.0) then
read(95,*) tfold,grold
read(95,*) ltd1,ltd2
read(95,*) fnstd
read(95,*) nfrd
if(nfrd.lt.0) then
nfrd=-nfrd
do ij=1,nfrd
read(95,*) frold
end do
endif
read(95,*) natold
if(natold.lt.0) natold=-natold
do ia=1,natold
read(95,*) iao,abnold
if(abnold.gt.1.e6) read(95,*) (vtbold(i),i=1,ndold)
end do
nlold=0
read(95,*) iato,izo,nlvo,ilasti,ilvi,instd
if(instd.ne.0) read(95,*) idui
do while (ilasti.ge.0)
n0old(iato,izo+1)=nlold+1
n1old(iato,izo+1)=nlold+nlvo
nlold=nlold+nlvo
read(95,*) iato,izo,nlvo,ilasti,ilvi,instd
if(instd.ne.0) read(95,*) idui
end do
else
read(95,*) tfold,grold,hmold
read(95,*) ltd1,ltd2,lcold,ispold,chmold
if(ispold.lt.0) read(95,*,err=203) iol1,iol2,iol3,iol4,
. iol5,iol6,iol7
if(iol6.ge.2) read(95,*) djmold
203 read(95,*,err=204) nitold,ndold,natold,niold,nlvold,
. iol1,iol2,iarold,iol4
204 continue
if(iol1.gt.10) then
read(95,*,err=205) iol1,iol2,iol3,iol4,iol5
read(95,*,err=205) iol1,iol2
read(95,*,err=205) iol1,iol2
end if
205 continue
if(niold.lt.0) then
niold=-niold
read(95,*,err=206) iol1,iol2,iol3
end if
206 continue
if(iarold.le.-100 .and. iarold.gt.-200) then
iarold=-iarold-100
read(95,*) iol1
endif
read(95,*) nfrd
if(nfrd.gt.0) then
nfrd=-nfrd
do ij=1,nfrd
read(95,*) frold
end do
else
nfrd=-nfrd
read(95,*) frold
end if
read(95,*,err=211) iol1,iol2,iol3
if(iol3.lt.0) read(95,*) pzold
211 continue
read(95,*) iol1,vtbol
if(iol1.ne.0) read(95,*) (vtbold(i),i=1,ndold)
read(95,*) natsold
if(natsold.lt.0) natsold=-natsold
iat=0
do ia=1,natsold
read(95,*) iol1,iol2,iol3,iol4,iol5,abnold
if(abnold.gt.1.e6) read(95,*) (vtbold(i),i=1,ndold)
if(iol1.eq.2) then
iat=iat+1
katold(1,iat)=iol2
katold(2,iat)=iol3
end if
end do
do ii=1,niold
read(95,*) k0old,k1old,k2old,izo
if(k0old.lt.0) then
k0old=-k0old
read(95,*) iol1
end if
do ia=1,iat
if(k0old.ge.katold(1,ia) .and. k1old.ge.katold(2,ia))
. iaol=ia
end do
n0old(iaol,izo)=k0old
n1old(iaol,izo)=k1old
n0old(iaol,izo+1)=k2old
n1old(iaol,izo+1)=k2old
end do
end if
C
WRITE(6,600)
600 FORMAT(' Levels: OLD model -> NEW model',/
. ' ------------------------------')
DO 300 II=1,NION
N0NEW=NFIRST(II)
N1NEW=NLAST(II)
IANEW=NUMAT(IATM(N0NEW))
IZNEW=IZ(IEL(N0NEW))
IF(N1OLD(IANEW,IZNEW).EQ.0) GO TO 300
KOLD=N1OLD(IANEW,IZNEW)-N0OLD(IANEW,IZNEW)
KNEW=NLAST(II)-NFIRST(II)
IF(KOLD.LT.KNEW) KNEW=KOLD
JL=N0OLD(IANEW,IZNEW)-1
DO IL=NFIRST(II),NFIRST(II)+KNEW
JL=JL+1
WRITE(6,601) JL,IL
601 FORMAT(10X,I8,5X,I8)
DO ID=1,ND
POPUL0(IL,ID)=POPUL(JL,ID)
END DO
END DO
300 CONTINUE
DO 310 II=1,NATOM
N0NEW=NKA(II)
IANEW=NUMAT(IATM(N0NEW))
IZNEW=IZ(IEL(N0NEW))+1
IF(N0OLD(IANEW,IZNEW).EQ.0) GO TO 310
WRITE(6,601) N0OLD(IANEW,IZNEW),N0NEW
DO ID=1,ND
POPUL0(N0NEW,ID)=POPUL(N0OLD(IANEW,IZNEW),ID)
END DO
310 CONTINUE
C
DO II=1,NLEVEL
DO ID=1,ND
POPUL(II,ID)=POPUL0(II,ID)
END DO
END DO
END IF
LTE=LTE0
C
END IF
RETURN
END