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