C./ ADD NAME=CHILC
C     PROGRAM TO CALCULATE, PRINT, AND PLOT TNE MAGNETIC                
C     SUSCEPTIBILITY OF HSP METALS.                                     
C                                                                       
      IMPLICIT REAL*8 (A-H,P-Z)
      REAL FUN(200,5),NAM(1000),NBOUND(10,2)
      DIMENSION S(1300,3),TM(2,2600),T(2),C(4),S1(3),TE(2,4),TO(2,5)
      COMMON/B/ CRIT, DE ,NTI                                       
      COMMON/A/NA, NS, KL, ISIG                                     
      COMMON/C/EPS, DLE                                             
      COMMON/D/ S12, S11                                            
      REAL*8 KL,KAPPA,KAPRM                                         
C                                                                   
      PI=3.1415927D0                                                
      open (21,FILE='FOR21.DAT')
      open (5,FILE='CHILC.OUT')
 7    READ (21,4100) KODE 
 400  FORMAT(' ENTER RESTART KODE ',I5,/)
      print   400, KODE                                                 
      IF (KODE.LT.1) KODE = 1                                           
      IF (KODE.GT.6) STOP                                               
      GO TO (2, 1, 5, 6, 22, 35), KODE                                  
 2    READ (21,600)  NTI, CRIT, DT, CN, I, NI, NR                       
 500  FORMAT (' ENTER NTI, CRIT, DT, CN, I, NI, NR',I5,3F10.5,3I5,/)    
      print   500, NTI, CRIT, DT, CN, I, NI, NR                         
 600  FORMAT (I5, 3F10.5, 3I5 )                                         
C     NTI IS NUMBER OF SIMPSON INTERVAL IN INTEGRATION FOR TM.          
C     CRIT IS WIDTH (IN KT) OF FERMI FUNCTION USED.                     
C     DT IS RANGE (IN KT) ABOVE FERMI LEVEL WHERE TM CALCULATED.        
C     CN CONTROLS SUBDIVIDING OF SIMPSON INTERVALS IN E INTEGRATION.    
C     IF CN = .0, ENTERED VALUES OF I, NI, NR ARE USEL.                 
C     I AND NI CONTROL RRNGE OF SUBDIVISION, FROM L=I TO L=I+2*NI.      
C     NR IS THE NUMWER OF SUBDIVISIONS PER INTERVAL, SDE=DE/(2NR).      
C     IF NTI = 0, STANDARD VALUES OF ABOVE ARE ENTERED, SEE BELOW.      
      IF (NTI.LT.0) GO TO 7                                             
      IF (NTI.GT.0) GO TO 1                                             
      NTI = 20                                                          
      CRIT = 10.                                                        
      DT = 10.                                                          
      CN = .0                                                           
      I = 0                                                             
      NI = 0                                                            
      NR = 1                                                            
C                                                                       
C     SET UP S VECTORS                                                  
 1    READ (21,2000)  NE, DELTA, EC, BET, ISIG                          
 1000 FORMAT (' ENTER NE, DELTA, EC, BETA, ISIG ',I5,3F10.5,I5,/)       
      print   1000, NE, DELTA, EC, BET, ISIG                            
