SUBROUTINE RTE
C
C     solution of the radiative transfer equation by Feautrier method
C
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      INCLUDE 'SYNTHP.FOR'
      INCLUDE 'LINDAT.FOR'
      DIMENSION D(3,3,MDEPTH),ANU(3,MDEPTH),AANU(MDEPTH),DDD(MDEPTH),
     *       AA(3,3),BB(3,3),CC(3,3),VL(3),AMU(3),WTMU(3),
     *       DT(MDEPTH),TAU(MDEPTH),
     *       RDD(MDEPTH),FKK(MDEPTH),ST0(MDEPTH),SS0(MDEPTH),
     *       RINT(MDEPTH,MMU)
      CHARACTER*4 TYPION(9) 
      COMMON/RTEOPA/CH(MFREQ,MDEPTH),ET(MFREQ,MDEPTH),
     *              SC(MFREQ,MDEPTH)
      COMMON/EMFLUX/FLUX(MFREQ),FLUXC(MFREQC)
      COMMON/BLAPAR/RELOP,SPACE0,CUTOF0,TSTD,DSTD,ALAMC
      COMMON/CTRFUN/CINT1(MDEPTH),CINT2(MDEPTH),
     *   CTRI(MDEPTH),CTRR(MDEPTH),XKAR(MDEPTH),
     *   ABXLI(MFREQ),EMXLI(MFREQ),IJCTR(MFREQ)
      COMMON/REFDEP/IREFD(MFREQ)
      COMMON/CENTRL/ZND,IFZ0
      PARAMETER (UN=1.D0, HALF=0.5D0)
      PARAMETER (THIRD=UN/3., QUART=UN/4., SIXTH=UN/6.D0)
      PARAMETER (TAUREF = 0.6666666666667)
      DATA AMU/.887298334620742D0,.5D0,.112701665379258D0/,
     1     WTMU/.277777777777778D0,.444444444444444D0,.277777777777778D0
     1         /
      DATA TYPION /' I  ',' II ',' III',' IV ',' V  ',
     *             ' VI ',' VII','VIII',' IX '/
C
      NMU=3
      ND1=ND-1
C
C     Overall loop over frequencies 
C
      DO IJ=1,NFREQ
      TAUMIN=CH(IJ,1)/DENS(1)*DM(1)*HALF
      TAU(1)=TAUMIN
      IREF=1
      DO I=1,ND1
         DT(I)=(DM(I+1)-DM(I))*(CH(IJ,I+1)/DENS(I+1)+CH(IJ,I)/DENS(I))*
     *         HALF
         ST0(I)=ET(IJ,I)/CH(IJ,I)
         SS0(I)=-SC(IJ,I)/CH(IJ,I)
         TAU(I+1)=TAU(I)+DT(I)
         IF(TAU(I).LE.TAUREF.AND.TAU(I+1).GT.TAUREF) IREF=I
      END DO
      IREFD(IJ)=IREF
      ST0(ND)=ET(IJ,ND)/CH(IJ,ND)
      SS0(ND)=-SC(IJ,ND)/CH(IJ,ND)
      FR=FREQ(IJ)
      BNU=BN*(FR*1.E-15)**3
      PLAND=BNU/(EXP(HK*FR/TEMP(ND  ))-UN)
      DPLAN=BNU/(EXP(HK*FR/TEMP(ND-1))-UN)
      DPLAN=(PLAND-DPLAN)/DT(ND1)
C
C   +++++++++++++++++++++++++++++++++++++++++
C   FIRST PART  -  VARIABLE EDDINGTON FACTORS
C   +++++++++++++++++++++++++++++++++++++++++
C
      ALB1=0.
      DO I=1,NMU
C
C   ************************
C   UPPER BOUNDARY CONDITION
C   ************************
C
         ID=1
         DTP1=DT(1)
         Q0=0.
         P0=0.
