SUBROUTINE PFCNO(IAT,IZI,T,ANE,PF)
c     ===================================

c     Partition functions for the high ions of CNO:

c       (a) H-like ions:  C VI, N VII, O VIII
c       (b) He-like ions: N VI, O VII
c       (c) O VI from Sparks & Fischel, 1971, NASA SP-3066

c     Output:  PF   partition function

      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'

      PARAMETER (P1=0.1402,P2=0.1285,P3=1.,P4=3.15,P5=4.)
      DIMENSION TT(35),PN(10)
      DIMENSION P6A(24),P6B(10,11) 

      DATA TT /
     * 18.,19.,20.,21.,22.,23.,24.,25.,26.,27.,28.,29.,30.,
     * 32.,34.,36.,38.,40.,42.,44.,46.,48.,
     * 50.,55.,60.,65.,70.,75.,80.,85.,90.,95.,100.,125.,150./

      DATA PN /-2.,-1.,0.,1.,2.,3.,4.,5.,6.,7./
c
      DATA P6A /
     * 0.302, 0.302, 0.302, 0.303, 0.303, 0.304, 0.305, 0.306,
     * 0.307, 0.308, 0.310, 0.312, 0.313, 0.318, 0.322, 0.327,
     * 0.333, 0.339, 0.346, 0.353, 0.360, 0.367, 0.375, 0.394/
c
      DATA P6B /
     * 0.414,9*0.413,
     * 0.436,0.433,0.433,7*0.432,
     * 0.472,0.458,0.453,7*0.451,
     * 0.560,0.499,0.478,0.471,0.469,5*0.468,
     * 0.762,0.593,0.522,0.497,0.489,0.486,4*0.485,
     * 1.090,0.782,0.611,0.539,0.513,0.505,0.502,3*0.501,
     * 1.478,1.070,0.775,0.615,0.550,0.527,0.519,0.517,0.516,0.516,
     * 1.867,1.408,1.018,0.749,0.612,0.557,0.539,0.533,0.531,0.530,
     * 2.233,1.752,1.306,0.944,0.713,0.604,0.564,0.550,0.545,0.544,
     * 3.665,3.166,2.668,2.176,1.700,1.269,0.934,0.735,0.648,0.616,
     * 4.633,4.133,3.633,3.134,2.636,2.146,1.674,1.254,0.942,0.763/     
     
     
      IF (IAT.LT.6 .OR. IAT.GT.8 .OR. IZI.LE.5) THEN
         PF=0.
         PRINT *,' Routine PFCNO does not provide U for this ion '
     *         ,IAT,IZI
         STOP
      END IF
      IF (IZI.GT.IAT) THEN
         PF=1.
         RETURN
      END IF
     
      TK = BOLK*T
      IZIT = IAT-IZI
     
C    1. H-like case

      IF(IZIT.EQ.0) THEN
        ANEL = LOG(ANE)
        AA = 0.09*EXP(ANEL/6.)/SQRT(T)
        ANEL23 = EXP(-2./3.*ANEL)
        Z  = FLOAT(IZI)
        Z2 = Z*Z
        CBZ = 2.*8.59D14*Z*Z*Z
        E0KT = EH*Z2/TK
        
        PF=0.
        DO II=1,NLMX
          XN  = FLOAT(II)
          XN2 = XN*XN
          
          IF(XN.LE.3.01) THEN
            XKN = 1.
          ELSE
            XN1 = 1./(XN+1.)
            XKN = 16./3.*XN*XN1*XN1
          END IF
          BETA = CBZ*XKN/XN2/XN2*ANEL23
          X=EXP(P4*LOG(1.+P3*AA))
          C1=P1*(X+P5*(Z-1.)*AA*AA*AA)
          C2=P2*X
          F=(C1*BETA*BETA*BETA)/(1.+C2*BETA*SQRT(BETA))
          WI=F/(1.+F)
          
          EE = EXP(-E0KT/XN2)
          PF = PF+XN2*WI*EE
        END DO
        PF = 2.*PF
      END IF

C    2. He-like case

      IF(IZIT.EQ.1) THEN
        PF=1.
      END IF

C    3. O VI

      IF(IZIT.EQ.2) THEN
     
        IF(T.LT.18000.) THEN
           PF=2.
        ELSE
          NA = 24
          NB = 11
          PNE = LOG10(ANE*TK)
          T0  = 0.001*T
          J = 1
          IF(PNE.LT.PN(1)) GO TO 15
          IF(PNE.GT.PN(10)) THEN
            J1 = 10
            J2 = 10
            GO TO 16
          END IF
          DO J=1,9
            IF(PNE.GE.PN(J) .AND. PNE.LT.PN(J+1)) GO TO 15
          END DO
   15     J1 = J
          J2 = J+1
          IF(PNE.LT.PN(1)) J2 = 1
   16     DO I=1,34
            IF(T0.GE.TT(I) .AND. T0.LT.TT(I+1)) GO TO 25
          END DO
   25     I1 = I
          I2 = I+1
          IF(T0.GT.TT(35)) THEN
            I1 = 35
            I2 = 35
          END IF
          IF(I2.LE.24) THEN
            PX1=P6A(I1)
            PX2=P6A(I1)
            PY1=P6A(I2)
            PY2=P6A(I2)
          ELSE IF (I1.EQ.24) THEN
            PX1=P6A(I1)
            PX2=P6A(I1)
            PY1=P6B(J1,I2-24)
            PY2=P6B(J2,I2-24)
          ELSE
            PX1=P6B(J1,I1-24)
            PX2=P6B(J2,I1-24)
            PY1=P6B(J1,I2-24)
            PY2=P6B(J2,I2-24)
          END IF
          DLGUNX=PX2-PX1
          PX=PX1+(PNE-PN(J1))*DLGUNX
          DLGUNY=PY2-PY1
          PY=PY1+(PNE-PN(J1))*DLGUNY
          DELT=TT(I2)-TT(I1)
          IF(DELT.NE.0.) THEN
            DLGUT=(PY-PX)/DELT
            PF=PX+(T0-TT(I1))*DLGUT
          ELSE
            PF=PX
          ENDIF
          PF=EXP(2.302585093*PF)
        END IF
      END IF

      RETURN
      END