C     DELTA, EC, AND BET=BETA=B/A ARE BENNETT-FALICOV BAND PARAMETERS.  
C     NE IS THE NUMBER OF SIMPSON INTERVALS, DE = (EMAX-EMIN)/(2NE).    
C     E2/(2DE) AND E3/(2DE) MUST BE INTEGERS, LATER IN THE              
C     PROGRAM, NE BECOMES 2NE+1.                                        
C     IN THIS SECTION IF ISIG GT 1, GET TRACE AND S1,S2,S;              
 2000 FORMAT (I5, 3F10.5, I5)                                           
      IF (NE.LE.0) GO TO 2                                              
      B2 = BET**.66666666D0                                             
      B4 = B2**2                                                        
      RP3 = (B4 + 2./B2)/3.D0                                           
      BE2 = BET**2                                                      
      E3 = EC + DELTA/3.D0                                              
      E2 = EC - DELTA/3.D0                                              
      EMIN = .0                                                         
      IF (E2.LT..0) EMIN = E2                                           
      EMAX = .0                                                         
      IF (E3.GT.0.) EMAX = E3                                           
      EMID = E2                                                         
      IF (EMID.EQ.EMIN) EMID = .0                                       
      IF (EMID.EQ.EMAX) EMID = E3                                       
      EL = .0                                                           
      IF (EC.EQ.0.) GO TO 4                                             
      EL=8.*EC**3/((((4.*RP3/(3.*B4))-1.)*EC**2+DELTA**2/9.D0)*27.*BE2) 
      IF (BET.EQ.1.) GO TO 4                                            
 21   KAPPA = EL*(EL-EC-DELTA/3.)*(EL-EC+DELTA/3.)                      
      KAPRM = 2.*EL*(EL-EC) + ((EL-EC)**2-(DELTA/3.)**2)                
      H3 = RP3*EL - 2.*EC/(3.*B2)                                       
      DEL = -(KAPPA-H3**3)/(KAPRM-3.*RP3*H3**2)                         
      EL = EL + DEL                                                     
      IF (DABS(DEL/(EL-EMID)).GT..001) GO TO 21                         
 4    print  1500, EL, E2, E3                                           
 1500 FORMAT (' EL =',E12.4,' E2 = ',E12.4,' E3 = ',E12.4,/ )           
      DE = (EMAX - EMIN)/(2.*NE)                                        
      NE = 2*NE + 1                                                     
 9    IF (NE.GT.1300) GO TO 90                                          
      IT=0                                                              
      NIT=0                                                             
      NB=0                                                              
      IF (CN.EQ.0.) GO TO 101                                           
C     PICK SUBDIVISION PARAMETERS                                       
      I=0                                                               
      NI=0                                                              
      NR=1                                                              
 101  C(1)=(EMAX-EMID)/DE                                               
      C(2)=(EMID-EMIN)/DE                                               
      C(3)=(EL-EMID)/DE                                                 
      TEMP=DABS(C(3))                                                   
      IF (C(1).GT.C(2)) GOTO 104                                        
      IF (C(1).LT.1.) GOTO 11                                           
      IF (C(1).GE.CN) GOTO 104                                          
      NR = CN/C(1)+1.                                                   
      NI = C(1)+.1                                                      
      I = NE-2*NI                                                       
      IF (I.GE.1) GOTO 107                                              
      NI = NI+(I-1)/2                                                   
      I = 1                                                             
      GOTO 107                                                          
 11   IF (TEMP.LE.2.) GOTO 95                                           
      EMAX=EMAX - 2.*DE                                                 
      NE = NE-2                                                         
      IF (CN.EQ.0.) GOTO 109                                            
      C(1) = C(1)-2.                                                    
      C(3) = C(3)+2.                                                    
      TEMP = TEMP-2.                                                    
      I = NE-2                                                          
      NI=TEMP/4.                                                        
      IF (NI.LT.2) NI=2                                                 
      NR = CN                                                           
      GOTO 107                                                          
 104  IF (C(2).LT.1.) GOTO 12                                           
      IF (C(2).GT.CN) GOTO 107                                          
      NR = CN/C(2)+1.                                                   
      NI = C(2)+.1                                                      
      I = 1                                                             
      GOTO 107                                                          
 12   IF (TEMP.LE.2.) GOTO 95                                           
      EMIN = EMIN+2.*DE                                                 
      NE = NE-2                                                         
      IF (CN.EQ.0.) GOTO 109                                            
      C(2) = C(2)-2.                                                    
      C(3) = C(3)-2.                                                    
      TEMP = TEMP-2.                                                    
      I = 1                                                             
      NI = TEMP/4.                                                      
      IF (NI.LT.2) NI=2                                                 
      NR = CN                                                           
 107  IF (TEMP.GE.CN.OR.CN.EQ.0.) GOTO 109                              
      IF (TEMP.LT..5/CN) GOTO 109                                       
      NRT = CN/TEMP+.999                                                
      NIT = IFIX(SNGL(TEMP+1.))+1                                       
      IT = -1-2*(IFIX(SNGL(TEMP))/4)                                    
      IF (C(3).LT.0.) IT = 2-2*NIT-IT                                   
      IT = IFIX(SNGL(C(2)+.1))+IT                                       
      IF (NRT.GT.NR) NR = NRT                                           
      IF (IT.LT.1) IT = 1                                               
      IF (I.EQ.0) I = IT                                                
      IF (I.LE.IT) GOTO 102                                             
      NI = NI+(I-IT)/2                                                  
      I = IT                                                            
 102  IF (I+2*NI.LT.IT+2*NIT) NI = (IT-I)/2+NIT                         
 109  print  1501, I, NI, NR, NB                                        