C
C       allowance for non-zero optical depth at the first depth point
C
         TAMM=TAUMIN/AMU(I)
         IF(TAMM.GT.0.01) THEN
            P0=UN-EXP(-TAMM)
          ELSE
            P0=TAMM*(UN-HALF*TAMM*(UN-TAMM*THIRD*(UN-QUART*TAMM)))
         END IF
         EX=UN-P0
         Q0=Q0+P0*AMU(I)*WTMU(I)
C
         DIV=DTP1/AMU(I)*THIRD
         VL(I)=DIV*(ST0(ID)+HALF*ST0(ID+1))+ST0(ID)*P0
         DO J=1,NMU
            BB(I,J)=SS0(ID)*WTMU(J)*(DIV+P0)-ALB1*WTMU(J)
            CC(I,J)=-HALF*DIV*SS0(ID+1)*WTMU(J)
         END DO
         BB(I,I)=BB(I,I)+AMU(I)/DTP1+UN+DIV
         CC(I,I)=CC(I,I)+AMU(I)/DTP1-HALF*DIV
         ANU(I,ID)=0.
      END DO
C
C     Matrix inversion: instead of calling MATINV, a very fast inlined 
C     routine MINV3 for a specific 3 x 3 matrix inversion
C
C     CALL MATINV(BB,NMU,3)
C
C     ******************************
      BB(2,1)=BB(2,1)/BB(1,1)
      BB(2,2)=BB(2,2)-BB(2,1)*BB(1,2)
      BB(2,3)=BB(2,3)-BB(2,1)*BB(1,3)
      BB(3,1)=BB(3,1)/BB(1,1)
      BB(3,2)=(BB(3,2)-BB(3,1)*BB(1,2))/BB(2,2)
      BB(3,3)=BB(3,3)-BB(3,1)*BB(1,3)-BB(3,2)*BB(2,3)
C
      BB(3,2)=-BB(3,2)
      BB(3,1)=-BB(3,1)-BB(3,2)*BB(2,1)
      BB(2,1)=-BB(2,1)
C
      BB(3,3)=UN/BB(3,3)
      BB(2,3)=-BB(2,3)*BB(3,3)/BB(2,2)
      BB(2,2)=UN/BB(2,2)
      BB(1,3)=-(BB(1,2)*BB(2,3)+BB(1,3)*BB(3,3))/BB(1,1)
      BB(1,2)=-BB(1,2)*BB(2,2)/BB(1,1)
      BB(1,1)=UN/BB(1,1)
C
      BB(1,1)=BB(1,1)+BB(1,2)*BB(2,1)+BB(1,3)*BB(3,1)
      BB(1,2)=BB(1,2)+BB(1,3)*BB(3,2)
      BB(2,1)=BB(2,2)*BB(2,1)+BB(2,3)*BB(3,1)
      BB(2,2)=BB(2,2)+BB(2,3)*BB(3,2)
      BB(3,1)=BB(3,3)*BB(3,1)
      BB(3,2)=BB(3,3)*BB(3,2)
C     ******************************
C
      DO I=1,NMU
         DO J=1,NMU
            S=0.
            DO K=1,NMU
               S=S+BB(I,K)*CC(K,J)
            END DO
            D(I,J,ID)=S
            ANU(I,1)=ANU(I,1)+BB(I,J)*VL(J)
         END DO
      END DO
