SUBROUTINE PRINC
C     ================
C
C     Auxiliary output routine, which enables printing 
C     more detailed information about chosen individual transitions
C
C     Input: NCT     - number of transitions for which information is
C                      printed
C            ICTR(I) - index of the I-the considered transition
C            INFR(I) - index of the characteristic frequency point for
C                      the I-th transition (i.e., standardly, the point just
C                      shortward of edge for continua; or the line center
C                      for lines)
C                    - if =0, then the program choses the characteristic
C                      frequency itself, in the standard manner
C
C     The printed table contains for each transition the following functions
C     of depth:
C
C     depth index;
C     optical depth at the characteristic frequency;
C     partial radiation force (actually, acceleration) - g(rad) -
C         produced by the transition
C     b-factor of the lower level;
C     b-factor of the upper level;
C     upward radiative rate;
C     downward radiative rate;
C     mean intensity at the characteristic frequency;
C     Planck function at the characteristic frequency;
C     total source function at the characteristic frequency;
C     net line (or continuum) source function at the characteristic frequency;
C     net heating rate (i.e heating minus cooling rate);
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ALIPAR.FOR'
      PARAMETER (NPTR=30,CCOR=0.09,SIXTH=UN/6.)
      DIMENSION ST(NPTR,MDEPTH),TAU(NPTR,MDEPTH),
     *          ABST(NPTR,MDEPTH),EMIT(NPTR,MDEPTH),SCTR(NPTR,MDEPTH),
     *          INFR(NPTR),ICTR(NPTR),ABM(NPTR),PRF(MFREQ),
     *          DWF(MFREQ)
C     DIMENSION RADM(MFREQ),ABSM(MFREQ),FAK0(MFREQ),RAD0(MFREQ)
c
      READ(44,*,END=200,ERR=200) NCT
      READ(44,*,END=200,ERR=200) (ICTR(I),I=1,NCT)
      READ(44,*,END=200,ERR=200) (INFR(I),I=1,NCT)
c
      DO IC=1,NCT
        ITR=ICTR(IC)
        IFR=INFR(IC)
        IF(LINE(ITR)) THEN
          IF(IFR.EQ.0)
     *      INFR(IC)=(IFR0(ITR)+IFR1(ITR))/2
        ELSE
          IF(IFR.EQ.0) INFR(IC)=IFR1(ITR)
        END IF
        IFR=INFR(IC)
        CALL OPACF1(IFR)
        DO ID=1,ND
          ABST(IC,ID)=ABSO1(ID)
          EMIT(IC,ID)=EMIS1(ID)
          SCTR(IC,ID)=SCAT1(ID)
        END DO
      END DO
C
      DO ID=1,ND
         T=TEMP(ID)
         ANE=ELEC(ID)
         SQT=SQRT(T)
         ANES=EXP(SIXTH*LOG(ANE))
         AACOR=CCOR*ANES/SQT
         CALL SABOLF(ID)
         DO I=1,NLEVEL
            POP(I)=POPUL(I,ID)
         END DO              
         DO IC=1,NCT
            ITR=ICTR(IC)
            IFR=INFR(IC)
            I=ILOW(ITR)
            J=IUP(ITR)
            IF(LINE(ITR)) THEN
               CALL LINPRO(ITR,ID,PRF)
               GG=G(I)/G(J)
               SG=PRF(IFR)
             ELSE
               GG=ANE*SBF(I)*EXP(-HK*FREQ(IFR)/T)
               QZ=IZ(IEL(I))
               MW=MCDW(ITR)
               CALL DWNFR(MW,NFREQ,FR0(ITR),AACOR,ANE,QZ,FREQ,DWF)
               SG=CROSS(ITRA(J,I),IFR)*DWF(IFR)
            END IF
            IJE=IJEX(IFR)
            ESCT=SCTR(IC,ID)*RAD(IFR,ID)
            ST(IC,ID)=(EMIT(IC,ID)+ESCT)/ABST(IC,ID)
            IF(ID.EQ.1) THEN
               TAU(IC,ID)=HALF *ABST(IC,ID)/DENS(ID)*DM(ID)
             ELSE
               TAU(IC,ID)=TAU(IC,ID-1)+(ABST(IC,ID)/DENS(ID)+ABM(IC))*
     *                 (DM(ID)-DM(ID-1))*HALF
            END IF
            ABM(IC)=ABST(IC,ID)/DENS(ID)
         END DO
      END DO
C
      DO IC=1,NCT
         ITR=ICTR(IC)
         I=ILOW(ITR)
         J=IUP(ITR)
         K=NNEXT(IEL(I))
         IFR=INFR(IC)
         IJE=IJEX(IFR)
         FR=FREQ(IFR)
         FR15=FR*1.D-15
         BNU=BN*FR15*FR15*FR15
         WRITE(16,600) ITR,IFR,FR,2.997925d18/fr
         DO ID=1,ND
            T=TEMP(ID)
            TK=BOLK*T
            SB=2.0706D-16/T/SQRT(T)*G(I)/G(K)*EXP(ENION(I)/TK)
            IF(J.LT.K)
     *      SJ=2.0706D-16/T/SQRT(T)*G(J)/G(K)*EXP(ENION(J)/TK)
            PI=POPUL(I,ID)
            X=EXP(-HK*FR/T)
            PLTE=SB*ELEC(ID)*POPUL(K,ID)
            BI=PI/PLTE
            BJ=1.
            IF(J.LT.K) BJ=POPUL(J,ID)/SJ/ELEC(ID)/POPUL(K,ID)
            PLANCK=BNU/(UN/X-UN)
            IF(LINE(ITR)) THEN
               RD=RRD(ITR,ID)*G(I)/G(J)
               GG=G(I)/G(J)*POPUL(J,ID)
             ELSE
               GG=PLTE*X
               RD=RRD(ITR,ID)
            END IF
            SL=BNU*GG/(PI-GG)
            ggrad=0.
            heat=0.
            WRITE(16,601) ID,TAU(IC,ID),GGRAD,
     *                    BI,BJ,RRU(ITR,ID),RD,
     *                    RAD(IFR,ID),PLANCK,ST(IC,ID),SL,
     *                    HEAT
         END DO
      END DO
C
  600 FORMAT(/
     *     '  PARAMETERS FOR TRANSITION',I5,'  IFR =' ,I5,
     *     '  FREQ =',1PD15.5,'  Wavelength =',0pf11.3/
     *        1H ,8X,'TAU',6X,'GR ',6X,'B-I',6X,'B-J',6X,' RU',
     *        6X,' RD',6X,'RAD',5X,'PLANCK',5X,'STOT',6X,'SL',
     *        5X,'HEAT'/)
  601 FORMAT(1H ,I3,1P13D9.2)
  200 RETURN
      END