C                                                                       
C     ENERGY INTEGRATION                                                
      E = EMIN                                                          
      DO 111 L = 1,NE                                                   
      DO 111 K = 1,3                                                    
 111  S(L,K)=.0                                                         
      IF (ISIG.GT.1) print  1501, I, NI, NR, NB                         
 1501 FORMAT(' I, NI, NR, NB = ',4I4,/)                                 
      NB = 1                                                            
      IF (I.EQ.1)  NB = NR                                              
      IF (ISIG.GT.1) print  1501, I, NI, NR, NB                         
      SDE = DE/NB                                                       
      DLE = .0                                                          
      CALL SS( E, SDE, S1, DELTA, EC, BET)                              
      NEM1= NE-1                                                        
      DO 3 J = 2, NEM1, 2                                               
      NB = 1                                                            
      IF (J.GT.I.AND.J.LT.I+2*NI) NB = NR                               
      IF (ISIG.GT.1) print  1502, I, NI, NR, NB, J                      
 1502 FORMAT (5I4)                                                      
      SDE = DE/NB                                                       
      DO 112 K = 1,3                                                    
 112  S(J-1,K) = S(J-1,K)+S1(K)/(2.*NB)                                 
      IF (ISIG.GT.1) print  6000, E, ( S(J-1,K),K = 1,3)                
      NOSC = 1                                                          
      A = -1.+1./NB                                                     
      NB2 = NB*2                                                        
      DO 115 M = 1,NB2                                                  
      E = E+SDE                                                         
C     IF (J.EQ.NE-1.AND.M.EQ.2*NB) E = EMAX                             
C     IF (DABS(E-EMID).LT..5*SDE.AND.M.EQ.2*NB) E = EMID                
C     DLE = .0                                                          
      EPS = E-EL                                                        
      IF (EPS.GT.-.5*SDE.AND.EPS.LE..5*SDE) DLE = SDE                   
      IF (DLE.EQ.0.D0) GOTO 114                                         
      DO 113 K = 1,3                                                    
 113  TM(1,K) = S1(K)                                                   
      IF (DLE.NE.0.D0.AND.DABS(E-EMID).LT..5*SDE) GOTO 17               
 114  CALL SS( E, SDE, S1, DELTA, EC, BET)                              
      IF (DLE.EQ.0.D0) GOTO 15                                          
C                                                                       
C     CALCULATE CONTRIBUTIONNOF SINGULARITY AT E = EL                   
      DO 13 K = 1,3                                                     
 13   C(K) = S1(K)                                                      
      DLE = .0D0                                                        
      CALL SS( E+SDE, SDE, S1, DELTA, EC, BET)                          
      IF (DABS(E-EMID).GT..1*SDE) GOTO 119                              
      S1(1) = S1(1)+S11                                                 
      S1(2) = S1(2)+S12                                                 
 119  TEMP = EPS/SDE                                                    
      DO 14 K = 1,3                                                     
      IF (MOD(M,2).NE.0) S1(K)=(2.*S1(K)+2.*TM(1,K)+(3.D0*TEMP          
     1 *DLOG((1.+TEMP)/(1.-TEMP))-6.D0)*C(K))/4.D0                      
      IF (MOD(M,2).EQ.0) S1(K) = (S1(K)+TM(1,K)+(3.D0*TEMP              
     1 *DLOG((1.+TEMP)/(1.-TEMP))+.5D0*DLOG((4.-TEMP**2)/(1.-TEMP**2)   
     2 )-6.D0)*C(K))/2.                                                 
 14   CONTINUE                                                          
      GOTO 15                                                           