C
C   *******************
C   NORMAL DEPTH POINTS
C   *******************
C
      DO ID=2,ND1
         DTM1=DTP1
         DTP1=DT(ID)
         DT0=HALF*(DTM1+DTP1)
         AL=UN/DTM1/DT0
         GA=UN/DTP1/DT0
         BE=AL+GA
         A=(UN-HALF*AL*DTP1*DTP1)*SIXTH
         C=(UN-HALF*GA*DTM1*DTM1)*SIXTH
         B=UN-A-C
         VL0=A*ST0(ID-1)+B*ST0(ID)+C*ST0(ID+1)
         DO I=1,NMU
            DO J=1,NMU
               AA(I,J)=-A*SS0(ID-1)*WTMU(J)
               CC(I,J)=-C*SS0(ID+1)*WTMU(J)
               BB(I,J)=B*SS0(ID)*WTMU(J)
            END DO
         END DO
         DO I=1,NMU
            DIV=AMU(I)**2
            VL(I)=VL0
            AA(I,I)=AA(I,I)+DIV*AL-A
            CC(I,I)=CC(I,I)+DIV*GA-C
            BB(I,I)=BB(I,I)+DIV*BE+B
         END DO
         DO I=1,NMU
            S1=0.
            DO J=1,NMU
               S=0.
               S1=S1+AA(I,J)*ANU(J,ID-1)
               DO K=1,NMU
                  S=S+AA(I,K)*D(K,J,ID-1)
               END DO
               BB(I,J)=BB(I,J)-S
            END DO
            VL(I)=VL(I)+S1
         END DO
C
C     Matrix inversion: instead of calling MATINV, a very fast inlined 
C     routine MINV3 for a specific 3 x 3 matrix inversion
C
C     CALL MATINV(BB,NMU,3)
C
C     ******************************
      BB(2,1)=BB(2,1)/BB(1,1)
      BB(2,2)=BB(2,2)-BB(2,1)*BB(1,2)
      BB(2,3)=BB(2,3)-BB(2,1)*BB(1,3)
      BB(3,1)=BB(3,1)/BB(1,1)
      BB(3,2)=(BB(3,2)-BB(3,1)*BB(1,2))/BB(2,2)
      BB(3,3)=BB(3,3)-BB(3,1)*BB(1,3)-BB(3,2)*BB(2,3)
C
      BB(3,2)=-BB(3,2)
      BB(3,1)=-BB(3,1)-BB(3,2)*BB(2,1)
      BB(2,1)=-BB(2,1)
C
      BB(3,3)=UN/BB(3,3)
      BB(2,3)=-BB(2,3)*BB(3,3)/BB(2,2)
      BB(2,2)=UN/BB(2,2)
      BB(1,3)=-(BB(1,2)*BB(2,3)+BB(1,3)*BB(3,3))/BB(1,1)
      BB(1,2)=-BB(1,2)*BB(2,2)/BB(1,1)
      BB(1,1)=UN/BB(1,1)
C
      BB(1,1)=BB(1,1)+BB(1,2)*BB(2,1)+BB(1,3)*BB(3,1)
      BB(1,2)=BB(1,2)+BB(1,3)*BB(3,2)
      BB(2,1)=BB(2,2)*BB(2,1)+BB(2,3)*BB(3,1)
      BB(2,2)=BB(2,2)+BB(2,3)*BB(3,2)
      BB(3,1)=BB(3,3)*BB(3,1)
      BB(3,2)=BB(3,3)*BB(3,2)
C     ******************************
C
         DO I=1,NMU
            ANU(I,ID)=0.
            DO J=1,NMU
               S=0.
               DO K=1,NMU
                  S=S+BB(I,K)*CC(K,J)
               END DO
               D(I,J,ID)=S
               ANU(I,ID)=ANU(I,ID)+BB(I,J)*VL(J)
            END DO
         END DO
      END DO
C
C   ************
C   LOWER BOUNDARY CONDITION
C   ************
C
      ID=ND
