C
C                    capote_micro.for v2.0
C
C
C     MICROSCOPICAL PARTICLE-HOLE STATE DENSITIES CALCULATION (May 1997)
C
C     Contact person: Dr.Roberto Capote
C     e-mail: rcapote@infomed.sld.cu
C     e-mail: rcapote@ceaden.edu.cu
C     FAX: (537)221518
C
C  This code is an improved version (almost completely rewritten) 
C  of the original ICAR code written by M.HERMAN and G.REFFO (1986-1988).
C  It was developed by R.CAPOTE and R.PEDROSA in 1995-1996(BOLOGNA-HAVANA).
C
C  REFERENCES:
C  1) Herman M., Reffo G., Phys. Rev.C36(1987) 1546
C  2) Herman M., Reffo G.,	Comp.Phys.Comm.47(1987) 103.
C  3) Herman M., Reffo G., "Combinatorial Approach to few quasiparticle
C     state densities", Proc.Worshop on Applied Nuclear Theory and
C     Nuclear Model Calculations for Nuclear Technology Applications,
C     15 Feb.-19 Mar.1988, ICTP Trieste, Italia, Edited by M.K.Mechta
C     and J.J.Schmidt, World Scientific, Singapore 1989.
C  4) Herman M., Reffo G., Rego R.A., Phys.Rev.C37(1988)797
C  5) Capote R., Pedrosa R., DENSIDAD code(1996), unpublished.
C
C  COMMENT:
C  Using Moller single-particle levels(FILE:moeller_levels.dat) you
C  will be able to calculate usually only STATE densities as the
C  corresponding single-particle level schemes are calculated for some
C  deformation. Therefore, using Moller single-particle levels it is
C  impossible to get level density RHO by the method employed in the code
C       (level density RHO(E,J)=OMEGA(M=J)-OMEGA(M=J+1))
C
C  R.Capote can provide single-particle levels at zero deformation for 
C  spherical nucleus on request.
C
C  OPERATING ENVIROMENT: LINUX workstation
C
C  The code was developed and tested on PC computer running LINUX.
C  with the F77 compiler.
C
C--------------------------------------------------------------------------
C  INPUT INFORMATION:
C
C  It is designed for interactive input, which is self explained.
C  You should prepare also some input files as well as files, containing
C  shell model single-particle levels. The single-particle levels files
C  are provided by the code moller_levels.for(included in the Supplement)
C
C  For a given nucleus XXXX all corresponding input files of the type:
C   nXXXX11.inp,nXXXX22.inp,nXXXX33.inp  are used if present.
C   pXXXX11.inp,pXXXX22.inp,pXXXX33.inp
C   (XXXX=fe58 in the test example).
C--------------------------------------------------------------------------
C
C  The main difference with original code is related to the particle-hole
C  state densities at low excitation energies. The original code do not
C  give a same value of the w(p,h,E) for the spherical and sligthly deformed
C  nucleus. But, from the physical point of view, you should get the same
C  results in both cases with small redistribution of the levels.
C  Differences were related to different treatment of the excited particles
C  and holes in the subshell corresponding to the Fermi level.
C
C  The code was tested up to 3p-3h configurations, checking by hand that
C  all possible configurations were generated.
C-----------------------------------------------------------------------------
C  SAMPLE INPUT FILE: pfe5811.inp
Cline#1:  P                     Proton || Neutron DENSITY CALCULATIONS
Cline#2:  1,1                   Particle/Hole Numbers
Cline#3:  58,26                 A,Z of the nucleus
Cline#4:  100,0.D0              Number of SPL orbitals to be considered,Bn
Cline#5:  1,4.D0,1              Calc.Key(0,1,2) ,Eog(ESM spacing),DEBUG key(0/1)
Cline#6:  13.85D0               G*A
Cline#7:  20                    Maximum excitation energy to be considered
Cline#8:  z26a58.spe            Shell Model input spectrum filename
C-----------------------------------------------------------------------------
C  COMMENTS: 1) The shell model input spectra file is provided by
C      moller_levels.for code, which extract single-particle levels from the
C      Moeller-Nix single-particle levels database(FILE:moller_levels.dat)
C            2) Bn = Binding energy is used only if you want to get bounded
C      p-h state densities(for MSC calc.for example)
C            3)  Calc.key is as follows:
C            0 : Equidistant Spacing Model with constant spacing Eog from the
C                input file is assumed
C            1 : SHELL MODEL SPECTRUM INPUT(obtained by moller_levels.for)
C            2 : Equidistant Spacing Model with constant spacing Eog from the
C                SHELL MODEL SPECTRUM
C ----------------------------------------------------------------------------
C   Using this input file you will get the following output files(LD=level density):
C pfe5811.fer - ASCII file containing p-h LD groupped into 1 MeV (Fermi model)
C pfe5811.dat - ASCII file containing p-h LD groupped into 1 MeV bins(BCS model)
C pfe5811.nav - ASCII file containing p-h LD in 100 KeV bins(BCS model).
C pfe5811.p_h - Debug file. Contains all generated configurations
C               with energy less than 17 MeV(only present when DEBUG key=1).
C pfe5811.out - Detailed ASCII output file of the level density calculations.
C              Contains parity distribution, level and state density for each
C              parity, momentum distribution, brief information about BCS
C              solution, etc.
C
      Program Densidad
      Implicit Logical*1(A-Z)
      Integer*4 Dim
      Parameter (Dim=500)        
      Integer*4 Ja,Jp,Nlv,NP,NH,Kind,Calc,Olv,MaxEner
      logical Defor,DEBUG
      Real*8 Gi,BE,Eog,Elv,Slv,EnGrGAS,EnGrBCS,SS
      Character*7 FileName   
      Character*4 Root
      Integer*4 NEN,NER,IERM,I
      Real*8 SFN,SFER
      Common/ERR/  NEN,NER,SFN,SFER,IERM(25)
      Common/File/FileName
      Common/Ground/EnGrGAS,EnGrBCS,DEBUG
      Common/Data/ Ja,Jp,Nlv,Gi,NP,NH,Kind,BE,Eog,Calc,MaxEner,defor,SS
      Common/SPLV/ Elv(Dim),Slv(Dim),Olv(Dim)
      Write(*,1000)
 1000 Format('  Nucleus(4 lett.:ni58,s32n,pu40) => ')
      Read(*,'(A)') Root
      
      do 120 i=1,6
            if(i.eq.1) then
              Filename='n'//root(1:4)//'11'
            elseif(i.eq.2) then
              Filename='p'//root(1:4)//'11'
            elseif(i.eq.3) then
              Filename='n'//root(1:4)//'22'
            elseif(i.eq.4) then
              Filename='p'//root(1:4)//'22'
            elseif(i.eq.5) then
              Filename='n'//root(1:4)//'33'
            elseif(i.eq.6) then
              Filename='p'//root(1:4)//'33'
            endif
           
            Open(1,File=FileName//'.inp',Status='OLD',Err=100)
            Open(3,File=FileName//'.out')
 130        Call SPL
              
              IF(NP.GT.3.OR.NH.GT.3) STOP 'MAX(P,H)>3'
              
              NEN=0
              NER=0
              SFER=0.
              SFN=0.
              Call PartHole
              Call OutPut
              If(NP+NH.Eq.0.or.Mod(NP+NH,2).NE.Mod(Jp,2)) goto 130
100           Close(1)
              Close(3)
120   continue
      End
      Subroutine Config
      Implicit Logical*1(A-Z)
      Integer*4 Dim
      Parameter (Dim=500)
C
C   Controls configuration generation, BCS solution and energy
C   determination; calculates parity and spin projection of the
C   states. Only one exciton configuration are considered.
C
      Integer*4 I,J,Ien,Ienn,K,KQ,JS,NE,NConfig
      Integer*4 Ja,Jp,Nlv,NP,NH,Kind,Calc,Fermi,Itmp(10)
      Integer*4 Olv,LvsHole(5),LvsPart(5),MaxEner,Ma,NewPP,NewHH
      Integer*4 Parity,Old,New,OldP,OldH,NewP,NewH,IPS,SHM,DenBCS,DenGas
      Real*8 Gi,BE,Eog,Elv,Slv,Emax,EHM,CHI,Energy,EPP,EHH
      Real*8 EGas,PMP(10),Gap,Lamda,EnState,PSH,SJ,EnOLD
      Logical START, Defor, Impar, DEBUG
      Real*8 EnGrGAS,EnGrBCS
      Integer*4 NEN,NER,IERM
      Real*8 SFN,SFER,SS
      Common/ERR/NEN,NER,SFN,SFER,IERM(25)
      Common/Cant_Conf/NConfig
      Common/Ground/EnGrGAS,EnGrBCS,DEBUG
      Common/Data/ Ja,Jp,Nlv,Gi,NP,NH,Kind,BE,Eog,Calc,MaxEner,defor,SS
      Common/SPLV/ Elv(Dim),Slv(Dim),Olv(Dim)
      Common/limit/ Emax,EHM,SHM(101),DenBCS(400,30,2),DenGas(400,30,2)
      Data START/.TRUE./
      Impar=.False.
      If(Mod(Jp,2).Eq.1) Impar=.TRUE.
      Fermi=Jp/2+Mod(Jp,2)
      
      Do 100 K=1,5
       LvsPart(K)=0
       LvsHole(K)=0
100   continue
      
      DO 110 I=1,NP
       Olv(Fermi+I)=1
       LvsPart(I)=I
110   continue
      DO 120 I=1,NH
       Olv(Fermi-I+1)=1
       LvsHole(I)=I
120   continue
      NewHH=-Fermi
      EHH=Elv(Fermi)
      
C-----GENERATION OF PARTICLE CONFIGURATION--------------------------
 1    If(NP.Ne.0)Then
       CALL PNXCB(LvsPart,NP,Nlv-Fermi+1,OldP,NewP)
       Old=OldP+Fermi-1
       New=NewP+Fermi-1
       NewPP=New  
       EPP=Elv(NewPP)
       Olv(Old)=0
       Olv(New)=1               
C      You can not fill FERMI level for EVEN NUCLEUS
       If(.NOT.Impar.AND.New.Eq.Fermi) GOTO 3        
      Endif
C-----GENERATION OF HOLE CONFIGURATION------------------------------
 2    If(NH.Ne.0)Then
       CALL PNXCB(LvsHole,NH,Fermi,OldH,NewH)
       Old=Fermi-OldH+1       
       New=Fermi-NewH+1               
       NewHH=-New 
       EHH=Elv(New)
       Olv(Old)=0
       Olv(New)=1                                 
       If(Impar.AND.New.Eq.Fermi) GOTO 3
      Endif           
      
      J=0
      DO 130 I=1,Nlv
       If(Olv(I).Eq.1) J=J+1                          
130   continue
      
      If(Impar) then      
        If(J.Eq.NP+NH-1.AND.Olv(Fermi).Eq.1) GOTO 3        
        If(J.Eq.NP+NH-1) then
         Olv(Fermi)=1
         J=J+1
        Endif 
      Endif                          
      
      If(J.Ne.(NP+NH)) GOTO 3              
      
C     Excitation energy calculation
      EGas=Energy(0.D0,1.D0)-EnGrGAS                          
C     MAXIMUM ENERGY CUT-OFF
      If(EGas.GT.Dble(MaxEner)) GOTO 3
C
C-------------------------------------------------------------------------------------                                 
C     ONLY FOR DEBUG PRINTING
      IF(DEBUG) THEN
       J=0
       DO 140 I=1,Nlv
        If(Olv(I).Eq.1) then
         J=J+1
         If(Impar) then
           If(I.Eq.Fermi) Itmp(J)=I
         Else
           If(I.Eq.Fermi) Itmp(J)=-I
         Endif
         If(I.Lt.Fermi) Itmp(J)=-I
         If(I.Gt.Fermi) Itmp(J)= I
        EndIf
140    continue
      ENDIF
      
C----------------------------------------------------
C----------Solution of BCS equations----------------
      If( START .or. (ABS(EGas-EnOLD).Gt.0.005) ) Then
       START=.FALSE.
       EnOLD=EGas
       Call Clamda(Gap,Lamda,CHI)
       EnState=Energy(Gap,Lamda)-EnGrBCS
       PSH=EnState-EGas
       IPS=PSH*10.D0+50.D0
       IF(IPS.LT.1 .OR. IPS.GT.100) IPS=101
      Endif
C-------------------------------------------------------------------------------------                                 
C     ONLY FOR DEBUG PRINTING
      IF(DEBUG .AND. EnState.Le.17) then
        Write(10,55863) (Itmp(I),I=1,J)
        WRITE(10,5587) NewHH,EHH,NewPP,EPP,EGas,Gap,Lamda,CHI,EnState
      EndIf            
C----------------------------------------------------
      If(EnState.Gt.0.) NConfig=NConfig+1
C------Determination of parity for the configuration---------
      Parity=1
      DO 150 I=1,NP
       Ma=LvsPart(I)+Fermi-1
       Parity=Parity*sign(1.d0,Slv(Ma))
       PMP(I)=Abs(Slv(Ma))
150   continue
      DO 160 I=1,NH
       Ma=Fermi+1-LvsHole(I)
       Parity=Parity*sign(1.d0,Slv(Ma))
       PMP(NP+I)=Abs(Slv(Ma))
160   continue
      IF(Parity.Eq.-1) Parity=2
C------Coupling of all spin projections to generate nuclear states----
      NE=NP+NH
      DO 170 I=1,2**(NE-1)
       KQ=I-1
       SJ=0.D0
       DO 180 J=1,NE
        SJ=SJ+(1-2*Mod(KQ,2))*PMP(J)
        KQ=KQ/2
180    continue
       JS=ABS(SJ)+1
       IF(JS.Le.30)Then
C
C       STORING OF THE NUCLEAR STATE IN THE APROPRIATE ELEMENT OF DMN
C
        If(EnState.Gt.0.) then
          NConfig=NConfig+1
          Ien=EnState*10+1
          DenBCS(Ien,JS,Parity)=DenBCS(Ien,JS,Parity)+2
          SHM(IPS)=SHM(IPS)+2.D0
          EHM=EHM+2.D0*PSH      
        EndIf 
        If(EGas.Gt.0.) then
          Ienn=EGas*10+1     
          DenGas(Ienn,JS,Parity)=DenGas(Ienn,JS,Parity)+2
        Endif
       Endif
170   continue
 3    If(NH.Ne.0)Then
       If(LvsHole(NH).NE.NH)Goto 2
      Endif                            
      Olv(Fermi)=0
      If(NP.Ne.0)Then
       If(LvsPart(NP).NE.NP)Goto 1
      Endif
      Do 190 I=1,Nlv
       Olv(I)=0   
190   continue
      
      Return                                         
55863 FORMAT(/1X,'Particles & Holes:',10(I3,2X))            
 5587 FORMAT(1X,'[',I3,':',F6.2,',',I2,':',F6.2,']','=> U(Gas)=',F6.2,
     > ' D=',F5.3,' L=',F6.2,' CHI=',G8.2,' U(BCS)=',F6.2)
      End
      Subroutine Promoted
C
C   Controls configuration generation, BCS solution and energy
C   determination; calculates parity and spin projection of the
C   states. Only excited pair configurations are considered.
      Implicit Logical*1(A-Z)
      Integer*4 Dim
      Parameter (Dim=500)
      Integer*4 Ipart,Ihole,I,Ien,Ienn
      Integer*4 Ja,Jp,Nlv,NP,NH,Kind,Calc,Fermi,MaxEner
      Integer*4 Olv,ParHole(2),ParPart(2),NewPP,NewHH,NConfig
      Integer*4 Old,New,OldP,OldH,NewP,NewH,IPS,SHM,DenBCS,DenGas
      Logical defor,Start,Impar,DEBUG
      Real*8 Gi,BE,Eog,Elv,Slv,Emax,EHM,CHI,Energy,EnOLD
      Real*8 EGas,Gap,Lamda,EnState,PSH,EPP,EHH
      Real*8 EnGrGAS,EnGrBCS
      Integer*4 NEN,NER,IERM
      Real*8 SFN,SFER,SS
      Common/ERR/NEN,NER,SFN,SFER,IERM(25)
      Common/Cant_Conf/NConfig
      Common/Ground/EnGrGAS,EnGrBCS,DEBUG
      Common/Data/ Ja,Jp,Nlv,Gi,NP,NH,Kind,BE,Eog,Calc,MaxEner,defor,SS
      Common/SPLV/ Elv(Dim),Slv(Dim),Olv(Dim)
      Common/limit/ Emax,EHM,SHM(101),DenBCS(400,30,2),DenGas(400,30,2)
      Impar=.False.
      Fermi=Jp/2+Mod(Jp,2)
      If(Mod(Jp,2).Eq.1) Impar=.TRUE.
      DO 20 Ipart=0,NP/2
      DO 10 Ihole=0,NH/2
        Do 100 I=1,2
         ParPart(I)=0
         ParHole(I)=0
100     continue
        do 110 i=1,nlv
          Olv(i)=0                     
110     continue
        
        Start=.TRUE.
        
        IF(Ipart.Eq.0.And.Ihole.Eq.0) Goto 10
          DO 120 I=1,Ipart
           Olv(Fermi+I)=2
           ParPart(I)=I               
120       continue
          
          Do 130 I=1,Ihole
           Olv(Fermi-I+1)=-2
           ParHole(I)=I
130       continue
          
          NewPP=0
          NewHH=0
          
C---------GENERATION OF PAIR CONFIGURATION (PARTICLES)
 1        IF(Ipart.Ne.0)Then
           CALL PNXCB(ParPart,Ipart,Nlv-Fermi+1,OldP,NewP)
           Old=OldP+Fermi-1
           New=NewP+Fermi-1
           NewPP=New 
           EPP=Elv(NewPP)           
           If(Old.Gt.Fermi) Olv(Old)=0
C          You can not fill FERMI level(EVEN)
           If(.NOT.Impar.AND.New.Eq.Fermi) GOTO 3 
           Olv(New)=2                                       
           
           EGas=Energy(0.D0,1.D0)-EnGrGAS
           IF(EGas.GT.Dble(MaxEner)) GOTO 3
           
          Endif
C---------GENERATION OF PAIR CONFIGURATION (HOLES)
 2        IF(Ihole.Ne.0)Then
           CALL PNXCB(ParHole,Ihole,Fermi,OldH,NewH)
           Old=Fermi-OldH+1
           New=Fermi-NewH+1
           NewHH=-New
           EHH=Elv(New)                      
           If(Impar) THEN
             If(Olv(Old).Eq.-2) Olv(Old)=0
           ELSE
             Olv(Old)=0
           Endif
           Olv(New)=-2
          Endif         
          
          IF(NP.Gt.2*Ipart .Or. NH.Gt.2*Ihole)Then
C---------GENERATION OF CONFIGURATIONS FOR EXCITONS REMAINING AFTER
C-------------------EXCLUSION OF THE PROMOTED PAIRS
           Call Config_PAR(Ipart,Ihole,ParPart,ParHole)
          Else
C----------SOLUTION OF BCS EQUATIONS FOR PURE PROMOTED PAIR STATES
           EGas=Energy(0.D0,1.D0)-EnGrGAS
           IF(EGas.GT.Dble(MaxEner).OR.EGas.Lt.0.005d0) GOTO 3
           
           If( START .or. (ABS(EGas-EnOLD).Gt.0.005) ) Then
                 START=.FALSE.
                 EnOLD=EGas
                 Call Clamda(Gap,Lamda,CHI)
                 EnState=Energy(Gap,Lamda)-EnGrBCS
                 PSH=EnState-EGas
                 IPS=PSH*10.D0+50.D0
                 IF(IPS.LT.1 .OR. IPS.GT.100) IPS=101
           Endif
C-------------------------------------------------------------------------------------                                 
C     ONLY FOR DEBUG PRINTING
           IF(DEBUG. AND. Ipart+Ihole.LE.6 .AND. EnState.Lt.17.00) then                      
             Do 140 I=1,Nlv
              If(Olv(I).Eq.-2) Write(10,55861) -I
              If(Olv(I).Eq.2) Write(10,55862) I
140          continue
             WRITE(10,5587) EGas,Gap,Lamda,CHI,EnState
           EndIf            
           
           IF(EnState.Gt.0.) then
             NConfig=NConfig+1
             Ien=EnState*10+1
             DenBCS(Ien,1,1)=DenBCS(Ien,1,1)+1
             SHM(IPS)=SHM(IPS)+1.D0
             EHM=EHM+PSH
           Endif          
 
           If(EGas.Gt.0.) then
             Ienn=EGas*10+1
             DenGas(Ienn,1,1)=DenGas(Ienn,1,1)+1
           Endif
             
          Endif
  3       If(IHole.Ne.0)Then
           If(ParHole(Ihole).Ne.Ihole)Goto 2
          Endif
          If(Ipart.Ne.0)Then
           If(ParPart(Ipart).Ne.Ipart)Goto 1
          Endif
  10  Continue
  20  Continue
      Return                                              
55861 FORMAT(/1X,'Hole     Pair:',I3)      
55862 FORMAT( 1X,'Particle Pair:',I3)      
 5587 FORMAT(1X,'U(Gas)=',F6.2,' D=',F5.3,' L=',F6.2,' CHI=',G8.2,
     >  ' U(BCS)=',F6.2)
      End
      Subroutine Config_PAR(Ipart,IHole,ParPart,ParHole)
      Implicit Logical*1(A-Z)
      Integer*4 Dim
      Parameter (Dim=500)
C
C   Controls configuration generation, BCS solution and energy
C   determination; calculates parity and spin projection of the
C   states. Configuration with two exciton in one and another in the second orbital
C   are considered.
C
      Integer*4 I,J,K,Ien,Ienn,KQ,JS,NE,NPares,JPP,JH,IP,IH
      Integer*4 Ja,Jp,Nlv,NP,NH,Kind,Calc,Fermi
      Integer*4 Olv,LvsHole(5),LvsPart(5),MaxEner,Itmp(10),ItmpP(4)
      Integer*4 IHole,IPart,ParPart(2),ParHole(2),Ma,NConfig
      Integer*4 Parity,Old,New,OldP,OldH,NewP,NewH,IPS,SHM,DenBCS,DenGas
      Real*8 Gi,BE,Eog,Elv,Slv,Emax,EHM,CHI,Energy,EnOLD
      Real*8 EGas,PMP(10),Gap,Lamda,EnState,PSH,SJ
      Logical START,Defor,Impar,DEBUG
      Real*8 EnGrGAS,EnGrBCS
      Integer*4 NEN,NER,IERM
      Real*8 SFN,SFER,SS
      Common/ERR/NEN,NER,SFN,SFER,IERM(25)
      Common/Cant_Conf/NConfig
      Common/Ground/EnGrGAS,EnGrBCS,DEBUG
      Common/Data/ Ja,Jp,Nlv,Gi,NP,NH,Kind,BE,Eog,Calc,MaxEner,defor,SS
      Common/SPLV/ Elv(Dim),Slv(Dim),Olv(Dim)
      Common/limit/ Emax,EHM,SHM(101),DenBCS(400,30,2),DenGas(400,30,2)
      START=.TRUE.                          
      Impar=.False.                  
      Fermi=Jp/2+Mod(Jp,2)
      If(Mod(Jp,2).Eq.1) Impar=.TRUE.
            
      Do 100 K=1,5
       LvsPart(K)=0
       LvsHole(K)=0
100   continue
      Do 110 K=1,2
       If(ParPart(K).Eq.0) ParPart(K)=1000
       If(ParHole(K).Eq.0) ParHole(K)=1000
110   continue
      
C-------------------------------------------------------------------------------------                                 
C     ONLY FOR DEBUG PRINTING
      IF(DEBUG) THEN
       I=0
       DO 120 Ien=1,2
        If(ParPart(Ien).Lt.1000) then
         J=(ParPart(Ien)+Fermi-1)
         I=I+1
         ItmpP(I)=J
        Endif  
        If(ParHole(Ien).Lt.1000) then
         K=(Fermi-ParHole(Ien)+1)
         I=I+1
         ItmpP(I)=-K
        Endif  
120    continue
       NPares=IPart+IHole      
      ENDIF
      
      J=0
      DO 130 I=1,NP-2*IPart
        J=J+1
135     continue
        IF(ABS(Olv(Fermi+J)).Eq.2) THEN
          J=J+1
          GOTO 135
        ENDIF
        Olv(Fermi+J)=1
        LvsPart(I)=I
130   continue
      
      J=0
      DO 140 I=1,NH-2*IHole
        J=J+1
145     continue
        IF(ABS(Olv(Fermi-J+1)).Eq.2) THEN
          J=J+1
          GOTO 145
        ENDIF
        Olv(Fermi-J+1)=1
        LvsHole(I)=I
140   continue
C-----GENERATION OF PARTICLE CONFIGURATION----------------
 1    If(NP-2*IPart.Ne.0)Then
       CALL PNXCB(LvsPart,NP-2*IPart,Nlv-Fermi+1-IPart,OldP,NewP)
       Old=OldP+Fermi-1
       New=NewP+Fermi-1
       Do 150 I=1,IPart
        If(New.Ge.ParPart(I)+Fermi-1) New=New+1
        If(Old.Ge.ParPart(I)+Fermi-1) Old=Old+1
        If(New.Gt.Nlv.OR.Old.Gt.Nlv) goto 3
150    continue
       
       If(Impar)then               
        If(Olv(Old).Eq.1) Olv(Old)=0
C       If(Old.Ge.Fermi)Olv(Old)=0                                          
       Else
        If(Old.Gt.Fermi)Olv(Old)=0                                         
C      You can not fill FERMI level(EVEN)         
        If(New.Eq.Fermi) GOTO 3
       Endif
       Olv(New)=1               
       
      Endif
C-----GENERATION OF HOLE CONFIGURATION------------------
 2    If(NH-2*IHole.Ne.0)Then
       CALL PNXCB(LvsHole,NH-2*IHole,Fermi-IHole,OldH,NewH)
       Old=Fermi-OldH+1
       New=Fermi-NewH+1
       Do 160 I=1,IHole
        If(New.Le.Fermi-ParHole(I)+1) New=New-1
        If(Old.Le.Fermi-ParHole(I)+1) Old=Old-1
        If(New.Lt.1.OR.Old.Lt.1) goto 3
160    continue
                                                                           
       If(Impar)then               
        If(Olv(Old).Eq.1) Olv(Old)=0
C      You can not get particle from the FERMI level(ODD)        
        If(New.Eq.Fermi) GOTO 3
       Else
        If(Old.Le.Fermi) Olv(Old)=0                                         
       Endif
       
       Olv(New)=1              
       
      Endif
      JH=0
      JPP=0
      IP=0
      IH=0
      IF(Impar) THEN
       DO 170 I=1,Fermi-1
        If(Olv(I).Eq.-2) IH=IH+1
        If(Olv(I).Eq.1) JH=JH+1                          
170    continue
       If(Olv(Fermi).Eq.-2) IH=IH+1              
       
       DO 180 I=Fermi,Nlv
        If(Olv(I).Eq.1) JPP=JPP+1                          
        If(Olv(I).Eq.2) IP=IP+1
180    continue
      ELSE
       DO 190 I=1,Fermi
        If(Olv(I).Eq.-2) IH=IH+1
        If(Olv(I).Eq.1) JH=JH+1                          
190    continue
       DO 200 I=Fermi+1,Nlv
        If(Olv(I).Eq.2) IP=IP+1       
        If(Olv(I).Eq.1) JPP=JPP+1                          
200    continue
      ENDIF                                  
      
      IF(Impar.AND.(2*IP+JPP.Eq.NP-1)) THEN
        If(Olv(Fermi).Eq.1) GOTO 3
        Olv(Fermi)=1
        JPP=JPP+1
      ENDIF
      
      If(JPP+2*IP.NE.NP.OR.JH+2*IH.NE.NH) GOTO 3                  
      
      EGas=Energy(0.D0,1.D0)-EnGrGAS                       
      If(EGas.Gt.Dble(MaxEner)) GOTO 3      
C-------------------------------------------------------------------------------------                                 
C     ONLY FOR DEBUG PRINTING
      IF(DEBUG) THEN   
       J=0                                         
       DO 210 Ien=1,Nlv
        If(Olv(Ien).Eq.1) then
         J=J+1                          
         If(Impar) then                         
           If(Ien.Eq.Fermi) Itmp(J)=Ien
         Else                     
           If(Ien.Eq.Fermi) Itmp(J)=-Ien
         Endif
         If(Ien.LT.Fermi) Itmp(J)=-Ien
         If(Ien.Gt.Fermi) Itmp(J)= Ien
        EndIf  
210    continue
      ENDIF         
        
C----------Solution of BCS equations----------------
      If( START .or. (ABS(EGas-EnOLD).Gt.0.005) ) Then
       START=.FALSE.
       EnOLD=EGas
       Call Clamda(Gap,Lamda,CHI)
       EnState=Energy(Gap,Lamda)-EnGrBCS
       PSH=EnState-EGas
       IPS=PSH*10.0+50
       IF(IPS.LT.1 .OR. IPS.GT.100) IPS=101
      Endif                 
      
C-------------------------------------------------------------------------------------                                 
C     ONLY FOR DEBUG PRINTING
      If(DEBUG .AND. EnState.Lt.17) then
        If(NPares.Gt.0) Write(10,55862) 
     >             (ItmpP(K),Elv(Abs(ItmpP(K))),K=1,NPares)
        If(J.Gt.0) Write(10,55863) (Itmp(I) ,I=1,J)
        WRITE(10,*) 'IP=',IP,' IH=',IH,' JPP=',JPP,' JH=',JH
        WRITE(10,5587) EGas,Gap,Lamda,CHI,EnState
      EndIf            
          
      If(EnState.Gt.0.) NConfig=NConfig+1
C------Determination of parity for the configuration---------
      Parity=1
      DO 220 I=1,NP-2*IPart
       Ma=LvsPart(I)+Fermi-1
       If(LvsPart(I).Ge.ParPart(1)) Ma=Ma+1
       If(LvsPart(I).Ge.ParPart(2)) Ma=Ma+1
       Parity=Parity*sign(1.d0,Slv(Ma))
       PMP(I)=Abs(Slv(Ma))
220   continue
      DO 230 I=1,NH-2*IHole
       Ma=Fermi+1-LvsHole(I)
       If(LvsHole(I).Ge.ParHole(1))Ma=Ma-1
       If(LvsHole(I).Ge.ParHole(2))Ma=Ma-1
       Parity=Parity*sign(1.d0,Slv(Ma))
       PMP(NP-2*IPart+I)=Abs(Slv(Ma))
230   continue
      IF(Parity.Eq.-1)Parity=2
C-----Coupling of all spin projections to generate nuclear states----
      NE=NP+NH-2*(IPart+IHole)
      DO 250 I=1,2**(NE-1)
       KQ=I-1
       SJ=0.D0
       DO 240 J=1,NE
        SJ=SJ+(1-2*Mod(KQ,2))*PMP(J)
        KQ=KQ/2
240    continue
       JS=ABS(SJ)+1
       IF(JS.Le.30)Then
        If(EnState.Gt.0.) then
         Ien=EnState*10+1    
         SHM(IPS)=SHM(IPS)+2.D0
         EHM=EHM+2.D0*PSH
         DenBCS(Ien,JS,Parity)=DenBCS(Ien,JS,Parity)+2
        EndIf
        If(EGas.Gt.0.) then
         Ienn=EGas*10+1     
         DenGas(Ienn,JS,Parity)=DenGas(Ienn,JS,Parity)+2
        Endif
       Endif
250   continue
 3    If(NH.Ne.2*IHole)Then
       If(LvsHole(NH-2*IHole).NE.NH-2*IHole)Goto 2
      Endif       
      If(Olv(Fermi).Eq.1) Olv(Fermi)=0      
      If(NP.Ne.2*IPart)Then
       If(LvsPart(NP-2*IPart).NE.NP-2*IPart)Goto 1
      Endif
      
      Do 260 I=1,Nlv
       If(Olv(I).Eq.1) Olv(I)=0   
260   continue
      
      Return
55862 FORMAT(/1X,'PAIRS:',4('[',I3,':',F6.2,'] ') )      
55863 FORMAT( 1X,'Particles & Holes:',10(I3,2X))      
 5587 FORMAT( 1X,' U(Gas)=',F6.2,
     > ' D=',F5.3,' L=',F6.2,' CHI=',G8.2,' U(BCS)=',F6.2)
      End
      Subroutine PartHole
      Implicit Logical*1(A-Z)
      Integer*4 Dim
      Parameter (Dim=500)
      Character*7 Text(2)
      Integer*4 I,J,NConfig
      Integer*4 Olv,SHM,DenBCS,DenGas,MaxEner
      Integer*4 Ja,Jp,Nlv,NP,NH,Kind,Calc,Cuasip,Num
      logical defor,Impar,DEBUG          
      Character*7 FileName         
      Real*8 Gi,BE,Eog,EHM,CHI
      Real*8 Elv,Slv,Gap,Lamda,EnGrBCS,Emax,EGapGAS,EGapBCS
      Real*8 EnGrGAS,Energy
      Integer*4 NEN,NER,IERM,Fermi
      Real*8 SFN,SFER,SS  
      Common/File/FileName
      Common/ERR/NEN,NER,SFN,SFER,IERM(25)
      Common/Cant_Conf/NConfig
      Common/Ground/EnGrGAS,EnGrBCS,DEBUG
      Common/Data/ Ja,Jp,Nlv,Gi,NP,NH,Kind,BE,Eog,Calc,MaxEner,defor,SS
      Common/SPLV/ Elv(Dim),Slv(Dim),Olv(Dim)
      Common/limit/ Emax,EHM,SHM(101),DenBCS(400,30,2),DenGas(400,30,2)
      Data Text/'NEUTRON','PROTON '/,Impar/.FALSE./
      Cuasip=NP*10+NH
      Fermi=Jp/2+Mod(Jp,2)           
      If(Mod(Jp,2).Eq.1) Impar=.TRUE.    
      
      Write(3,1001)Text(Kind)
      DO 100 I=1,400
        DO 100 J=1,30
         DenBCS(I,J,1)=0
         DenBCS(I,J,2)=0
         DenGas(I,J,1)=0
         DenGas(I,J,2)=0
100   continue
      
      NConfig=0
      EHM=0.D0
      Do 110 I=1,101
       SHM(I)=0
110   continue
      
      Do 120 I=1,Nlv
       Olv(I)=0
120   continue
      If(Impar) Olv(Fermi)=1
      
C---------BCS Solution for the Ground State----------------------
      Call Clamda(Gap,Lamda,CHI)
      EnGrBCS=Energy(Gap,Lamda)
      EnGrGAS=Energy(0.D0,1.D0)
      Write(3,1002)EnGrGAS,Gap,Lamda,CHI,EnGrBCS,EnGrGAS-EnGrBCS
      If(Cuasip.Eq.0) return
      If(Mod(NP+NH,2).ne.Mod(Jp,2))Then
       Write(3,1004)
       return
      Endif
C----Calculation of the first excited state (without promoted pairs)----
      If(Impar) then
        
        If(NP.Gt.1) then
          
          Num=NP-1
115       continue
          IF(Num.Gt.0) then
           J=1
125        continue
           IF(Elv(Fermi+J).EQ.Elv(Fermi)) then
            J=J+1
            goto 125
           Endif
135        continue
           if(Olv(Fermi+J).Eq.1) then
              J=J+1
              goto 135
           Endif
           Olv(Fermi+J)=1
           Num=Num-1
           goto  115
          Endif
      
          Do 137 I=1,NH
           Olv(Fermi-I)=1
137       continue
          
        Else
        
          J=1
145       continue
          IF (Elv(Fermi+J).EQ.Elv(Fermi)) then
              J=J+1
              goto 145
          Endif
          
          Olv(Fermi)=0
          Olv(Fermi+J)=1
          
        EndIf
        
        
      Else
        
        Num=NP
155     continue
        IF(Num.Gt.0) then
         J=1
165      continue
         IF(Elv(Fermi+J).EQ.Elv(Fermi)) THEN
          J=J+1
          goto 165
         Endif
175      continue
         IF(Olv(Fermi+J).Eq.1) then
            J=J+1
            goto 175
         Endif
         Olv(Fermi+J)=1
         Num=Num-1
         goto 155
        Endif
        Do 160 I=1,NH
         Olv(Fermi-I+1)=1
160     continue
        
      Endif
C
C     BCS solution to calculate Gap and Lamda 
C
      Call Clamda(Gap,Lamda,CHI)
      EGapBCS=Energy(Gap,Lamda)-EnGrBCS
      EGapGAS=Energy(0.D0,1.D0)-EnGrGAS
      Write(3,1005)EGapGAS,Gap,Lamda,CHI,EGapBCS,EGapBCS-EGapGAS
C-----Generation of the configurations and densities calculation-------
      Do 140 I=1,Nlv
       Olv(I)=0
140   continue
                     
      NEN=0        
      NConfig=0       
C-------------------------------------------------------------------------------------                                 
C     ONLY FOR DEBUG PRINTING
      IF(DEBUG) open(10,file=filename//'.p_h')
      IF(DEBUG) write(10,1006)
      Call Config()
      Call Promoted 
      IF(DEBUG) close(10)      
      
      Return
 1001 Format('  Solution for ',A7,' System')
 1002 Format(/'-------------Solution for BCS Ground State----------',/
     >     '  Energy of Free Gas Ground State = ',F12.5,' MeV',/
     >     '  Gap= ',F10.5,'  Lamda= ',F10.5,' MeV','  CHI= ',E10.3/,
     >     '    Energy of BCS Ground State = ',F12.5,' MeV',/
     >     '           Condensation Energy = ',F10.5,' MeV')
 1004 Format('   Only an IDIOT can look for odd configurations '/
     >       '        in even system')
 1005 Format(/'-----------Solution for First excited state----------',/
     >        '   First Excited Energy of Free Gas'/
     >        '     (without promoted pairs)=',F12.5,' MeV',/
     >     '  Gap= ',F10.5,'  Lamda= ',F10.5,' MeV','  CHI= ',E10.3/,
     >        '   First Excited Energy with Pairing Force'/
     >        '     (without promoted pairs)=',F12.5,' MeV',/
     >        '             Pairing Gap=',F10.5,' MeV')
 1006 Format(1x,
     > ' Here is an enumeration of all generated configurations'/
     > '  U(Gas) is the configuration excitation energy in the FGM'/
     > '  U(BCS) is the configuration excitation energy in the BCSM'/
     > '  D and L are DELTA and LAMBDA of the BCS solution'/ )
      End
      Subroutine CLamda(Gap,Lamda,CHI)
      Implicit Logical*1(A-Z)
      Integer*4 Dim,IterMax,Ni,Nf,Ncalc
      Real*8 Eps
      Parameter (Dim=500,IterMax=400,Eps=1.D-7)
      Integer*4 Ja,Jp,Nlv,NP,NH,Kind,Calc
      Integer*4 Olv,I,MaxEner
      Integer*4 Iter,Fermi,Num
      logical defor
      Real*8 Gi,BE,Eog,Elv,Slv,Gap,Lamda,Gap2,LamdaI,Lambda
      Real*8 DenE,Ex,CHI,A1,A2,B1,B2,C1,C2,EFermi
      Real*8 DelLam,DelGap,G,RelE,Epsj,Det,DLLast,DelMax
      Integer*4 NEN,NER,IERM,IER
      Real*8 SFN,SFER,SS
      Common/ERR/  NEN,NER,SFN,SFER,IERM(25)
      Common/Data/ Ja,Jp,Nlv,Gi,NP,NH,Kind,BE,Eog,Calc,MaxEner,defor,SS
      Common/SPLV/ Elv(Dim),Slv(Dim),Olv(Dim)
      NEN=NEN+1
      Gap=0.D0
      
      Fermi=Jp/2+Mod(Jp,2)
105   continue
      IF(Elv(Fermi).Eq.Elv(Fermi+1)) then
       Fermi=Fermi+1
       goto 105
      Endif
      
      EFermi=Elv(Fermi)
      
      Ni=1
      Nf=Nlv
       
      Ncalc=Nf-Ni+1
      
C------ Initial Guess ------------
      Lambda=Fermi                                         
      LamdaI=0.5D0*(Elv(Fermi)+Elv(Fermi+1))       
      
      If(Gi.Eq.0.D0) Return
      
      Num=0
      Do 100 I=Ni,Nf
       If(Olv(I).Eq.1)Num=Num+1
       If(Olv(I).Eq.2)Num=Num+2
100   continue
      DenE=(Elv(Nf)-Elv(Ni))/Ncalc
      DelMax=2.D0*DenE
      G=Gi/Dble(Ja)
      CHI=0.D0
      Iter=0
      C1=0
      C2=(G/2.D0*Ncalc)**2
      DO 110 Iter=1,30
       Gap2=(C1+C2)/2.D0
       IF(DABS(C1-C2).LE.0.1D0)Goto 2
       B1=-2.D0/G                        
        DO 120 I=Ni,Nf
        IF(Olv(I).Eq.0)Then
         RelE=Elv(I)-LamdaI
         Epsj=Sqrt(Gap2+RelE*RelE)
         If(Epsj.Lt.0.005) goto 125
         B1=B1+1.D0/Epsj
        Endif
120     continue
125     IF(B1.Le.0.D0)Then
        C2=Gap2
       Else
        C1=Gap2
       Endif
110   continue
      Ex=EXP(DenE/G)
      Gap2=Sqrt(Dble(Jp*(2*Ncalc-Jp)))*Ex/(Ex**2-1.)*DenE
      Gap2=Gap2*Gap2
 2    DLLast=DelMax
 1    Iter=Iter+1
      A1=0.D0
      A2=0.D0
      B1=2.D0/G
      B2=Dble(Jp-Num)
      Do 130 I=Ni,Nf
       If(Olv(I).Eq.0)Then
        RelE=Elv(I)-LamdaI
        Epsj=Sqrt(Gap2+RelE*RelE)
        
        IF(Epsj.Lt.0.005) GOTO 20
        
        B1=B1-1.D0/Epsj
        B2=B2-(1.D0-RelE/Epsj)
        A1=A1+1.D0/Epsj**3
        A2=A2+1.D0*RelE/Epsj**3
       Endif
130   continue
      CHI=Abs(B1)+Abs(B2)
      If(CHI.Le.Eps) THEN
        SFN=SFN+CHI
        Goto 20
      ENDIF
      
      Det=0.5D0*(Gap2*A1*A1+A2*A2)
      If(Abs(Det).Lt.1.D-10) Det=Sign(1.D-3,Det)
      DelLam=0.5D0*(A1*B2+A2*B1)/Det
      DelGap=(A2*B2-Gap2*A1*B1)/Det
      If(DLLast*DelLam.Lt.0.D0.And.Abs(DelLam).Gt.Abs(DLLast))Then
       DelLam=-0.5D0*DLLast
       DelGap=-0.5D0*DelGap*DLLast/DelLam
      Endif
      DLLast=DelLam
      If(Abs(DelGap).Gt.DelMax) DelGap=Sign(DelMax,DelGap)
      If(Abs(DelLam).Gt.DelMax) DelLam=Sign(DelMax,DelLam)
      Gap2=Gap2+DelGap
      LamdaI=LamdaI+DelLam              
      If(Iter.Gt.IterMax) Goto 20
      If(Gap2.Lt.0.D0) Gap2=0.5D0*(Gap2-DelGap)
      If(Iter.Le.IterMax) Goto 1
 20   Continue
      If(Gap2.Lt.0.05)Then
        IF(CHI.Gt.0.001) then
         I=I+1
        Endif 
        NER=NER+1
        IER=ABS(3.0-LOG10(CHI))
        IF(IER.LT.1) IER=1
        IF(IER.GT.25) IER=25
        IERM(IER)=IERM(IER)+1
        SFER=SFER+CHI
        Gap=0.D0
        Lamda=0.5D0*(Elv(Fermi)+Elv(Fermi+1))
      Else
        Gap=Sqrt(Abs(Gap2))
        Lamda=LamdaI
      Endif                             
      
      Return
      End
      Function Energy(Gap,Lamda)
      Implicit Logical*1(A-Z)
      Integer*4 Dim
      Parameter (Dim=500)
      Integer*4 Ja,Jp,Nlv,NP,NH,Kind,Calc,Fermi
      Integer*4 Olv,I,Num,MaxEner
      logical defor
      Real*8 Gi,BE,Eog,Elv,Slv,Gap,Lamda,BCS
      Real*8 G,RelE,Epsj,Energy,SS
      Common/Data/ Ja,Jp,Nlv,Gi,NP,NH,Kind,BE,Eog,Calc,MaxEner,defor,SS
      Common/SPLV/ Elv(Dim),Slv(Dim),Olv(Dim)
      
      Fermi=Jp/2+Mod(Jp,2)        
      G=Gi/Dble(Ja)
      
      If(Gap.Eq.0.D0)Then
C
C     FERMI GAS MODEL
C       
       Num=Jp
       BCS=0.D0
       Do 100 I=1,Nlv
        If(Olv(I).Eq.1)Then
         BCS=BCS+Elv(I)
         Num=Num-1
        Endif
        
        If(Olv(I).Eq.2)Then
         BCS=BCS+2.D0*Elv(I)
         Num=Num-2
        Endif
100    continue
       I=0
105    continue
       IF(Num.NE.0) then
          I=I+1
          
          If(Olv(I).Ne.0) goto 105
          
          If(Num.Eq.1)Then
              BCS=BCS+Elv(I)
              Num=0
          Else
              BCS=BCS+2.D0*Elv(I)
              Num=Num-2
          Endif
          goto 105
       ENDIF
      Else
C
C     SUPERFLUID MODEL
C
       BCS=-Gap*Gap/G
       Do 120 I=1,Nlv
        If(Olv(I).Eq.0) Then
         RelE=Elv(I)-Lamda
         Epsj=Sqrt(Gap*Gap+RelE*RelE)
         BCS=BCS+Elv(I)*(1.D0-RelE/Epsj)            
        Endif                               
        If(Olv(I).Eq.1) BCS=BCS+Elv(I)
        If(Olv(I).Eq.2) BCS=BCS+2.D0*Elv(I)
120    continue
      Endif
      
      Energy=BCS
      Return
      End
      Subroutine SPL
      Implicit Logical*1(A-Z)
      Integer*2 Dim
      Parameter (Dim=500)
      Character*2 Smat,Type*1
      Character*10 Text1(2),Name*12,Title*80
      Character*69 Text2(-1:1)
      Integer*2 I,J,K,ITMP
      Integer*4 Ja,Jp,Jz,Nlv,NP,NH,Kind,Calc,Fermi,NlvTmp,FermiT
      Integer*4 Olv,MaxEner
      logical defor,DEBUG
      Real*8 Gi,BE,Eog,Elv,Slv,STMP,PAR,SS
      Real*8 EnGrGAS,EnGrBCS
      Common/Ground/EnGrGAS,EnGrBCS,DEBUG
      
      Common/Data/ Ja,Jp,Nlv,Gi,NP,NH,Kind,BE,Eog,Calc,MaxEner,defor,SS
      Common/SPLV/ Elv(Dim),Slv(Dim),Olv(Dim)
      Data Text1/' Neutrons ',' Protons  '/
      Data Text2/' Equidistant spacing model(ESM)            ',
     >           ' Shell Model spectra calculation           ',
     >           ' ESM with spacing from SHELL model         '      /
      
      Read(1,'(A1)')Type
      If(Type.eq.'N'.or.Type.eq.'n')Then
       Kind=1
      else
       Kind=2
      endif
      Read(1,*)NP,NH
      Read(1,*)Ja,Jz
      Read(1,*)Nlv,BE
c     CALC = 0 =>ESM with Eog from input
c            1 =>SHELL spectrum 
c            2 =>ESM with Eog from SHELL spectrum
c
c     Eog  =  INPUT SINGLE-PARTCILE SPACING
c     ITMP = 0 (NO DEBUG)
      Read(1,*,ERR=123)Calc,Eog,ITMP
123   Read(1,*)Gi
      Read(1,*)MaxEner          
      
      DEBUG=.FALSE.
      If(ITMP.Eq.1) DEBUG=.TRUE.
      
      If(Kind.eq.1)Then
       Jp=Ja-Jz
      Else
       Jp=Jz
      EndIf
      K=Calc-1
      IF(K.GT.1) K=1
      IF(K.LT.-1) K=-1
c      K=-1 => ' Equidistant model( Spacing from input file)',
c      K=0  => ' Shell Model without equidistant levels          '
c      K=+1 => ' Equidistant model with spacing from SHELL model '
      Write(3,1000)Text2(K),Text1(Kind),Smat(Jz),Ja,Ja,Jz
      Write(3,1007)NH,NP            
      Write(*,1000)Text2(K),Text1(Kind),Smat(Jz),Ja,Ja,Jz
      Write(*,1007)NH,NP                  
      
      
      If(Calc.Ne.-1)Then
C
C     SHELL-SPECTR input
C
         Read(1,'(A)')Name
         Open(2,File=Name,Status='OLD',err=126)
         Write(3,1006)
         Read(2,1005)Title
         read(2,1086)nlvtmp
         
         IF(Kind.EQ.2) THEN
          DO 100 I=1,nlvtmp
              Read(2,1005)Title
100       continue
          read(2,1086)nlvtmp                        
         ENDIF
         
         if(nlvtmp.Lt.Nlv)then
           Write(*,*)'  You ask more states than exist in spectr file'
           stop
         Endif
         Do 110 I=1,nlvtmp
          Read(2,1111) Elv(I),STMP,PAR
          Slv(I)=PAR*STMP
110      continue
         Close(2)
         defor=.true.
         Do 120 I=2,Nlv
          If(Elv(I).eq.Elv(I-1)) THEN
              Defor=.false.
              goto 125
          EndIf    
120      continue
125      continue
      Else
C
C     EQUIDISTANT SPACING MODEL
C
       If(Eog.Eq.0.D0) Stop 'ESM selected and Eog = 0'
      EndIf          
      
      If(Calc.Ne.-1) Eog=Dble(Nlv)/(Elv(Nlv)-Elv(1))
C
C     D = 2./Eog 
C
      Write(3,1004) 2.D0/Eog
      Write(3,1008)
      If(Calc.gt.1)Then
       Elv(1)=0.D0
       Do 130 I=2,nlv
        Elv(I)=Elv(I-1)+2.D0/Eog
        Slv(I)=0.5
130    continue
      EndIf
      
      Fermi=Jp/2+Mod(Jp,2)                                          
      FermiT=Fermi
105   continue
      IF(Elv(Fermi).Eq.Elv(Fermi+1)) then
       Fermi=Fermi+1
       goto 105
      Endif 
      
      If(BE.Eq.0.D0)Then
       BE=Elv(Nlv)-Elv(Fermi)
      Else
       I=Fermi
 2     I=I+1
       If(Elv(I)-Elv(Fermi).Le.BE.And.I.Lt.Nlv) Goto 2
       Nlv=I
      Endif
C
C     CHECK FOR THE SPIN (HALF OR INTEGER); SS
C
      SS=1.0
      IF(MOD((NP+NH),2).EQ.1) SS=0.5
      
      Do 140 I=1,Nlv,4
       Write(3,1002)(J,Elv(J),Slv(J),J=I,Min(I+3,Nlv))
140   continue
      
      Write(3,1008)                                    
      Write(3,1010) FermiT,Fermi,Elv(Fermi)
      
      if(defor) then
        write(3,10031)
      else
        write(3,10032)
      endif
      Write(3,1009)BE,Nlv
      
C     MaxEner=1.2*BE
      
      Return                                   
 126  stop 'SPECTRUM FILE NOT FOUND'    
 1000 Format(/'  CALCULATION FOR ',A69,/
     >  '  FOR',A10,' of NUCLEUS ',A2,'-',i3,'    A=',I3,' Z=',I3/)
 1001 Format(3(1x,I3,'!',3x,F8.3,2x,I3,4x))
 1002 Format(4(1x,I3,1x,F7.3,1x,F4.1,' |'))
 1086 FORMAT(3I10)
 1111 FORMAT(E15.5,F10.1,F5.0)
10031 FORMAT(1X,/,10X,' NUCLEUS IS TREATED AS DEFORMED'/)
10032 FORMAT(1X,/,10X,' NUCLEUS IS TREATED AS SPHERICAL'/)
 1004 Format('      Spacing = ',f10.3,' MeV')
 1005 Format(A80)
 1006 Format(15('-'),' Input Data of Shell Model Spectr ', 15('-'))
 1007 Format('  Configuration with ',i2,' Holes and ',i2,' Particles'//)
 1008 Format(80('-'))
 1009 Format(' Gap of energy relative to Fermi level',/,
     >       '    below which orbitals are considered=',F10.5,' MeV',/,
     >       ' Number of orbitals wich are considered=',i4,//)
 1010 Format(/10x,'True Fermi Level Number=',I3,/,
     >        10x,'Shifted Fermi Level Number=',I3,/,
     >        10x,'Fermi Level Energy=',F10.5,' MeV'/)
      End
      Subroutine OutPut
C     Implicit Logical*1(A-Z)
      Integer*2 Dim
      Parameter (Dim=500)
      Integer*4 I,J,K,Fermi,Ja,Jp,Nlv,NP,NH,Kind,Calc,NC
      Integer*4 Olv,MaxEner,SHM,DenBCS,DenGas,Jz,NConfig
C     REAL*4 GIFFF
      logical defor,DEBUG
      Real*8 Gi,BE,Eog,Elv,Slv,Emax,EHM,EnGrGAS,EnGrBCS
      REAL*4 ERS(101),TOTSD,TOTD,ASA,EOM(40,2),PARM(40),EOP,EOMS
      REAL*4 SIG(40,2),SPIN(30),ENERG(40),ERO(40,2)
      DIMENSION SAVE(2),SAVEN(2),SAVENA(2),TRC(2),SLOPE(2)
      INTEGER*4 IEL,IEM,IE,IK,ILI,KMAX2
      Character*7 FileName
      Character*2 SMat
      Integer*4 NEN,NER,IERM
      Real*8 SFN,SFER,SS
      Common/ERR/  NEN,NER,SFN,SFER,IERM(25)            
      Common/Cant_Conf/NConfig
      Common/File/FileName
      Common/Ground/EnGrGAS,EnGrBCS,DEBUG
      Common/Data/ Ja,Jp,Nlv,Gi,NP,NH,Kind,BE,Eog,Calc,MaxEner,defor,SS
      Common/SPLV/ Elv(Dim),Slv(Dim),Olv(Dim)
      Common/limit/ Emax,EHM,SHM(101),DenBCS(400,30,2),DenGas(400,30,2)
      SFN=SFN/(NEN-NER)
      IF(NER.NE.0) SFER=SFER/NER
      WRITE(3,798) NEN,SFN,NER,SFER
  798 FORMAT(1X/1X,I8,' MINIMIZATION ENTRIES'/
     X          9X,   ' AVERAGE CHI.SQ. FOR NORMAL SOLUTIONS=',G15.5,/
     X         /1X,I8,' ENTRIES WITHOUT REAL SOLUTIONS',/
     X          9X,   ' AVERAGE CHI.SQ. FOR UNREAL SOLUTIONS',G15.5)
      IF(NER.GT.0) THEN
        WRITE(3,814)
  814   FORMAT(1X,/,3X,'CHI.SQUARE SPECTRUM OF UNREAL SOLUTIONS'/)
        DO 816 I=1,25
        J=3-I
  815   FORMAT(15X,'1.0E',I3,I15)
  816   WRITE(3,815) J,IERM(I)
      ENDIF
      TOTSD=0.0                     
      KMAX=40
      IEL=0
      IEM=0
C
C   PRINTOUT OF BCS ENERGY SHIFTS
C
      DO 666 I=1,101
      IF(I.GT.100) GO TO 666
      ERS(I)=-4.95+0.1*FLOAT(I)
      IF(IEL.EQ.0 .AND. SHM(I).NE.0.0) IEL=I
      IF(SHM(I).NE.0.0) IEM=I
      TOTSD=TOTSD+SHM(I)
  666 CONTINUE
      TOTD=TOTSD+SHM(101)
      IF(IEL.GT.1) IEL=IEL-1
      IF(IEL.GT.1) IEL=IEL-1
      IF(IEM.LT.100) IEM=IEM+1
      IE=IEM-IEL+1
      DO 667 I=1,IE
      ERS(I)=ERS(I+IEL-1)
  667 SHM(I)=SHM(I+IEL-1)
      ERS(IE+1)=0.0
      SHM(IE+1)=0.0
      WRITE(6,698)
      WRITE(3,698)
  698 FORMAT(1X,//,3X,'HISTOGRAM OF ENERGY SHIFTS ',/,
     X11X,'MEV',8X,'BCS-FG',/)
      WRITE(6,662) (ERS(I),SHM(I),I=1,IE)
      WRITE(3,662) (ERS(I),SHM(I),I=1,IE)
  662 FORMAT(3X,F6.2,2X,I9)
      DO 668 I=1,IE
  668 ERS(I)=ERS(I)+0.05
      IF(TOTSD.GT.0) EHM=EHM/TOTSD
      IF(GI.EQ.0.0) EHM=0.0
      WRITE(6,773) EHM
      WRITE(3,773) EHM
  773 FORMAT(1X,/,3X,'AVERAGE ENERGY SHIFT BCS-FG=',F6.2,/)
      WRITE(6,949) TOTD
      WRITE(3,949) TOTD
 949  FORMAT(3X,'TOTAL NUMBER OF STATES(BCS)=',G15.5,/)
      WRITE(6,9491) NConfig
      WRITE(3,9491) NConfig
 9491 FORMAT(3X,'TOTAL NUMBER OF CONFIGURATIONS=',I10/)
 
      IF(TOTD.Eq.0) GOTO 54321 ! NO BCS STATES , only FERMI GAS CALC.
      If(Kind.eq.1)Then
       Jz=Ja-Jp
      Else
       Jz=Jp
      EndIf               
      
      Fermi=Jp/2+Mod(Jp,2)
C      These lines are intended to prepare an input binary file
C      for convolution code.
C      open(21,file=filename//'.con',form='unformatted')
C      if(defor)kind=-kind
C      GIFFF=GI
C      write(21)Ja,Jz,NP,NH,Kind,KMAX*10,GIFFF,FLOAT(Nlv-Fermi),
C     >     (((Float(DenBCS(I,J,K)),J=1,30),I=1,KMAX*10),K=1,2)
C      Close(21)
C
C      open(21,file=filename//'.cof',form='unformatted')
C      if(defor)kind=-kind
C      GIFFF=GI
C      write(21)Ja,Jz,NP,NH,Kind,KMAX*10,GIFFF,FLOAT(Nlv-Fermi),
C     >     (((Float(DenGas(I,J,K)),J=1,30),I=1,KMAX*10),K=1,2)
C      Close(21)
      Open(18,File=FileName//'.nav')
      Do 115 I=1,MAXENER*10
       D3=0.D0
       D4=0.D0
       DO 120 J=1,30
         DO 120 K=1,2
           D3=D3+DenGas(I,J,K)
           D4=D4+DenBCS(I,J,K)
120    continue
       If(D3.eq.0.D0)D3=0.1D0
       If(D4.eq.0.D0)D4=0.1D0
       Write(18,'(1x,F6.2,1X,2(F14.2,1x))') (I-1)*0.1,D3,D4
       Write(18,'(1x,F6.2,1X,2(F14.2,1x))') I*0.1,D3,D4
115   continue
      CLOSE(18)
      NTOT=0
      DO 333 K=1,2
      Neg_Lev=0
      WRITE(3,334)
  334 FORMAT(1X,//)
      IF(K.EQ.1) WRITE(3,335)
  335 FORMAT(/20X,'P O S I T I V E   P A R I T Y')
      IF(K.EQ.2) WRITE(3,336)
  336 FORMAT(/20X,'N E G A T I V E   P A R I T Y')
      WRITE(3,334)
C
C     TRANSFORMATION FROM THE 0.1 MEV TO 1.0 MEV ENERGY BIN
C
      DO 302 J=1,30
      DO 302 I=1,40
      If(J.Eq.1) then
       EOM(I,K)=0
       ERO(I,K)=0
       SIG(I,K)=0
      Endif
      IK=I-1
      ASA=0.0
      ASA1=0.0      
      DO 303 ILI=1,10
      ASA1=ASA1+DenGas(IK*10+ILI,J,K)      
  303 ASA=ASA+DenBCS(IK*10+ILI,J,K)
c      DenGas(I,J,K)=ASA1*0.1  ! CORRECTION 03/97
c  302 DenBCS(I,J,K)=ASA *0.1  !
      DenGas(I,J,K)=ASA1
  302 DenBCS(I,J,K)=ASA
C
C    STATE DENSITY CALCULATION
C     
      NC=0
      DO 174 I=1,KMAX
      EOM(I,K)=0.
      DO 174 J=1,30                                     
      IF(DenBCS(I,J,K).Gt.0) THEN
       EOM(I,K)=EOM(I,K)+DenBCS(I,J,K)
       NC=NC+DenBCS(I,J,K)
      ENDIF            
  174 CONTINUE             
      NTOT=NTOT+NC
C
C    MAXIMUM ENERGY OF THE CALCULATIONS
C 
      IF(K.EQ.1) KMAX2=40
      IF(K.EQ.2) KMAX2=MIN(KMAX2,40)
      DO 130 J=KMAX,1,-1
        IF(EOM(J,K).EQ.0.) goto 130
        IF(K.EQ.1) KMAX2=J
        IF(K.EQ.2) KMAX2=MAX(KMAX2,J)
        goto 135
130   continue
135   continue
  337  FORMAT(/20X,'NO STATES WITH THIS PARITY'/)    
      
      IF(NC.EQ.0) THEN
       WRITE(3,337) 
       GOTO 333
      ENDIF
      
      IF(K.NE.2) GO TO 461
      
C
C  CALCULATION OF THE +/ALL PARITY RATIO
C                         
      DO 462 I=1,KMAX2
      PARM(I)=100.0
      EOP=EOM(I,1)
C     IF(EOP.EQ.1.0) EOP=0.0
      EOMS=EOP+EOM(I,2)
      IF(EOMS.EQ.0.0) GO TO 462
      PARM(I)=EOP/EOMS
  462 CONTINUE
  461 CONTINUE
      DO 401 I=1,30
  401 SPIN(I)=I-SS
      DO 402 I=1,KMAX
      SIG(I,K)=0.0
  402 ENERG(I)=I-0.5
      DO 503 I=1,KMAX
      IF(EOM(I,K).EQ.0.0) GO TO 503
      
      SUM=0.0
C
C   CALCULATION OF THE SPIN CUT-OFF FACTORS
C
      DO 504 J=1,30
      IF(DenBCS(I,J,K).GE.0.0) GO TO 519
      SIG(I,K)=0.0
      GO TO 503
  519 SUM=SUM+SPIN(J)**2*DenBCS(I,J,K)
  504 CONTINUE     
      IF(EOM(I,K).GT.0.) THEN 
      SIG(I,K)=SUM/EOM(I,K)
      SIG(I,K)=SQRT(SIG(I,K))
      ELSE
      SIG(I,K)=0
      ENDIF
  503 CONTINUE
C
C    S T A T E    D E N S I T Y     P R I N T O U T
C
      WRITE(3,397)
 397  FORMAT(3X/30X,' ENERGY AND M DEPENDENCE OF STATE DENSITIES'//
     X3X,1HE,4X,'OMEGA(E,PI)  ','SIGMA',18X,' OMEGA(E,M,PI)')
      WRITE(3,94) (ENERG(I),EOM(I,K),SIG(I,K),(DenBCS(I,J,K),J=1,15),I=1
     X,KMAX2)
      WRITE(3,342) ((DenBCS(I,J,K),J=16,30),I=1,KMAX2)
      IF(DEFOR) GO TO 214
C
C   CALCULATION OF THE LEVEL DENSITY FROM STATE DENSITY
C
      DO 93 J=1,KMAX2
      ERO(J,K)=0.0
      IF(SS.EQ.1.0) DenBCS(J,1,K)=DenBCS(J,1,K)*2.
      DO 93 I=1,29                                                   
      DenBCS(J,I,K)=(DenBCS(J,I,K)-DenBCS(J,I+1,K))/2.0
      If(DenBCS(J,I,K).LT.0) then
          Neg_Lev=Neg_Lev+1
          DenBCS(J,I,K)=0
      EndIf
  93  ERO(J,K)=ERO(J,K)+DenBCS(J,I,K)  
  
C
C    L E V E L    D E N S I T Y     P R I N T O U T
C
      WRITE(3,398)
 398  FORMAT(3X//30X,'ENERGY AND SPIN DEPENDENCE OF LEVEL DENSITIES'//
     X3X,1HE,4X,'  RHO(E,PI)  ','SIGMA',18X,          '   RHO(E,J,PI)')
      WRITE(3,94) (ENERG(I),ERO(I,K),SIG(I,K),(DenBCS(I,J,K),J=1,15)
     X,I=1,KMAX2)
  94  FORMAT(F5.1,1X,E13.6,F6.2,1X,15I7)
      WRITE(3,342) ((DenBCS(I,J,K),J=16,30),I=1,KMAX2)
  342 FORMAT(1X,/,(26X,15I7))
      WRITE(3,345) Neg_Lev
  345 FORMAT(/5X,' NUMBER OF LEVELS WITH NEGATIVE DENSITY=',I8/)   
  214 IF(K.NE.2) GO TO 463
C
C   +/ALL PARITY RATIO PRINTOUT
C
      WRITE(3,334)
      WRITE(3,471) (ENERG(I),PARM(I),I=1,KMAX2)
  471 FORMAT(1X,'POS. PARITY STATES/ALL STATES',//,(5X,F4.1,2X,F5.3))
  463 CONTINUE
C
C  SPIN CUT-0FF FACTOR ENERGY DEPENDENCE ANALYSIS
C
      NN=0
      ENES2=0.0
      ENES=0.0
      SIGS=0.0
      SMIX=0.0
      DO 756 I=1,KMAX2
      IF(SIG(I,K).EQ.0.0) GO TO 756
      NN=NN+1
      ENES2=ENES2+ENERG(I)**2
      ENES=ENES+ENERG(I)
      SIGS=SIGS+SIG(I,K)
  756 SMIX=SMIX+SIG(I,K)*ENERG(I)
      IF(NN.GT.1) THEN
       SAVE(K)=SIGS/NN
       EAVE=ENES/NN
       STDE=ENES2/NN-EAVE**2
       IF(STDE.GT.0.0) THEN
        SLOPE(K)=SMIX/NN-EAVE*SAVE(K)
        SLOPE(K)=SLOPE(K)/STDE
        TRC(K)=SAVE(K)-SLOPE(K)*EAVE
        SAVEN(K)=SAVE(K)**2/(NP+NH)
        SAVENA(K)=SAVEN(K)/(FLOAT(Ja)**0.666666)
        WRITE(3,757)
        WRITE(3,758) SAVE(K),SAVEN(K),SAVENA(K),TRC(K),SLOPE(K)
  757   FORMAT(1X,//,10X,'SPIN CUT-OFF ENERGY DEPENDENCE ANALYSIS')      
  758   FORMAT(1X,/,10X,'SIG AVE.=',F6.3,' SIG**2/N=',
     X F6.3,' SIG**2/N/A**2/3=',F6.3,/,10X,'SIG**2=',F6.3,'+',F6.3,'*E')
       ENDIF
      ENDIF
  333 CONTINUE
      If(Kind.eq.1) then
        WRITE(3, 20 ) SMat(Ja-Jp),Ja,NLV,NLV
      Else
        WRITE(3, 20 ) SMat(Jp),Ja,NLV,NLV
      EndIf
      WRITE(3,117)GI/Ja,GI
      
      IF(NTOT.GT.0) THEN
      
      IPP=KMAX2
      IF(IPP.LT.IE) IPP=IE
      IF(IPP.GT.KMAX) IPP=KMAX
      DO 565 I=IE,100
      SHM(I)=0.0
  565 ERS(I)=0.0
  
      WRITE(3,467)
      WRITE(3,466) (ENERG(I),EOM(I,1),EOM(I,2),ERO(I,1),
     X  ERO(I,2),SIG(I,1),SIG(I,2),PARM(I),ERS(I),SHM(I),I=1,IPP)
      Open(8,File=FileName//'.dat')
      Open(18,File=FileName//'.nav')
      NTOT=0
      NTOTL=0
      Do 170 I=1,Kmax2
C      D3 = STATE DENSITY 
C      D4 = LEVEL DENSITY
       D3=EOM(I,1)+EOM(I,2)                                               
       D4=ERO(I,1)+ERO(I,2)                                               
       NTOT=NTOT+NINT(D3)
       NTOTL=NTOTL+NINT(D4)
       If(D3.eq.0.D0)D3=0.1D0                                             
       If(D4.eq.0.D0)D4=0.1D0                                             
       Write(8,1001) (I-1),D4,D3,MAX(REAL(NTOT),1.)
       Write(8,1001) I,D4,D3,MAX(REAL(NTOT),1.)
170   continue
C***************************************************
      close(8)          
      ELSE
       NTOTL=0
      ENDIF
      
      WRITE(3,1003) NTOT,NTOTL
      WRITE(6,1003) NTOT,NTOTL
      
54321 WRITE(3,434)
  434 FORMAT(1X,//
     > 10x,'************************************************'/
     > 10x,'************* FERMI GAS RESULTS ****************'/
     > 10x,'************************************************')
     
      DO 433 K=1,2
      Neg_Lev=0
      IF(K.EQ.1) WRITE(3,335)
      IF(K.EQ.2) WRITE(3,336)
      WRITE(3,334)
      
      NC=0
      DO 474 I=1,40
      EOM(I,K)=0.
      ERO(I,K)=0.
      DO 474 J=1,30
      IF(DenGas(I,J,K).GT.0) THEN
       EOM(I,K)=EOM(I,K)+DenGas(I,J,K)                                
       NC=NC+DenGas(I,J,K)
      ENDIF 
  474 CONTINUE    
          
      NTOT=NTOT+NC  
      
      If(NTOT.Eq.0.AND.K.Eq.2) 
     >        STOP 'NO STATES at all FOR THESE PARAMETERS '
C
C    MAXIMUM ENERGY OF THE CALCULATIONS
C 
      IF(K.EQ.1) KMAX2=40
      IF(K.EQ.2) KMAX2=MIN(KMAX2,40)
      DO 180 J=KMAX,1,-1
        IF(EOM(J,K).EQ.0.) goto 180
        IF(K.EQ.1) KMAX2=J
        IF(K.EQ.2) KMAX2=MAX(KMAX2,J)
        goto 185
180   continue
185   continue
      
      IF(NC.EQ.0) THEN
       WRITE(3,337) 
       GOTO 433
      ENDIF       
  
      IF(K.NE.2) GO TO 561
      
C
C  CALCULATION OF THE +/ALL PARITY RATIO
C
      DO 562 I=1,KMAX2
      PARM(I)=100.0
      EOP=EOM(I,1)
C     IF(EOP.EQ.1.0) EOP=0.0
      EOMS=EOP+EOM(I,2)
      IF(EOMS.EQ.0.0) GO TO 562
      PARM(I)=EOP/EOMS
  562 CONTINUE
  561 CONTINUE
      DO 601 I=1,30
  601 SPIN(I)=I-SS
      DO 602 I=1,40
      SIG(I,K)=0.0
  602 ENERG(I)=I-0.5
      DO 603 I=1,40
      IF(EOM(I,K).EQ.0.0) GO TO 603
      SUM=0.0
C
C   CALCULATION OF THE SPIN CUT-OFF FACTORS
C
      DO 604 J=1,30
      IF(DenGas(I,J,K).GE.0.0) GO TO 619
      SIG(I,K)=0.0
      GO TO 603
  619 SUM=SUM+SPIN(J)**2*DenGas(I,J,K)
  604 CONTINUE
      IF(EOM(I,K).GT.0.) THEN 
      SIG(I,K)=SUM/EOM(I,K)
      SIG(I,K)=SQRT(SIG(I,K))
      ELSE
      SIG(I,K)=0
      ENDIF
  603 CONTINUE
C
C    S T A T E    D E N S I T Y     P R I N T O U T
C
      WRITE(3,397)
      WRITE(3,94) (ENERG(I),EOM(I,K),SIG(I,K),(DenGas(I,J,K),J=1,15),I=1
     X,KMAX2)
      WRITE(3,342) ((DenGas(I,J,K),J=16,30),I=1,KMAX2)
      IF(DEFOR) GO TO 314
C
C   CALCULATION OF THE LEVEL DENSITY FROM STATE DENSITY
C            
      DO 193 J=1,KMAX2
      ERO(J,K)=0.0
      IF(SS.EQ.1.0) DenGas(J,1,K)=DenGas(J,1,K)*2.
      DO 193 I=1,29         
      DenGas(J,I,K)=(DenGas(J,I,K)-DenGas(J,I+1,K))/2.0
      If(DenGas(J,I,K).LT.0) then
          Neg_Lev=Neg_Lev+1
          DenGas(J,I,K)=0
      EndIf
 193  ERO(J,K)=ERO(J,K)+DenGas(J,I,K)
  
C
C    L E V E L    D E N S I T Y     P R I N T O U T
C
      WRITE(3,398)
      WRITE(3,94) (ENERG(I),ERO(I,K),SIG(I,K),(DenGas(I,J,K),J=1,15)
     X,I=1,KMAX2)
      WRITE(3,342) ((DenGas(I,J,K),J=16,30),I=1,KMAX2)
      WRITE(3,345) Neg_Lev
  314 IF(K.NE.2) GO TO 563
C
C   +/ALL PARITY RATIO PRINTOUT
C
      WRITE(3,334)
      WRITE(3,471) (ENERG(I),PARM(I),I=1,KMAX2)
  563 CONTINUE
C
C  SPIN CUT-0FF FACTOR ENERGY DEPENDENCE ANALYSIS
C
      NN=0
      ENES2=0.0
      ENES=0.0
      SIGS=0.0
      SMIX=0.0
      DO 856 I=1,KMAX2
      IF(SIG(I,K).EQ.0.0) GO TO 856
      NN=NN+1
      ENES2=ENES2+ENERG(I)**2
      ENES=ENES+ENERG(I)
      SIGS=SIGS+SIG(I,K)
  856 SMIX=SMIX+SIG(I,K)*ENERG(I)
      IF(NN.GT.1) THEN
       SAVE(K)=SIGS/NN
       EAVE=ENES/NN
       STDE=ENES2/NN-EAVE**2
       IF(STDE.GT.0.) THEN
        SLOPE(K)=SMIX/NN-EAVE*SAVE(K)
        SLOPE(K)=SLOPE(K)/STDE
        TRC(K)=SAVE(K)-SLOPE(K)*EAVE
        SAVEN(K)=SAVE(K)**2/(NP+NH)
        SAVENA(K)=SAVEN(K)/(FLOAT(Ja)**0.666666)
        WRITE(3,757)
        WRITE(3,758) SAVE(K),SAVEN(K),SAVENA(K),TRC(K),SLOPE(K)
       ENDIF
      ENDIF
  433 CONTINUE
      Open(8,File=FileName//'.fer')
      NTOT=0
      NTOTL=0
      Do 190 I=1,Kmax2
C      D3 = STATE DENSITY 
C      D4 = LEVEL DENSITY
       D3=EOM(I,1)+EOM(I,2)                                               
       D4=ERO(I,1)+ERO(I,2)                                               
       NTOT=NTOT+NINT(D3)
       NTOTL=NTOTL+NINT(D4)
       If(D3.eq.0.D0)D3=0.1D0                                             
       If(D4.eq.0.D0)D4=0.1D0                                             
       Write(8,1001) (I-1),D4,D3
       Write(8,1001) I,D4,D3
190   continue
C***************************************************
      close(8)     
      
      WRITE(3,1004) NTOT,NTOTL
      WRITE(6,1004) NTOT,NTOTL
      
      CLOSE(6)
      If(Kind.eq.1) then
        WRITE(3, 20 ) SMat(Ja-Jp),Ja,NLV,NLV
      Else
        WRITE(3, 20 ) SMat(Jp),Ja,NLV,NLV
      EndIf
      WRITE(3,117)GI/Ja,GI
      
      IPP=KMAX2
      IF(IPP.LT.IE) IPP=IE
      IF(IPP.GT.40) IPP=40
      DO 465 I=IE,100
      SHM(I)=0.0
  465 ERS(I)=0.0
  
      WRITE(3,467)
  467 FORMAT(2X,'MEV ',3X,
     X'+ST. DENS.',3X,'-ST. DENS.',3X,'+LEV. DENS.',2X,'-LEV. DENS.',
     X1X,'+SIG.',1X,'-SIG.',1X,'+/ALL  ','  MEV ','   SHIFT',/)
      WRITE(3,466) (ENERG(I),EOM(I,1),EOM(I,2),ERO(I,1),
     XERO(I,2),SIG(I,1),SIG(I,2),PARM(I),ERS(I),SHM(I),I=1,IPP)
  466 FORMAT(1X,F4.1,1X,4E13.5,2F6.2,2X,F5.3,2X,F5.2,I9.0)
      WRITE(3,758) SAVE(1),SAVEN(1),SAVENA(1),TRC(1),SLOPE(1)
      WRITE(3,758) SAVE(2),SAVEN(2),SAVENA(2),TRC(2),SLOPE(2)
 1003 FORMAT(/1X,'TOTAL NUMBER OF STATES(BCS) = ',I10/
     X        1X,'TOTAL NUMBER OF LEVELS(BCS) = ',I10)
 1004 FORMAT(/1X,'TOTAL NUMBER OF STATES(Fermi Gas) = ',I10/
     X        1X,'TOTAL NUMBER OF LEVELS(Fermi Gas) = ',I10)
   20 FORMAT(//25X,' NUCLEUS ',A2,I3/,25X,15(1H-),//10X,I3,
     X' NEUTRON OR',I4,' PROTON LEVELS TAKEN INTO ACCOUNT')
  117 FORMAT(/10X,'G=',F10.5,' CORRESPONDS TO',F6.2,'/A'//)
   24 FORMAT (10X,I1,' PARTICLES',I4,' HOLES CONFIGURATIONS',//)
 1000 Format(6(1x,E12.5))                                                         
 1001 Format(1x,I4,4(1x,F14.2))
      Return
      End
      Function Smat(I)
      Character*2 Mat(99),Smat
      DATA MAT/' H','He','Li','Be',' B',' C',' N',' O',' F','Ne','Na',
     >         'Mg','Al','Si',' P',' S','Cl','Ar',' K','Ca','Sc','Ti',
     >         ' V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As',
     >         'Se','Br','Kr','Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru',
     >         'Rh','Pd','Ag','Cd','In','Sn','Sb','Te',' I','Xe','Cs',
     >         'Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy',
     >         'Ho','Er','Tm','Yb','Lu','Hf','Ta',' W','Re','Os','Ir',
     >         'Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra',
     >         'Ac','Th','Pa',' U','Np','Pu','Am','Cm','Bk','Cf','Es'/
      if(i.lt.100)Then
       Smat= Mat(I)
      else
       Smat= 'NO'
      endif
      Return
      End
      SUBROUTINE PNXCB(M,K,N,MOLD,MNEW)      
      Implicit Logical*1(A-Z)
      INTEGER*4 M(*),K,N,MOLD,MNEW,J
C-----------------------------------------------------------------------
C-----FROM:
C-----COMBINATION GENERATION
C-----W.H.PAYNE,F.M.IVES - ACM TRANS.OF MATH.SOFT. 5,2(1979) P.163
C-----------------------------------------------------------------------
C LIU-TANG SEQUENCE POINTER PROGRAMMED COMBINATION GENERATOR.
C EXCHANGE MARK LOCATED AT MOLD WITH MARK LOCATED AT MNEW.
C N MARKS TAKEN K AT A TIME. N>1. 0<K<=N. SEPTEMBER 1975.
C M IS A POINTER ARRAY OF DIMENSION MIN(K,N-K) POINTING TO THE
C POSITIONS OF THE MIN(K,N-K) IDENTICAL MARKS.
      J=1
C FOR 'DOWN-UP' SEQUENCE INSERT GO TO 20 HERE.
C     J-TH MARK IS MOVING RIGHT.
10    MOLD=M(J)
      MNEW=MOLD+1
      IF(J.NE.K)GO TO 11
      IF(MNEW.GT.N)GO TO 31
      GO TO 32
11    IF(MNEW.LT.M(J+1))GO TO 32
C     J-TH MARK CAN'T MOVE. CONSIDER NEXT MARK.
      J=J+1
C     J-TH MARK IS MOVING LEFT
20    MOLD=M(J)
      MNEW=MOLD-1
      IF(MNEW.GE.J)GO TO 42
      IF(J.EQ.K)GO TO 41
C     J-TH MARK CAN'T MOVE. CONSIDER NEXT MARK.
      J=J+1
      GO TO 10
C     MOVE MARK J.
C     PREVIOUS SECTION BRANCHED TO 31,32,41,OR 42.
C     DEPENDING ON TYPE OF MOVE.
C     31 MOVING RIGHT, J=K,AT END OF SEQUENCE
C     32 MOVING RIGHT, NOT AT END OF SEQUENCE
31    MNEW=K
      M(J)=MNEW
      RETURN
32    M(J)=MNEW
      IF(J.EQ.1)RETURN
      M(J-1)=MOLD
      MOLD=J-1
      RETURN
C     41 MOVING LEFT, J=K, AT END OF SEQUENCE
C     42 MOVING RIGHT, NOT AT END OF SEQUENCE
41    MNEW=N
      M(J)=MNEW
      RETURN
42    M(J)=MNEW
      IF(J.EQ.1)RETURN
      MNEW=J-1
      M(J-1)=MNEW
      RETURN
      END