C                                                                       
C     CALCULATE CONTRIBUTION OF SINGULARITY AT E = EL = EMID            
 17   DLE = .0D0                                                        
      DTE = DABS(EMID-2.*EC/(3.*RP3*B2))                                
      IF (DTE.GT..5*SDE) DTE = .5*SDE                                   
      IF (DTE.LT..05*SDE) DTE = .05*SDE                                 
      TW3 = 2.**.33333333D0                                             
      TEMP = 3.D0/(2.D0*TW3)-1.                                         
      TEMP3 = (SDE/DTE)**.33333333D0-2.+DTE/SDE+(1.-DTE/SDE)/TW3        
      TEMP1 = 3.D0*TEMP/TEMP3                                           
      TEMP2 = (1.+3.D0*TEMP*(DTE/SDE-2.)/TEMP3)                         
      TEMP3 = -.5D0+3.*TEMP*(1.-DTE/SDE)/TEMP3                          
      CALL SS ( E-2.*SDE, SDE, C, DELTA, EC, BET)                       
      DO 117 K = 1,3                                                    
 117  TM(1,K) = TEMP2*S1(K)+TEMP3*C(K)                                  
      CALL SS ( E+SDE, SDE, S1, DELTA, EC, BET)                         
      CALL SS ( E+2.*SDE, SDE, C, DELTA, EC, BET)                       
      DO 121 K = 1,3                                                    
 121  TM(1,K) = TM(1,K)+TEMP2*S1(K)+TEMP3*C(K)                          
      CALL SS ( E-DTE, DTE, S1, DELTA, EC, BET)                         
      CALL SS ( E+DTE, DTE, C, DELTA, EC, BET)                          
      DO 118 K = 1,3                                                    
 118  S1(K) = TM(1,K)+TEMP1*(S1(K)+C(K))                                
C                                                                       
 15   TEMP = (3.D0+NOSC)/(4.D0*NB)                                      
      IF ( M.EQ.2*NB) TEMP = TEMP/2.                                    
      IF (ISIG.GT.1) print  1503, M, NOSC, A, TEMP                      
 1503 FORMAT (' M, NOSC, A, TEMP ', 2I4, 2E12.4)                        
      DO 116 K = 1,3                                                    
      S(J-1,K) = S(J-1,K)-TEMP*A*(1.-A)*S1(K)                           
      S(J,K) = S(J,K)+TEMP*(1.-A**2)*S1(K)                              
 116  S(J+1,K) = S(J+1,K)+TEMP*A*(1.+A)*S1(K)                           
      IF (ISIG.GT.1) print  6000, E, ( S(J-1,K), K = 1,3)               
      IF (ISIG.GT.1) print  6000, E, ( S(J,K), K = 1,3)                 
      IF (ISIG.GT.1) print  6000, E, ( S(J+1,K), K = 1,3)               
      A = A+1./NB                                                       
  115 NOSC = -NOSC                                                      
      IF (DABS(E-EMID).GT..5*SDE.OR.DABS(EPS).LT..5*SDE) GOTO 3         
C                                                                       
C     STEP FUNCTION CORRECTOR                                           
      TEMP = S11*DSIGN(1.D0,EL-EMID)                                    
      S(J+1,1) = S(J+1,1)-TEMP/(2.*NB)                                  
      S1(1) = S1(1)+TEMP                                                
      TEMP = S12*DSIGN(1.D0,EL-EMID)                                    
      S(J+1,2) = S(J+1,2)-TEMP/(2.*NB)                                  
      S1(2) = S1(2)+TEMP                                                
C                                                                       
 3    CONTINUE                                                          
      DO 120 K = 1,3                                                    
      S(1,K) = 2.*S(1,K)                                                
 120  S(NE,K) = 2.*S(NE,K)                                              
C                                                                       
C     PLOT OR PRINT S; ISIG = 0, PLOT S(M); ISIG > 0. PRINT;            
C     ISIG < 0. PRINT AND PLOT.                                         
      READ (21,4100)  ISIG,M                                            
 3000 FORMAT (' ENTER ISIG,M ',2I5,/)                                   
      print   3000, ISIG,M                                              
 4100 FORMAT ( 2I5 )                                                    
      IF (ISIG) 5,5,6
 5    IF (M.LT.1.OR.M.GT.3) GOTO 20                                     
      SMIN = S(1,M)                                                     
      SMAX = SMIN                                                       
      DO 16 L = 2,NE                                                    
      SMIN = DMIN1(SMIN,S(L,M))                                         
 16   SMAX = DMAX1(SMAX,S(L,M))                                         
      print  4200, SMIN, SMAX                                           
 4200 FORMAT (' SMIN =',E12.4,' SMAX =',E12.4/' ENTER SMIN,SMAX '/)     
C     READ (21,4000), SMIN, SMAX                                        
 4000 FORMAT (2F10.5)                                                   