C
C     First option:
C     b.c. is different from stellar atmospheres; expresses symmetry
C     at the central plane   I(taumax,-mu,nu)=I(taumax,+mu,nu)
C
      IF(IFZ0.EQ.0) THEN
         B=DTP1*HALF
         A=0.
         DO I=1,NMU
            BI=B/AMU(I)
            AI=A/AMU(I)
            VL(I)=ST0(ID)*BI+ST0(ID-1)*AI
            DO J=1,NMU
               AA(I,J)=-AI*SS0(ID-1)*WTMU(J)
               BB(I,J)=BI*SS0(ID)*WTMU(J)
            END DO
            AA(I,I)=AA(I,I)+AMU(I)/DTP1-AI
            BB(I,I)=BB(I,I)+AMU(I)/DTP1+BI
         END DO
         DO I=1,NMU
            S1=0.
            DO J=1,NMU
               S=0.
               S1=S1+AA(I,J)*ANU(J,ID-1)
               DO K=1,NMU
                  S=S+AA(I,K)*D(K,J,ID-1)
               END DO
               BB(I,J)=BB(I,J)-S
            END DO
            VL(I)=VL(I)+S1
         END DO
C
C     Second option:
C     b.c. is the same as in stellar atmospheres - the last depth point
C     is not at the central plane
C
      ELSE
         DO I=1,NMU
            AA(I,I)=AMU(I)/DTP1
            VL(I)=PLAND+AMU(I)*DPLAN+AA(I,I)*ANU(I,ID-1)
            DO J=1,NMU
               BB(I,J)=-AA(I,I)*D(I,J,ID-1)
            END DO
            BB(I,I)=BB(I,I)+AA(I,I)+UN
         END DO
      END IF
C
C     Matrix inversion: instead of calling MATINV, a very fast inlined 
C     routine MINV3 for a specific 3 x 3 matrix inversion
C
C     CALL MATINV(BB,NMU,3)
C
C     ******************************
      BB(2,1)=BB(2,1)/BB(1,1)
      BB(2,2)=BB(2,2)-BB(2,1)*BB(1,2)
      BB(2,3)=BB(2,3)-BB(2,1)*BB(1,3)
      BB(3,1)=BB(3,1)/BB(1,1)
      BB(3,2)=(BB(3,2)-BB(3,1)*BB(1,2))/BB(2,2)
      BB(3,3)=BB(3,3)-BB(3,1)*BB(1,3)-BB(3,2)*BB(2,3)
C
      BB(3,2)=-BB(3,2)
      BB(3,1)=-BB(3,1)-BB(3,2)*BB(2,1)
      BB(2,1)=-BB(2,1)
C
      BB(3,3)=UN/BB(3,3)
      BB(2,3)=-BB(2,3)*BB(3,3)/BB(2,2)
      BB(2,2)=UN/BB(2,2)
      BB(1,3)=-(BB(1,2)*BB(2,3)+BB(1,3)*BB(3,3))/BB(1,1)
      BB(1,2)=-BB(1,2)*BB(2,2)/BB(1,1)
      BB(1,1)=UN/BB(1,1)
C
      BB(1,1)=BB(1,1)+BB(1,2)*BB(2,1)+BB(1,3)*BB(3,1)
      BB(1,2)=BB(1,2)+BB(1,3)*BB(3,2)
      BB(2,1)=BB(2,2)*BB(2,1)+BB(2,3)*BB(3,1)
      BB(2,2)=BB(2,2)+BB(2,3)*BB(3,2)
      BB(3,1)=BB(3,3)*BB(3,1)
      BB(3,2)=BB(3,3)*BB(3,2)
C     ******************************
C
      DO I=1,NMU
         ANU(I,ID)=0.
         DO J=1,NMU
            D(I,J,ID)=0.
            ANU(I,ID)=ANU(I,ID)+BB(I,J)*VL(J)
         END DO
      END DO
C
C   ************
C   BACKSOLUTION
C   ************
C
      ID=ND
      FKK(ND)=THIRD
      AJ=0.
      AK=0.
      DO I=1,NMU
         RMU=WTMU(I)*ANU(I,ID)
         AJ=AJ+RMU
         AK=AK+RMU*AMU(I)*AMU(I)
      END DO
      RDD(ID)=AJ
      FKK(ND)=AK/AJ
      DO ID=ND-1,1,-1
         DO I=1,NMU
            DO J=1,NMU
               ANU(I,ID)=ANU(I,ID)+D(I,J,ID)*ANU(J,ID+1)
            END DO
         END DO
         AJ=0.
         AK=0.
         DO I=1,NMU
            DIV=WTMU(I)*ANU(I,ID)
            AJ=AJ+DIV
            AK=AK+DIV*AMU(I)**2
         END DO
         FKK(ID)=AK/AJ
      END DO
