SUBROUTINE BRTEZ(ID) C ==================== C C The part of matrices A,B,C corresponding to the linearized C radiative transfer equation C i.e. the first NFREQE rows C INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'ATOMIC.FOR' INCLUDE 'MODELQ.FOR' INCLUDE 'ALIPAR.FOR' INCLUDE 'ARRAY1.FOR' PARAMETER (XCON=8.0935D-21,YCON=1.68638E-10) PARAMETER (SIXTH=UN/6.D0, * THIRD=UN/3.D0) C IF(NFREQE.LE.0) RETURN ispl=isplin if(isplin.ge.5) isplin=isplin-5 NHE=NFREQE+INHE NRE=NFREQE+INRE NPC=NFREQE+INPC NSE=NFREQE+INSE-1 NMP=NFREQE+INMP C GP=0. GN=UN IF(INMP.GT.0) THEN GP=UN GN=0. END IF c c in the case of Compton scattering - boundary condition c for the highest frequency C IJ1=1 if(icompt.gt.0.and.icombc.gt.0.and.ijex(1).gt.0) then IJ1=2 ij=1 iji=nfreq zj1=exp(-hk*freq(ij)/temp(id)) zj2=exp(-hk*freq(ij+1)/temp(id)) dlt=delj(iji-1,id) if(ichcoo.eq.0) then zj0=un/(hk*sqrt(freq(ij)*freq(ij+1))/temp(id)) zxx=un-3.*zj0+(un-dlt)*zj1+dlt*zj2 combid=zj0/dlnfr(iji-1)+(un-dlt)*zxx comaid=-zj0/dlnfr(iji-1)+dlt*zxx else e2=ycon*temp(id) zxx0=xcon*freq(ij)*(un+zj1)-3.*e2 zxxm=xcon*freq(ij+1)*(un+zj2)-3.*e2 zxx=(un-dlt)*zxx0+dlt*zxxm combid=e2/dlnfr(iji-1)+(un-dlt)*zxx comaid=-e2/dlnfr(iji-1)+dlt*zxx end if b(ij,ij)=combid b(ij,ij+1)=comaid vecl(ij)=-b(ij,ij)*rad(iji,id)-b(ij,ij+1)*rad(iji-1,id) end if C C C ---------------------------------------- C For ID = 1 - upper boundary condition C ---------------------------------------- C IF(ID.GT.1) GO TO 50 DDP=(ZD(1)-ZD(2))*HALF DO IJ=IJ1,NFREQE IJT=IJFR(IJ) OMEG0=ABSO0(IJ) OMEGP=ABSOP(IJ) DZP=OMEG0+OMEGP DTAUP=DZP*DDP ALF1=(FK0(IJ)*RAD0(IJ)-FKP(IJ)*RADP(IJ))/DTAUP CHIEL0=SCAT0(IJ) CHIELP=SCATP(IJ) S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ) BS=HALF*DTAUP CS=0. C2=0. GAM2=0. BET2=0. SP=0. c c additional terms for Compton scattering c if(icompt.gt.0) then call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd) s0=s0+cms end if C IF(MOD(ISPLIN,3).GT.0) THEN C C Spline collocation and/or Hermitian method (ISPLIN=1 or 2) - C both give the same expression for the boundary conditions C BS=DTAUP*THIRD CS=HALF *BS SP=(EMISP(IJ)+CHIELP*RADP(IJ))/ABSOP(IJ) C2=CS/ABSOP(IJ) GAM2=CS*(RADP(IJ)-SP) END IF C C auxiliary quantities C ALF2=BS*(RAD0(IJ)-S0) BET2=ALF2+GAM2 X1=(ALF1-BET2)/DZP B2=(BS+Q0(IJ))/ABSO0(IJ) B1=X1 B1=B1+UU0(IJ)*S0*DM(1)/DENS(1) C1=X1 B1=B1-B2*S0 C1=C1-C2*SP C C *** elements of the IJ-th row of matrices B and C C B(IJ,NRE)=B1*DABT0(IJ)+B2*(DEMT0(IJ)+DST*RAD0(IJ)) C(IJ,NRE)=C1*DABTP(IJ)+C2*(DEMTP(IJ)+DST*RADP(IJ)) B(IJ,NPC)=B1*DABN0(IJ)+ * B2*(DEMN0(IJ)+(DSN+SIGEC(IJT))*RAD0(IJ)) C(IJ,NPC)=C1*DABNP(IJ)+ * C2*(DEMNP(IJ)+(DSN+SIGEC(IJT))*RADP(IJ)) B(IJ,NMP)=B1*DABM0(IJ)+B2*DEMM0(IJ)-GP*RTN C(IJ,NMP)=C1*DABMP(IJ)+C2*DEMMP(IJ)-GP*RTNC DO II=1,NLVEXP B(IJ,NSE+II)=B(IJ,NSE+II)+ * B1*DRCH0(II,IJ)+B2*DRET0(II,IJ) C(IJ,NSE+II)=C(IJ,NSE+II)+ * C1*DRCHP(II,IJ)+C2*DRETP(II,IJ) END DO B(IJ,NFREQE)=0. B(IJ,IJ)=-FK0(IJ)/DTAUP-FH(IJT)-BS*(UN-CHIEL0/ABSO0(IJ))+ * Q0(IJ)*CHIEL0/ABSO0(IJ) C(IJ,NFREQE)=0. C(IJ,IJ)=FKP(IJ)/DTAUP-CS*(UN-CHIELP/ABSOP(IJ)) C C *** the IJ-th element of the rhs vector C VECL(IJ)=ALF1+BET2+FH(IJT)*RAD0(IJ)-S0*Q0(IJ) IF(IWINBL.LT.0) VECL(IJ)=VECL(IJ)-HEXTRD(IJT) c c additional terms for Compton scattering c if(icompt.gt.4) then iji=nfreq-kij(ijt)+1 b(ij,ij)=b(ij,ij)+bs*(cmb+cme) if(iji.gt.1) then ijm=ijex(ijorig(iji-1)) if(ijm.gt.0) b(ij,ijm)=b(ij,ijm)+bs*cma end if if(iji.lt.nfreq) then ijp=ijex(ijorig(iji+1)) if(ijp.gt.0) b(ij,ijp)=b(ij,ijp)+bs*cmc end if if(inre.gt.0) b(ij,nre)=b(ij,nre)+cmd*bs if(inpc.gt.0) b(ij,npc)=b(ij,npc)+cms*bs/elec(id) end if c END DO isplin=ispl go to 500 C C --------------------------------------- C For 1 < ID < ND - normal depth point C --------------------------------------- C 50 DDM=(ZD(ID-1)-ZD(ID))*HALF IF(ID.EQ.ND) GO TO 150 DDP=(ZD(ID)-ZD(ID+1))*HALF DO IJ=IJ1,NFREQE IJT=IJFR(IJ) OMEG0=ABSO0(IJ) OMEGP=ABSOP(IJ) OMEGM=ABSOM(IJ) DZP=OMEG0+OMEGP DZM=OMEG0+OMEGM DTAUP=DZP*DDP DTAUM=DZM*DDM DTAU0=HALF *(DTAUP+DTAUM) FRD=FK0(IJ)*RAD0(IJ) ALF1=(FRD-FKP(IJ)*RADP(IJ))/DTAUP/DTAU0 GAM1=(FRD-FKM(IJ)*RADM(IJ))/DTAUM/DTAU0 BET1=ALF1+GAM1 X1=HALF *BET1/DTAU0 A1=(GAM1+X1*DTAUM)/DZM C1=(ALF1+X1*DTAUP)/DZP B1=A1+C1 BS=UN CHIELM=SCATM(IJ) CHIEL0=SCAT0(IJ) CHIELP=SCATP(IJ) S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ) AS=0. CS=0. A2=0. C2=0. A3=0. C3=0. BET2=0. SM=0. SP=0. c c additional terms for Compton scattering c if(icompt.gt.0) then call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd) s0=s0+cms end if C IF(MOD(ISPLIN,3).EQ.0) GO TO 60 SM=(EMISM(IJ)+RADM(IJ)*CHIELM)/ABSOM(IJ) SP=(EMISP(IJ)+RADP(IJ)*CHIELP)/ABSOP(IJ) IF(ISPLIN.EQ.1) THEN C C spline collocation (ISPLIN=1) C AS=DTAUM/DTAU0*SIXTH CS=DTAUP/DTAU0*SIXTH BS=0.666666666666667D0 ALF2=AS*(RADM(IJ)-SM) GAM2=CS*(RADP(IJ)-SP) BET2=ALF2+GAM2 X =HALF *BET2/DTAU0 A2=(GAM2-X*DTAUM)/DZM C2=(ALF2-X*DTAUP)/DZP ELSE C C Hermitian method (ISPLIN=2) C AS=DTAUP*DTAUP/DTAUM/DTAU0 CS=DTAUM*DTAUM/DTAUP/DTAU0 AL3=(RADP(IJ)-SP-RAD0(IJ)+S0)*SIXTH GA3=(RADM(IJ)-SM-RAD0(IJ)+S0)*SIXTH AV=AL3*CS CV=GA3*AS AS=(UN-HALF *AS)*SIXTH CS=(UN-HALF *CS)*SIXTH BS=UN-AS-CS X=(AV+CV)/DTAU0/4.D0 A2=(X*DTAUM+HALF *CV-AV)/DZM C2=(X*DTAUP+HALF *AV-CV)/DZP BET2=AS*(RADM(IJ)-SM)+CS*(RADP(IJ)-SP) END IF C C auxiliary quantities C B1=B1-(A2+C2) A1=A1-A2 C1=C1-C2 A2=AS/ABSOM(IJ) C2=CS/ABSOP(IJ) A3=A2*SM C3=C2*SP 60 B2=BS/ABSO0(IJ) B3=B2*S0 A1=A1-A3 B1=B1-B3 C1=C1-C3 C C *** elements of the IJ-th row of matrices A, B, and C C A(IJ,NRE)= A1*DABTM(IJ)+A2*(DEMTM(IJ)+DST*RADM(IJ)) B(IJ,NRE)= B1*DABT0(IJ)+B2*(DEMT0(IJ)+DST*RAD0(IJ)) C(IJ,NRE)= C1*DABTP(IJ)+C2*(DEMTP(IJ)+DST*RADP(IJ)) A(IJ,NPC)= A1*DABNM(IJ)+ * A2*(DEMNM(IJ)+(DSN+SIGEC(IJT))*RADM(IJ)) B(IJ,NPC)= B1*DABN0(IJ)+ * B2*(DEMN0(IJ)+(DSN+SIGEC(IJT))*RAD0(IJ)) C(IJ,NPC)= C1*DABNP(IJ)+ * C2*(DEMNP(IJ)+(DSN+SIGEC(IJT))*RADP(IJ)) A(IJ,NMP)= A1*DABMM(IJ)+A2*DEMMM(IJ)-GP*RTNA B(IJ,NMP)= B1*DABM0(IJ)+B2*DEMM0(IJ)-GP*RTN C(IJ,NMP)= C1*DABMP(IJ)+C2*DEMMP(IJ)-GP*RTNC DO II=1,NLVEXP A(IJ,NSE+II)=A(IJ,NSE+II)+ * A1*DRCHM(II,IJ)+A2*DRETM(II,IJ) B(IJ,NSE+II)=B(IJ,NSE+II)+ * B1*DRCH0(II,IJ)+B2*DRET0(II,IJ) C(IJ,NSE+II)=C(IJ,NSE+II)+ * C1*DRCHP(II,IJ)+C2*DRETP(II,IJ) END DO A(IJ,NFREQE)=0. A(IJ,IJ)=FKM(IJ)/DTAUM/DTAU0-AS*(UN-CHIELM/ABSOM(IJ)) B(IJ,NFREQE)=0. B(IJ,IJ)=-FK0(IJ)/DTAU0*(UN/DTAUP+UN/DTAUM)- * BS*(UN-CHIEL0/ABSO0(IJ)) C(IJ,NFREQE)=0. C(IJ,IJ)=FKP(IJ)/DTAUP/DTAU0-CS*(UN-CHIELP/ABSOP(IJ)) C C *** the IJ-th element of the rhs vector C VECL(IJ)=BET1+BET2+BS*(RAD0(IJ)-S0) c c additional terms for Compton scattering c if(icompt.gt.4) then iji=nfreq-kij(ijt)+1 b(ij,ij)=b(ij,ij)+bs*(cmb+cme) if(iji.gt.1) then ijm=ijex(ijorig(iji-1)) if(ijm.gt.0) b(ij,ijm)=b(ij,ijm)+bs*cma end if if(iji.lt.nfreq) then ijp=ijex(ijorig(iji+1)) if(ijp.gt.0) b(ij,ijp)=b(ij,ijp)+bs*cmc end if if(inre.gt.0) b(ij,nre)=b(ij,nre)+cmd*bs if(inpc.gt.0) b(ij,npc)=b(ij,npc)+cms*bs/elec(id) end if c END DO isplin=ispl go to 500 C C -------------------------------------- C For ID=ND - lower boundary condition C -------------------------------------- C 150 CONTINUE IF(IDISK.EQ.0.OR.IFZ0.LT.0) THEN T=TEMP(ID) TM=TEMP(ID-1) HKT=HK/T HKTM=HK/TM C C auxiliary quantites for both options C DO IJ=1,NFREQE IJT=IJFR(IJ) CHIELM=SCATM(IJ) CHIEL0=SCAT0(IJ) OMEGM=ABSOM(IJ) OMEG0=ABSO0(IJ) DZM=OMEG0+OMEGM DTAUM=DZM*DDM FRD=FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ) GAM1=FRD/DTAUM A1=GAM1/DZM AS=0. BS=0. A2=0. B2=0. A3=0. B3=0. ALF2=0. BET2=0. GAM2=0. C C second-order boundary condition C IF(IBC.GT.0.AND.IBC.LT.4) THEN BS=DTAUM*HALF S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ) c c additional terms for Compton scattering c if(icompt.gt.0) then call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd) s0=s0+cms end if C GAM2=BS*(RAD0(IJ)-S0) BET2=GAM2 X1=BET2/DZM A1=A1-X1 B2=BS/ABSO0(IJ) B3=B2*S0 END IF C C auxiliary parameters C FR=FREQ(IJT) FR15=FR*1.D-15 X=HKT*FR EX=EXP(X) XM=HKTM*FR EXM=EXP(XM) PLAN=BN*FR15*FR15*FR15/(EX-UN) IF(INRE.EQ.0.OR.ID.GE.NDRE) THEN PLANM=BN*FR15*FR15*FR15/(EXM-UN) GAM3=(PLAN-PLANM)/DTAUM*THIRD A1=A1-GAM3/DZM GAM1=GAM1-GAM3 END IF C1=A1 B1=C1 A1=A1-A3 B1=B1-B3 C C *** elements of the IJ-th row of matrices A and B C IF(INRE.EQ.0.OR.ID.GE.NDRE) THEN DPLANM=PLANM*XM/TM/(UN-UN/EXM) A(IJ,NRE)=A1*DABTM(IJ)+A2*(DEMTM(IJ)+DST*RADM(IJ))- * DPLANM/DTAUM*THIRD A(IJ,NPC)=A1*DABNM(IJ)+ * A2*(DEMNM(IJ)+(DSN+SIGEC(IJT))*RADM(IJ)) BB=HALF+THIRD/DTAUM DPLAN=PLAN*X/T/(UN-UN/EX) B(IJ,NRE)= B1*DABT0(IJ)+B2*(DEMT0(IJ)+DST*RAD0(IJ))+ * BB*DPLAN B(IJ,NPC)= B1*DABN0(IJ)+ * B2*(DEMN0(IJ)+(DSN+SIGEC(IJT))*RAD0(IJ)) A(IJ,NMP)= A1*DABMM(IJ)+A2*DEMMM(IJ)-GP*RTNA B(IJ,NMP)= B1*DABM0(IJ)+B2*DEMM0(IJ)-GP*RTN DO II=1,NLVEXP A(IJ,NSE+II)=A(IJ,NSE+II)+ * A1*DRCHM(II,IJ)+A2*DRETM(II,IJ) B(IJ,NSE+II)=B(IJ,NSE+II)+ * B1*DRCH0(II,IJ)+B2*DRET0(II,IJ) END DO A(IJ,NFREQE)=0. A(IJ,IJ)=FKM(IJ)/DTAUM-AS*(UN-CHIELM/ABSOM(IJ)) B(IJ,NFREQE)=0. END IF C C *** the IJ-th element of the rhs vector C IF(IBC.EQ.0.OR.IBC.EQ.4) THEN B(IJ,IJ)=B(IJ,IJ)-FK0(IJ)/DTAUM- * BS*(UN-CHIEL0/ABSO0(IJ))-HALF VECL(IJ)=GAM1+BET2-HALF*(PLAN-RAD0(IJ)) ELSE B(IJ,IJ)=B(IJ,IJ)-FK0(IJ)/DTAUM- * BS*(UN-CHIEL0/ABSO0(IJ))-FHD(IJT) VECL(IJ)=GAM1+BET2-HALF*PLAN+FHD(IJT)*RAD0(IJ) END IF c c additional terms for Compton scattering c if(icompt.gt.4) then iji=nfreq-kij(ijt)+1 b(ij,ij)=b(ij,ij)+bs*(cmb+cme) if(iji.gt.1) then ijm=ijex(ijorig(iji-1)) if(ijm.gt.0) b(ij,ijm)=b(ij,ijm)+bs*cma end if if(iji.lt.nfreq) then ijp=ijex(ijorig(iji+1)) if(ijp.gt.0) b(ij,ijp)=b(ij,ijp)+bs*cmc end if if(inre.gt.0) b(ij,nre)=b(ij,nre)+cmd*bs if(inpc.gt.0) b(ij,npc)=b(ij,npc)+cms*bs/elec(id) end if c END DO C ELSE C C -------------------------------------- C For ID=ND - lower boundary condition C -------------------------------------- C C for disks - C lower b.c. expresses just I(taumax,-mu,nu)=I(taumax,+mu,nu) C DO IJ=IJ1,NFREQE IJT=IJFR(IJ) CHIELM=SCATM(IJ) CHIEL0=SCAT0(IJ) OMEGM=ABSOM(IJ) OMEG0=ABSO0(IJ) DZM=OMEG0+OMEGM DTAUM=DZM*DDM FRD=FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ) GAM1=FRD/DTAUM A1=GAM1/DZM AS=0. A2=0. A3=0. ALF2=0. GAM2=0. BS=DTAUM*HALF S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ) c c additional terms for Compton scattering c if(icompt.gt.0) then call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd) s0=s0+cms end if C GAM2=BS*(RAD0(IJ)-S0) BET2=ALF2+GAM2 X1=BET2/DZM A1=A1-X1 B2=BS/ABSO0(IJ) B3=B2*S0 B1=A1 A1=A1-A3 B1=B1-B3 C C *** elements of the IJ-th row of matrix A C A(IJ,NRE)=A1*DABTM(IJ)+A2*(DEMTM(IJ)+DST*RADM(IJ)) A(IJ,NPC)=A1*DABNM(IJ)+ * A2*(DEMNM(IJ)+(DSN+SIGEC(IJT))*RADM(IJ)) A(IJ,NMP)= A1*DABMM(IJ)+A2*DEMMM(IJ)-GP*RTNA DO I=1,NLVEXP A(IJ,NSE+I)=A1*DRCHM(I,IJ)+A2*DRETM(I,IJ) END DO A(IJ,NFREQE)=0. A(IJ,IJ)=FKM(IJ)/DTAUM-AS*(UN-CHIELM/ABSOM(IJ)) C C *** elements of the IJ-th row of matrix B C B(IJ,NRE)=B1*DABT0(IJ)+B2*(DEMT0(IJ)+DST*RAD0(IJ)) B(IJ,NPC)=B1*DABN0(IJ)+ * B2*(DEMN0(IJ)+(DSN+SIGEC(IJT))*RAD0(IJ)) B(IJ,NMP)= B1*DABM0(IJ)+B2*DEMM0(IJ)-GP*RTN DO I=1,NLVEXP B(IJ,NSE+I)=B1*DRCH0(I,IJ)+B2*DRET0(I,IJ) END DO B(IJ,NFREQE)=0. B(IJ,IJ)=-FK0(IJ)/DTAUM-BS*(UN-CHIEL0/ABSO0(IJ)) C C *** the IJ-th element of the rhs vector C VECL(IJ)=GAM1+BET2 c c additional terms for Compton scattering c if(icompt.gt.4) then iji=nfreq-kij(ijt)+1 b(ij,ij)=b(ij,ij)+bs*(cmb+cme) if(iji.gt.1) then ijm=ijex(ijorig(iji-1)) if(ijm.gt.0) b(ij,ijm)=b(ij,ijm)+bs*cma end if if(iji.lt.nfreq) then ijp=ijex(ijorig(iji+1)) if(ijp.gt.0) b(ij,ijp)=b(ij,ijp)+bs*cmc end if if(inre.gt.0) b(ij,nre)=b(ij,nre)+cmd*bs if(inpc.gt.0) b(ij,npc)=b(ij,npc)+cms*bs/elec(id) end if c END DO END IF isplin=ispl 500 CONTINUE c c zeroing radiation field for very low intensities (if required) c if(radzer.gt.0) then c C find the peak in nu*rad_nu: c radsum=0. DO IJ=IJ1,NFREQE radsum=max(freq(ij)*radex(ij,id),radsum) END DO C C if much smaller than peak in nu*rad_nu, then set to zero: C DO IJ=IJ1,NFREQE if(freq(ij)*radex(ij,id).lt.radzer*radsum) then do ii=1,nn0 a(ij,ii)=0. b(ij,ii)=0. c(ij,ii)=0. end do vecl(ij)=0. b(ij,ij)=un end if end do end if isplin=ispl c RETURN END