C     CALL PLOT (EMIN,.0,12)                                            
C     CALL PLOT (.0,SMIN,13)                                            
C     CALL PLOT (.0,SMIN,12)                                            
C     CALL PLOT (.0,SMAX,12)                                            
C     CALL TTY                                                          
 6    IF (ISIG.EQ.3) print  5000                                        
 5000 FORMAT(' ',5X,'E',8X,'S(1)',8X,'S(2)',8X,'N(E)')                  
      E = EMIN                                                          
      DO 10 L = 1,NE                                                    
      IF (ISIG.EQ.3) print  6000, E, ( S(L,K), K = 1,3)                 
 6000 FORMAT (F12.6, 3E12.4)                                            
C     IF (ISIG.LE.0.AND.L.EQ.1) CALL PLOT(EMIN,S(L,M),13)               
C     IF (ISIG.LE.0) CALL PLOT ( E, S(L,M), 12)                         
C     IF (ISIG.LT.0) CALL TTY                                           
 10   E = E+DE                                                          
C     CALL TTY                                                          
      IF (ISIG.LE.0) GOTO 5                                             
C                                                                       
C     PREPARE S WITH SIMPSON FACTORS                                    
 20   OSC = DE/3.D0                                                     
      DO 8 L = 1,NE                                                     
      OSC = -OSC                                                        
      SIMP = DE+OSC                                                     
      IF (L.EQ.1.OR.L.EQ.NE) SIMP = DE/3.D0                             
      DO 8 K = 1,3                                                      
 8    S(L,K) = S(L,K)*SIMP                                              
C                                                                       
      LGR = 1                                                           
C     SET UP TM'S                                                       
 22   READ (21,8000)  EFMIN, EFMAX, TRY, A, B                           
 7000 FORMAT (' ENTER  EFMIN, EFMAX, TRY, ABAR, B ',5F10.6,/)           
      print   7000, EFMIN, EFMAX, TRY, A, B                             
 8000 FORMAT (5F10.6)                                                   
C     EFMIN, EFMAX ARE MINIMUM AND MAXIMUM FERMI LEVELS                 
C     TRY IS THE TEMPERATURE (IN ENERGY UNITS)                          
C     ABAR AND B ARE BAND PARAMETERS (NOT THE B/A =BETA)                
      IF (B.EQ.0.) GOTO 1                                               
C     FMUL = 1.71732D-6*A**2/DSQRT(B)                                   
      FMUL = 1.71732D-6*A**2/DSQRT(B)   /7.14                           
      FNE = 1.8769D4*RP3/A**4                                           
      NM = -(EFMIN-EMIN)/DE                                             
      EFMIN = EMIN-NM*DE                                                
      NM = (EFMAX-EMIN)/DE                                              
      EFMAX = EMIN+NM*DE                                                
      NF = (EFMAX-EFMIN)/DE+1.1                                         
      NT = (EFMAX-EMIN+DT*TRY)/DE+2.1                                   
      IF (NT.GT.NE+NF) NT = NE+NF                                       
      IF (NT.GT.2600) GOTO 93                                           
      IF (NT.LT.1) GOTO 35                                              
      EMU = EMIN-EFMAX                                                  
      DO 30 L = 1,NT                                                    
      CALL TN( 1.D0, TRY, EMU, T)                                       
      IF (ISIG.EQ.4) print  6000, EMU, ( T(K), K = 1,2)                 
      DO 25 K = 1,2                                                     
 25   TM(K,L) = T(K)*FMUL                                               
 30   EMU = EMU+DE                                                      
      IF (TRY.NE.0.) GOTO 35                                            
      CALL TONE(2.D0,TO)                                                  
      SM = DSQRT(DE)*FMUL                                               
      DO 32 K = 1,2                                                     
      DO 31 N = 1,4                                                     
 31   TO(K,N) = TO(K,N)*SM                                              
      IF (ISIG.GT.2) print  6000, SM, (TO(K,L), L=1,4)                  
 32   SM = SM/DE                                                        
      OSC = 1.D0/3.D0                                                   
      DO 34 N = 1,4                                                     
      DO 33 K = 1,2                                                     
      TE(K,N) = TO(K,N)/(1.-OSC)                                        
 33   TO(K,N) = TO(K,N)/(1.+OSC)                                        
 34   OSC = -OSC                                                        
      DO 28 K = 1,2                                                     
      TO(K,5) = TO(K,4)                                                 
      TO(K,4) = TO(K,3)                                                 
      TO(K,3) = TO(K,2)                                                 
      TO(K,2) = TO(K,1)+3.D0*TM(K,NT-3)/8.D0                            
      TO(K,1) = 5.D0*TM(K,NT-4)/4.D0                                    
 28   IF (ISIG.GT.2) print  6000, OSC, (TO(K,L), L = 1,5)               
      DO 29 K = 1,2                                                     
      TM(K,NT-3) = TE(K,1)+TM(K,NT-3)/2.                                
      TM(K,NT-2) = TE(K,2)                                              
      TM(K,NT-1) = TE(K,3)                                              
 29   TM(K,NT) = TE(K,4)                                                