C
C     surface Eddington actor
C
      AH=0.
      DO I=1,NMU
         AH=AH+WTMU(I)*AMU(I)*ANU(I,1)
      END DO
      FH=AH/AJ-HALF*ALB1
C
c     FKK(ND)=THIRD
C
C
C   +++++++++++++++++++++++++++++++++++++++++
C   SECOND PART  -  DETERMINATION OF THE MEAN INTENSITIES
C   RECALCULATION OF THE TRANSFER EQUATION WITH GIVEN EDDINGTON FACTORS
C   +++++++++++++++++++++++++++++++++++++++++
C
      DTP1=DT(1)
      DIV=DTP1*THIRD
      BBB=FKK(1)/DTP1+FH+DIV+SS0(1)*(DIV+Q0)
      CCC=FKK(2)/DTP1-HALF*DIV*(UN+SS0(2))
      VLL=DIV*(ST0(1)+HALF*ST0(2))+ST0(1)*Q0
      AANU(1)=VLL/BBB
      DDD(1)=CCC/BBB
      DO ID=2,ND1
         DTM1=DTP1
         DTP1=DT(ID)
         DT0=HALF*(DTP1+DTM1)
         AL=UN/DTM1/DT0
         GA=UN/DTP1/DT0
         A=(UN-HALF*DTP1*DTP1*AL)*SIXTH
         C=(UN-HALF*DTM1*DTM1*GA)*SIXTH
         AAA=AL*FKK(ID-1)-A*(UN+SS0(ID-1))
         CCC=GA*FKK(ID+1)-C*(UN+SS0(ID+1))
         BBB=(AL+GA)*FKK(ID)+(UN-A-C)*(UN+SS0(ID))
         VLL=A*ST0(ID-1)+C*ST0(ID+1)+(UN-A-C)*ST0(ID)
         BBB=BBB-AAA*DDD(ID-1)
         DDD(ID)=CCC/BBB
         AANU(ID)=(VLL+AAA*AANU(ID-1))/BBB
      END DO
C
C     Lower boundary condition
C     1.option -  different from stellar atmospheres
C
      IF(IFZ0.EQ.0) THEN
         B=DTP1*HALF
         BBB=FKK(ND)/DTP1+B*(UN+SS0(ND))
         AAA=FKK(ND-1)/DTP1
         VLL=B*ST0(ND)
       ELSE
C
C     Lower boundary condition
C     2.option - stellar atmospheric
C
         BBB=FKK(ND)/DTP1+HALF
         AAA=FKK(ND1)/DTP1
         VLL=HALF*PLAND+DPLAN*THIRD
      END IF
      BBB=BBB-AAA*DDD(ND1)
      RDD(ND)=(VLL+AAA*AANU(ND1))/BBB
      DO IID=1,ND1
         ID=ND-IID
         RDD(ID)=AANU(ID)+DDD(ID)*RDD(ID+1)
      END DO
      FLUX(IJ)=FH*RDD(1)
