203 lines
5.0 KiB
Fortran
203 lines
5.0 KiB
Fortran
SUBROUTINE RYBSOL
|
|
C =================
|
|
C
|
|
C driver of the complete-linearization solution in the Rybicki
|
|
C formalism
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
INCLUDE 'ARRAY1.FOR'
|
|
INCLUDE 'ITERAT.FOR'
|
|
COMMON/RYBMTX/RA(MDEPTH),RB(MDEPTH),RC(MDEPTH),VR(MDEPTH),
|
|
* UA(MDEPTH),UB(MDEPTH),UC(MDEPTH),
|
|
* VA(MDEPTH),VB(MDEPTH),VC(MDEPTH),WR(MDEPTH),
|
|
* WM(MDEPTH,MDEPTH)
|
|
common/imodlc/imodl0(mlevel)
|
|
DIMENSION AL(MDEPTH),AU(MDEPTH),VAL(MDEPTH),UCOL(MDEPTH),
|
|
* VAU(MDEPTH,MDEPTH),CHANGT(MDEPTH)
|
|
dimension pop1(mlevel),babs(mlevel,mdepth)
|
|
C
|
|
C zeroing needed quantities
|
|
C
|
|
DO ID=1,ND
|
|
REIT(ID)=0.
|
|
REIN(ID)=0.
|
|
AREIT(ID)=0.
|
|
AREIN(ID)=0.
|
|
CREIT(ID)=0.
|
|
CREIN(ID)=0.
|
|
REDT(ID)=0.
|
|
REDTM(ID)=0.
|
|
REDTP(ID)=0.
|
|
REDN(ID)=0.
|
|
REDNM(ID)=0.
|
|
REDNP(ID)=0.
|
|
FCOOLI(ID)=0.
|
|
FLFIX(ID)=0.
|
|
FLRD(ID)=0.
|
|
ABROSD(ID)=0.
|
|
SUMDPL(ID)=0.
|
|
END DO
|
|
C
|
|
C zero vectors and matrices of the Rybicki formalism
|
|
C
|
|
DO ID=1,ND
|
|
RA(ID)=0.
|
|
RB(ID)=0.
|
|
RC(ID)=0.
|
|
UA(ID)=0.
|
|
UB(ID)=0.
|
|
UC(ID)=0.
|
|
VA(ID)=0.
|
|
VB(ID)=0.
|
|
VC(ID)=0.
|
|
VR(ID)=0.
|
|
WR(ID)=0.
|
|
AL(ID)=0.
|
|
AU(ID)=0.
|
|
VAL(ID)=0.
|
|
DO ID1=1,ND
|
|
WM(ID1,ID)=0.
|
|
VAU(ID1,ID)=0.
|
|
END DO
|
|
END DO
|
|
C
|
|
if(ioptab.lt.-1) call setdrt
|
|
c
|
|
DO IJ=1,NFREQ
|
|
FR=FREQ(IJ)
|
|
W0=W0E(IJ)
|
|
c IF(IOPTAB.GE.0) THEN
|
|
CALL OPACTR(IJ)
|
|
c ELSE
|
|
c CALL OPACFD(IJ)
|
|
c END IF
|
|
CALL RTEFR1(IJ)
|
|
CALL ALIFR1(IJ)
|
|
CALL ROSSTD(IJ)
|
|
C
|
|
CALL RYBMAT(IJ)
|
|
CALL TRIDAG(RA,RB,RC,VR,AL,ND)
|
|
ID=1
|
|
VAL(ID)=VAL(ID)+VB(ID)*AL(ID)+VC(ID)*AL(ID+1)
|
|
DO ID=2,ND-1
|
|
VAL(ID)=VAL(ID)+VA(ID)*AL(ID-1)+VB(ID)*AL(ID)+
|
|
* VC(ID)*AL(ID+1)
|
|
END DO
|
|
ID=ND
|
|
VAL(ID)=VAL(ID)+VA(ID)*AL(ID-1)+VB(ID)*AL(ID)
|
|
C
|
|
DO IDC=1,ND
|
|
DO ID=1,ND
|
|
UCOL(ID)=0.
|
|
END DO
|
|
UCOL(IDC)=UB(IDC)
|
|
IF(IDC.GT.1) UCOL(IDC-1)=UC(IDC-1)
|
|
IF(IDC.LT.ND) UCOL(IDC+1)=UA(IDC+1)
|
|
CALL TRIDAG(RA,RB,RC,UCOL,AU,ND)
|
|
ID=1
|
|
VAU(ID,IDC)=VAU(ID,IDC)+VB(ID)*AU(ID)+VC(ID)*AU(ID+1)
|
|
DO ID=2,ND-1
|
|
VAU(ID,IDC)=VAU(ID,IDC)+
|
|
* VA(ID)*AU(ID-1)+VB(ID)*AU(ID)+VC(ID)*AU(ID+1)
|
|
END DO
|
|
ID=ND
|
|
VAU(ID,IDC)=VAU(ID,IDC)+VA(ID)*AU(ID-1)+VB(ID)*AU(ID)
|
|
END DO
|
|
END DO
|
|
C
|
|
DO ID=1,ND
|
|
ABROSD(ID)= SUMDPL(ID)/ABROSD(ID)
|
|
FCOOL(ID)=REINT(ID)*FCOOLI(ID)*DENS(ID)-REDIF(ID)*FLFIX(ID)
|
|
END DO
|
|
if(ioptab.lt.0) CALL ROSSTD(0)
|
|
C
|
|
C final evaluation of matrices and the global inversion
|
|
C
|
|
CALL RYBENE
|
|
C
|
|
do id=1,70,34
|
|
write(6,603) id,wm(id,id),wm(id,id+1),wr(id)
|
|
end do
|
|
603 format('rybene ',8x,i4,1p4e11.3)
|
|
DO ID=1,ND
|
|
DO IDC=1,ND
|
|
WM(ID,IDC)=WM(ID,IDC)-VAU(ID,IDC)
|
|
END DO
|
|
WR(ID)=WR(ID)-VAL(ID)
|
|
END DO
|
|
CALL LINEQS(WM,WR,CHANGT,ND,MDEPTH)
|
|
C
|
|
IF(.NOT.LTE.and.ifryb.gt.1) THEN
|
|
lte=.true.
|
|
iflev0=iflev
|
|
iflev=1
|
|
do i=1,nlevel
|
|
imodl(i)=imodl0(i)
|
|
end do
|
|
call levset
|
|
do id=1,nd
|
|
call steqeq(id,pop1,0)
|
|
do i=1,nlevel
|
|
babs(i,id)=un
|
|
if(pop1(i).gt.0.) babs(i,id)=popul(i,id)/pop1(i)
|
|
end do
|
|
end do
|
|
lte=.false.
|
|
iflev=iflev0
|
|
do i=1,nlevel
|
|
imodl(i)=imodl0(i)
|
|
end do
|
|
call levset
|
|
end if
|
|
C
|
|
IF(.NOT.LTE.and.ifryb.gt.2) THEN
|
|
DO IJ=1,NFREQ
|
|
CALL OPACTR(IJ)
|
|
CALL RTEFR1(IJ)
|
|
CALL ALIFR1(IJ)
|
|
CALL ROSSTD(IJ)
|
|
CALL RYBMAT(IJ)
|
|
ID=1
|
|
WR(ID)=VR(ID)-UB(ID)*CHANGT(ID)-UC(ID)*CHANGT(ID+1)
|
|
DO ID=2,ND-1
|
|
WR(ID)=VR(ID)-UA(ID)*CHANGT(ID-1)-UB(ID)*CHANGT(ID)-
|
|
* UC(ID)*CHANGT(ID+1)
|
|
END DO
|
|
ID=ND
|
|
WR(ID)=VR(ID)-UA(ID)*CHANGT(ID-1)-UB(ID)*CHANGT(ID)
|
|
CALL TRIDAG(RA,RB,RC,WR,AL,ND)
|
|
RAD(IJ,ID)=RAD1(ID)+AL(ID)
|
|
END DO
|
|
END IF
|
|
C
|
|
CALL RYBCHN(CHANGT)
|
|
c
|
|
if(.not.lte.and.ifryb.gt.1) then
|
|
lte=.true.
|
|
iflev0=iflev
|
|
iflev=1
|
|
do i=1,nlevel
|
|
imodl(i)=imodl0(i)
|
|
end do
|
|
call levset
|
|
do id=1,nd
|
|
call steqeq(id,pop1,0)
|
|
do i=1,nlevel
|
|
popul(i,id)=pop1(i)*babs(i,id)
|
|
end do
|
|
end do
|
|
lte=.false.
|
|
iflev=iflev0
|
|
do i=1,nlevel
|
|
imodl(i)=imodl0(i)
|
|
end do
|
|
call levset
|
|
end if
|
|
c
|
|
RETURN
|
|
END
|