C                                                                       
C     PLOT OR PRINT C'S AND CHI = C(3), N(E) = C(4)                     
C     IF ISIG = 0, PLOT C(M)                                            
C     IF ISIG > 0, PRINT ALL C'S                                        
C     IF ISIG < 0, PRINT ALL C'S AND PLOT C(M)                          
 35   READ (21,8800)  ISIG, M,NPR                                       
 8500 FORMAT(' ENTER ISIG, M, NPR',3I5,/)                               
      print   8500, ISIG, M, NPR                                        
 8800 FORMAT (3I5)                                                      
      IF (M.LE.0.OR.M.GT.4) GOTO 22                                     
      IF (ISIG) 36,36,37                                                
 36   print  9000                                                       
 9000 FORMAT (' ENTER CMIN, CMAX '/)                                    
      READ (21,4000)  CMIN, CMAX                                        
      DO 52 K=1,5                                                       
      NBOUND(K,1) = CMIN                                                
 52   NBOUND(K,2) = CMAX                                                
 9700 FORMAT(9X,'EF',4X,'C(1)',7X,'C(2)',7X,'CHI',8X,'N(E)',/)          
 37   IF (ISIG.GE.1) print  9700                                        
C                                                                       
C     COMRUTE C'S                                                       
      EF = EFMIN                                                        
      DO 50 N = 1,NF                                                    
      DO 38 K = 1,4                                                     
 38   C(K) = .0D0                                                       
C                                                                       
C     END CORRECTIONS                                                   
      NTI = 2*NTI                                                       
      CRIT = 2.*CRIT                                                    
      DEN = .0D0                                                        
      IF (DABS(E2).LT.DE) DEN = DE                                      
      CALL SS ( EMIN, DEN, S1, DELTA, EC, BET)                          
      IF (DABS(E2).LT.DE) S12 = S1(2)                                   
      CALL TN(1.D0, TRY, EMIN-EF, T)                                    
      C(1) = -FMUL*T(1)*S12                                             
      DEN = .0D0                                                        
      IF (DABS(E3).LT.DE) DEN = DE                                      
      CALL SS ( EMAX, DEN, S1, DELTA, EC, BET)                          
      IF (DABS(E3).LT.DE) S12 = S1(2)                                   
      CALL TN( 1.D0, TRY, EMAX-EF, T)                                   
      C(1) = C(1)+FMUL*T(1)*S12                                         
      IF (DABS(E2).GT.DE.AND.DABS(E3).GT.DE) GOTO 45                    
C                                                                       
C     END CORRECTIONS IF E2 OR E3 = 0                                   
      CALL TN( 1.D0, TRY, -EF, T)                                       
      C(2) = C(2)-PI*FMUL*T(2)/(3.*RP3*B2)                              
      C(1) = C(1)+FMUL*T(1)*PI*(2.*B4/3.D0+.5D0*(RP3+1./B2))/(RP3*EC)   
      TEMP = -2.D0*FMUL*PI*B4*DE/(3.D0*RP3*EC*DABS(EC))                 
      C(1) = C(1)+TEMP*T(1)                                             
      CALL TN (1.D0, TRY, DSIGN(DE,EC)-EF, T)                           
      C(1) = C(1)+4.D0*TEMP*T(1)                                        
      CALL TN( 1.D0, TRY, 2.*DSIGN(DE,EC)-EF, T)                        
      C(1) = C(1)+TEMP*T(1)                                             
 45   NTI = NTI/2                                                       
      CRIT = CRIT/2.D0                                                  
      NL = NT+N-NF                                                      
      IF (NL.LT.3) GOTO 41                                              
      IF (NL.GT.NE) NL = NE                                             
      DO 40 L = 1,NL                                                    
      NTL = NF-N+L                                                      
      IF (NTL.GT.NT) GOTO 41                                            
      TMM = TM(2,NTL)                                                   
      IF(TRY.NE.0) GOTO 39                                              
      IF (MOD(NL,2).NE.0.AND.NTL.GE.NT-4) TMM = TO(2,NTL-NT+5)          
 39   C(4) = C(4)+S(L,3)*TMM                                            
      DO 40 K = 1,2                                                     
      TMM = TM(K,NTL)                                                   
      IF (TRY.NE.0) GOTO 40                                             
      IF (MOD(NL,2).NE.0.AND.NTL.GE.NT-4) TMM = TO(K,NTL-NT+5)          
 40   C(K) = C(K)+S(L,K)*TMM                                            
 41   DO 42 K = 1,2                                                     
 42   C(3) = C(3)+C(K)                                                  
      C(4) = FNE*C(4)                                                   
      NN=N/NPR
      IF (ISIG.GE.1.AND.NN*NPR.EQ.N) print  9800, EF, ( C(K), K = 1,4)  
 9800 FORMAT (F12.8,4D11.3)                                             
      IF (NN*NPR.EQ.N) NAM((LGR-1)*(NF/NPR)+NN) = C(M)                  
      IF (LGR.NE.1.OR.N.NE.1) GO TO 48                                  
      CMIN = C(M)                                                       
      CMAX = C(M)                                                       
 48   CMIN = DMIN1(CMIN, C(M))                                          
      CMAX = DMAX1(CMAX, C(M))                                          
 49   CONTINUE                                                          
 50   EF = EF+DE                                                        