C
C     if needed (if iprin.ge.3), output of interesting physical
C     quantities at the monochromatic optical depth  tau(nu)=2/3
C
      IF(IPRIN.ge.3) THEN
      T0=LOG(TAU(IREF+1)/TAU(IREF))
      X0=LOG(TAU(IREF+1)/TAUREF)/T0
      X1=LOG(TAUREF/TAU(IREF))/T0
      DMREF=EXP(LOG(DM(IREF))*X0+LOG(DM(IREF+1))*X1)
      TREF=EXP(LOG(TEMP(IREF))*X0+LOG(TEMP(IREF+1))*X1)
      STREF=EXP(LOG(ST0(IREF))*X0+LOG(ST0(IREF+1))*X1)
      SCREF=EXP(LOG(-SS0(IREF))*X0+LOG(-SS0(IREF+1))*X1)
      SSREF=EXP(LOG(-SS0(IREF)*RDD(IREF))*X0+
     *           LOG(-SS0(IREF+1)*RDD(IREF+1))*X1)
      SREF=STREF+SSREF
      ALM=2.997925E18/FREQ(IJ)
      WRITE(96,636) IJ,ALM,IREF,DMREF,TREF,SCREF,STREF,SSREF,SREF
  636 FORMAT(1H ,I3,F10.3,I4,1PE10.3,0PF10.1,1X,1P3E10.3,E11.3)
      END IF
C
C   THIRD PART  -  DETERMINATION OF THE SPECIFIC INTENSITIES
C   RECALCULATION OF THE TRANSFER EQUATION WITH GIVEN SOURCE FUNCTION
C
      if(iflux.eq.0) return   
      DO IMU=1,NMU0
      ANX=ANGL(IMU)
      DTP1=DT(1)
      DIV=DTP1*THIRD/ANX
C
      TAMM=TAUMIN/ANX
      IF(TAMM.LT.0.01) THEN
         P0=TAMM*(UN-HALF*TAMM*(UN-TAMM*THIRD*(UN-QUART*TAMM)))
       ELSE
         P0=UN-EXP(-TAMM)
      END IF
C
      BBB=ANX/DTP1+UN+DIV
      CCC=ANX/DTP1-HALF*DIV
      VLL=(DIV+P0)*(ST0(1)-SS0(1)*RDD(1))
     *    +HALF*DIV*(ST0(2)-SS0(2)*RDD(2))
      AANU(1)=VLL/BBB
      DDD(1)=CCC/BBB
      DIV=ANX*ANX
      DO ID=2,ND1
         DTM1=DT(ID-1)
         DTP1=DT(ID)
         DT0=HALF*(DTP1+DTM1)
         AL=UN/DTM1/DT0
         GA=UN/DTP1/DT0
         A=(UN-HALF*DTP1*DTP1*AL)*SIXTH
         C=(UN-HALF*DTM1*DTM1*GA)*SIXTH
         AAA=DIV*AL-A
         CCC=DIV*GA-C
         BBB=DIV*(AL+GA)+UN-A-C
         VLL=A*(ST0(ID-1)-SS0(ID-1)*RDD(ID-1))+
     *       C*(ST0(ID+1)-SS0(ID+1)*RDD(ID+1))+
     *       (UN-A-C)*(ST0(ID)-SS0(ID)*RDD(ID))
         BBB=BBB-AAA*DDD(ID-1)
         DDD(ID)=CCC/BBB
         AANU(ID)=(VLL+AAA*AANU(ID-1))/BBB
      END DO
C
C     Lower boundary condition
C     1.option -  different from stellar atmospheres
C
      IF(IFZ0.EQ.0) THEN
         B=DTP1*HALF/ANX
         BBB=ANX/DTP1+B*(UN+SS0(ND))
         AAA=ANX/DTP1
         VLL=B*ST0(ND)
       ELSE
C
C     Lower boundary condition
C     2.option - stellar atmospheric
C
         AAA=ANX/DTP1
         BBB=AAA+UN
         VLL=PLAND+ANX*DPLAN
      END IF
C
      RINT(ND,IMU)=(VLL+AAA*AANU(ND1))/(BBB-AAA*DDD(ND1))
      DO IID=1,ND1
         ID=ND-IID
         RINT(ID,IMU)=AANU(ID)+DDD(ID)*RINT(ID+1,IMU)
      END DO
      END DO
c
      FLX=0.
      DO IMU=1,NMU0
         RINT(1,IMU)=RINT(1,IMU)/HALF
         FLX=FLX+ANGL(IMU)*WANGL(IMU)*RINT(1,IMU)
      END DO
      FLX=FLX*HALF
