C
      SUBROUTINE MOLFYS(T,XKH2,XKH2P,DEH2,DEH2P,deh2nodis,deh2pnodis)
C
C        THIS ROUTINE GIVES DISSOCIATION CONSTANTS XKH2 (=N(H I)*N(H I)/(NH2))
C        AND XKH2P (=N(H I)*N(H II)/N(H2+)), EXPRESSED IN NUMBER PER CM3 AND
C        THE SUM OF ROTATION AND VIBRATION ENERGIES, DEH2 AND
C        DEH2P FOR H2 AND H2+, RESPECTIVELY (EXPRESSED IN ERGS PER MOLECULE)
C        THE DATA ARE FROM VARDYA, M.N.R.A.S. 129, 205 (1965) AND EARLIER
C        REFERENCES. THE DISSOCIATION CONSTANT FOR H2 IS FROM TSUJI,
C        ASTRON. ASTROPHYS. 1973.
C
      DIMENSION A1(5),A2(4),B1(5),B2(5),TE(5)
      DATA A1/12.739,-5.1172,1.2572E-1,-1.4149E-2,6.3021E-4/,
     *A2/11.20699 ,-2.794276 ,-0.079196   ,0.024790   /,
     *B1/2.6757,-1.4772,0.60602,-0.12427,0.009750/,
     *B2/2.9216,-2.0036,1.7231,-0.82685,0.15253/
      TEX=5040./T
      TE(1)=1.
      DO1 K=1,4
    1 TE(K+1)=TE(K)*TEX
      XKH2=0.
      XKH2P=0.
      DEH2=0.
      DEH2P=0.
      DO2 K=1,4
      XKH2=A1(K)*TE(K)+XKH2
      XKH2P=A2(K)*TE(K)+XKH2P
      DEH2=B1(K)*TE(K)+DEH2
    2 DEH2P=B2(K)*TE(K)+DEH2P
      XKH2=A1(5)*TE(5)+XKH2
      DEH2nodis=(B1(5)*TE(5)+DEH2)*8.617E-5*T
      DEH2Pnodis=(B2(5)*TE(5)+DEH2P)*8.617E-5*T
      DEH2=deh2nodis-4.476
      DEH2P=deh2pnodis*T-2.648
      XKH2=10.**XKH2
      XKH2P=10.**XKH2P
      RETURN
      END