C     DO 52 K=1,LGR                                                     
C     NBOUND(K,1) = CMIN                                                
C52   NBOUND(K,2) = CMAX                                                
      LGR=LGR+1                                                         
      IF (ISIG.GE.0) GO TO 22                                           
      LGR=LGR-1                                                         
      CALL GRAF (EFMIN,EFMAX,NAM,NBOUND,NF/NPR,LGR  ,FUN)               
      print  9900, CMIN, CMAX                                           
 9900 FORMAT (' CMIN = ', E12.4,' CMAX = ', E12.4)                      
C                                                                       
      IF (ISIG.EQ.2) GO TO 1                                            
      IF (ISIG.EQ.-1) STOP                                              
C                                                                       
C     ERROR MESSAGES                                                    
   90 print  1100                                                       
 1100 FORMAT (' NE TOO LARGE')                                          
      GOTO 1                                                            
 93   print  1200                                                       
 1200 FORMAT(' NT TOO LARGE')                                           
      GOTO 22                                                           
 95   print  1300                                                       
 1300 FORMAT (' NE TOO SMALL')                                          
      GOTO 1                                                            
      END                                                               
C                                                                       
      SUBROUTINE TN(B, TRY, EMU, T )                                    
C                                                                       
      IMPLICIT  REAL *8 (A-H,O-Z)                                       
C     SUBROUTINE TO CALCULATE THE INTEGRAL OF THE FERMI FUNCTION        
C     OF E = B*Z**2, T(1), AND ITS DERIVATIVE WITH RESPECT TO           
C     FERMI LEVEL, T(2), TRY IS THE TEMPERATURE (IN ENERGY UNITS)       
C     AND EMU = E - MU,WHERE MU IS THE FERMI LEVEL,                     
      DIMENSION T(2)                                                    
      COMMON/B/     CRIT, DE ,NT                                        
C                                                                       
      PI = 3.1415927D0                                                  
      IF (DABS(EMU).LT..1*DE) EMU = .0D0                                
      T(1) = .0D0                                                       
      T(2) = .0D0                                                       
      N = NT                                                            
      IF (TRY.EQ.0.) GO TO 40                                           
      BETA = 1.D0/TRY                                                   
      BEMU = BETA*EMU                                                   
      IF (BEMU.GT.CRIT) GO TO 70                                        
      IF (BEMU.LT.-CRIT) GO TO 80                                       
      Z0 = 0.D0                                                         
      ARF = ((-CRIT)/BETA-EMU)/B                                        
C     IF (ARF.GT.0D0)  Z0 = DSQRT(ARF)                                  
      ZL = DSQRT((((CRIT+5.D0)/BETA) -EMU)/B)                           