c     FLUX(IJ)=FLX
C
C     output of emergent specific intensities to Unit 10 
C     and 18 (continuum)
C
      IF(IJ.GT.2) THEN
         WRITE(10,641) WLAM(IJ),FLX,(RINT(1,IMU),IMU=1,NMU0)
       ELSE
         WRITE(18,641) WLAM(IJ),FLX,(RINT(1,IMU),IMU=1,NMU0)
      END IF
c
      if(iprin.eq.4) then
c
c     compute contribution function C_i (ctri) and C_r (ctrr)
c     following Magain (1986, A&A 163, 135)
c
       if(ijctr(ij).gt.0) then
         xfr0=(freq(ij)-freq(2))/(freq(1)-freq(2))
         tauc=ch(1,1)/dens(1)*dm(1)*half
         do id=1,nd
           chc1=ch(1,id)
           chc2=ch(2,id)
           chcc=chc2+xfr0*(chc1-chc2)
           etc1=et(1,id)
           etc2=et(2,id)
           etcc=etc2+xfr0*(etc1-etc2)
           stcc=etcc/chcc
           cint=cint2(id)+xfr0*(cint1(id)-cint2(id))
           avx=(chc1+chc2)*0.5*relop
           call linop(id,abxli,emxli,avx)
           sli0=emxli(ij)/abxli(ij)
           abt0=ch(ij,id)
           emt0=et(ij,id)
           stt0=emt0/abt0
           Xkar(id)=abxli(ij)+chcc*stcc/cint
           ctri(id)=tauc*abt0/chc1*stt0*exp(-tau(id))
           if(tau(id).gt.70.) ctri(id)=0.
           ctrr(id)=tauc/chc1*abxli(ij)*(un-sli0/cint)
           if(id.lt.nd) then
             dtc=(ch(1,id+1)/dens(id+1)+ch(1,id)/dens(id))
             tauc=tauc+half*dtc*(dm(id+1)-dm(id))
           endif
         end do
         taurs=Xkar(1)/dens(1)*dm(1)*half
         xcti=ctri(1)*half*(dm(2)-dm(1))
         xctr=ctrr(1)*half*(dm(2)-dm(1))
         do i=1,nd-1
           ctrr(i)=ctrr(i)*exp(-taurs)
           if(i.eq.1) xctr=xctr*exp(-taurs)
           if(i.gt.1) then
             xcti=xcti+ctri(i)*half*(dm(i+1)-dm(i-1))
             xctr=xctr+ctrr(i)*half*(dm(i+1)-dm(i-1))
           endif
           if(taurs.gt.70.) ctrr(i)=0.
           dtrs=(dm(i+1)-dm(i))*(Xkar(i+1)/dens(i+1)+Xkar(i)/dens(i))
           taurs=taurs+half*dtrs
         end do
         ctrr(nd)=0.
         alam=2.997925d18/freq(ij)
         il0=ijctr(ij)
         il=indlin(il0)
         iat=indat(il)/100
         ion=mod(indat(il),100)
         write(97,376) il,alam,typat(iat),typion(ion),iref,dmref,tref
  376    format(i5,f11.4,2x,2a4,i8,1pe12.4,0pf10.1)
         do id=1,nd
           ctrip=ctri(id)/xcti
           ctrrp=ctrr(id)/xctr
           write(97,377) id,dm(id),tau(id),ctrip,ctrrp
  377      format(i4,1p4e12.4)
         end do
       else if(ij.eq.1) then
         do id=1,nd
           cint1(id)=rint(id,nmu0)
         end do
       else if(ij.eq.2) then
         do id=1,nd
           cint2(id)=rint(id,nmu0)
         end do
       endif
      endif
  641 FORMAT(1H ,f10.3,1pe15.5/(1P5E15.5))
c
c     end of the global loop over frequencies
c
      END DO
      RETURN
      END