C                                                                       
      DZ = (ZL-Z0) /(2*N )                                              
      OSC = DZ/3.D0                                                     
      Z = Z0                                                            
      N2P1  = 2*N  + 1                                                  
      DO 50 L = 1, N2P1                                                 
      OSC = -OSC                                                        
      EX = .0                                                           
      SIMP = DZ + OSC                                                   
      IF (L.EQ.1.OR.L.EQ.2*N +1) SIMP = DZ/3.D0                         
      BZ2 = B*Z*Z                                                       
      ARG = BETA * ( EMU + BZ2 )                                        
      IF (DABS(ARG).GT.4.*CRIT) GO TO 50                                
      EX = DEXP(ARG)                                                    
      D = 1.D0 + EX                                                     
      DFDM = BETA*EX/D**2                                               
      T(1) = T(1) + 2.D0*SIMP*BZ2*DFDM                                  
      T(2) = T(2) + SIMP*DFDM                                           
 50   Z = Z + DZ                                                        
      IF (EX.EQ.0D0) RETURN                                             
      TEMP = .5/(EX*BETA*B*(Z-DZ))                                      
      DO 60 L = 1, 2                                                    
      T(L) = T(L) + TEMP                                                
 60   TEMP = BETA*TEMP                                                  
      RETURN                                                            
C                                                                       
 40   IF (EMU.GE.0.D0) RETURN                                           
      T(1) = DSQRT(-EMU/B)                                              
      T(2) = -.5D0*T(1)/EMU                                             
      RETURN                                                            
C                                                                       
 70   IF (BEMU.GT.50) RETURN                                            
      EX = DEXP(-BEMU)                                                  
      TEMP = .5D0*EX*DSQRT(PI/(BETA*B))                                 
      DO 75 L = 1, 2                                                    
      T(L) = TEMP                                                       
 75   TEMP = TEMP*BETA                                                  
      RETURN                                                            
C                                                                       
C                                                                       
 80   T(1) =DSQRT(-EMU/B)                                               
      T(2) = -.5D0*T(1)/EMU                                             
      TEMP = (.5D0*T(2)/EMU)*((PI/BETA)**2)/6.D0                        
      T(1) = T(1) + TEMP                                                
      T(2) = T(2) + 1.5*TEMP/EMU                                        
      RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
      SUBROUTINE TONE(X,T)                                              
      IMPLICIT  REAL *8 (A-H,O-Z)                                       
      DIMENSION T(2,4)                                                  
      X2 = X*X                                                          
      R = DSQRT(X)                                                      
      T(1,1) = X2*R*(X2/9.D0-.2D0 )/3.D0                                
      T(2,1) = X*R*(X2/6.D0-X2/7.D0+1.D0/9.D0-1.D0/6.D0)                
      T(1,2) = X2*R*(-X2/9.D0+X/7.D0+.4D0 )                             
      T(2,2) = X*R*(3.D0*X2/7.D0-2.D0*X/5.D0-2.D0/3.D0-(X+1.)*(X-2.)/2.)
      T(1,3) = X*R*(X*X2/9.D0-2.D0*X2/7.D0-X/5.D0+2.D0/3.D0)            
      T(2,3) = -X*R*(3.D0*X2/7.D0-4.D0*X/5.D0-1.D0/3.D0)                
      T(1,4) = -X2*R*(X2/9.D0-3.D0*X/7.D0+2.D0/5.D0)/3.D0               
      T(2,4) = X*R*(3.D0*X2/7.D0-6.D0*X/5.D0+2.D0/3.D0)/3.D0            
      RETURN                                                            
      END                                                               
      SUBROUTINE GRAF (EFMIN,EFMAX,NAM,NBOUND,NF,LGR,FUN)
      LOGICAL*1 FIELD,EQSMBL,SMBLS(5)
      LOGICAL LEQ,LALLX
      INTEGER GR,DIGY,GRINTX,FSTNX
      REAL *8 LBX,STX,LBY,SCY,EFMIN,EFMAX
      REAL NAM(310),NBOUND(10,2),FUN(NF,LGR)
      DATA FIELD/1H /,EQSMBL/1H*/,SMBLS/1H*,1H+,1H=,1H^,1H-/
      LALLX=.FALSE.
      LEQ=.FALSE.
      DO 22 L=1,LGR
      DO 22 N=1,NF
 22   FUN(N,L)=NAM((L-1)*NF+N)
      LBX=EFMIN
      STX=(EFMAX-EFMIN)/(NF-1)
      LBY=NBOUND(1,1)
      SCY=NBOUND(1,2)
      CALL MAZIL2(FIELD,EQSMBL,LEQ,LGR,NF,FUN,SMBLS,LBX,STX,
     *            LBY,SCY,3,3,10,1,LALLX,0)
      RETURN
      END
