hdiff output

r31608/Cv.BS2.f90 2016-12-09 15:30:08.721664234 +0000 r31607/Cv.BS2.f90 2016-12-09 15:30:08.973667478 +0000
  1: !  1: !
  2: !  Calculate Cv curve from weights and visit statistics  2: !  Calculate Cv curve from weights and visit statistics
  3: !  for replicas at different temperatures REPT(J1).  3: !  for replicas at different temperatures REPT(J1).
  4: !  Minimise chi^2 statistic to extract best probabilities.  4: !  Minimise chi^2 statistic to extract best probabilities.
  5: !  5: !
  6: PROGRAM CVBS  6: PROGRAM CVBS
  7: IMPLICIT NONE  7: IMPLICIT NONE
  8: DOUBLE PRECISION TMIN, TMAX, TINT, Z0, Z1, Z2, TEMPERATURE, DUMMY, DELTA, ONEMEXP, PSUM, DP, PSUM2, TANAL, DPLUS, DMINUS  8: DOUBLE PRECISION TMIN, TMAX, TINT, Z0, Z1, Z2, TEMPERATURE, DUMMY, DELTA, ONEMEXP, PSUM, DP, PSUM2, TANAL, DPLUS, DMINUS
  9: DOUBLE PRECISION COLOURTHRESH, HSHIFT, DUMMY2, DWEIGHT, DMIN, DUMMY3, DUMMY4, LN2DSUM, PEMIN, LN1DSUM, DQPLUS, DQMINUS  9: DOUBLE PRECISION COLOURTHRESH, HSHIFT, DUMMY2, DWEIGHT, DMIN, DUMMY3, DUMMY4, LN2DSUM
 10: DOUBLE PRECISION, ALLOCATABLE :: ENERGY(:), LNWEIGHT(:), MEANE(:), ZSAVE(:), PERMINP(:), PERMINM(:), EMIN(:), LNW2D(:,:) 10: DOUBLE PRECISION, ALLOCATABLE :: ENERGY(:), LNWEIGHT(:), MEANE(:), ZSAVE(:), PERMINP(:), PERMINM(:), EMIN(:), LNW2D(:,:)
 11: DOUBLE PRECISION, ALLOCATABLE :: QENERGY(:), CQP(:), CQM(:), CHECK1D(:) 11: DOUBLE PRECISION, ALLOCATABLE :: QENERGY(:), CQP(:), CQM(:)
 12: INTEGER, ALLOCATABLE :: CLOSEST(:), QBIN(:) 12: INTEGER, ALLOCATABLE :: CLOSEST(:), QBIN(:)
 13: INTEGER, ALLOCATABLE :: MINP(:), MINM(:) 13: INTEGER, ALLOCATABLE :: MINP(:), MINM(:)
 14: DOUBLE PRECISION, ALLOCATABLE :: PERMINQP(:), PERMINQM(:) 14: INTEGER J1, NTEMP, J2, KAPPA, NPLUS, NMINUS, NDUMMY, NMIN, MAXBIN, NQBINS, NCOUNT
 15: INTEGER, ALLOCATABLE :: MINQP(:),MINQM(:) 
 16: INTEGER J1, NTEMP, J2, KAPPA, NPLUS, NMINUS, NDUMMY, NMIN, MAXBIN, NQBINS, NQPLUS, NQMINUS 
 17: LOGICAL YESNO, DO2DT 15: LOGICAL YESNO, DO2DT
 18: LOGICAL, ALLOCATABLE :: NO1D(:),NO2D(:,:) 16: LOGICAL, ALLOCATABLE :: NO1D(:),NO2D(:,:)
 19:  17: 
 20: OPEN(UNIT=10,FILE='Cv.BS.data',STATUS='OLD') 18: OPEN(UNIT=10,FILE='Cv.BS.data',STATUS='OLD')
 21: READ(10,*) TMIN, TMAX, NTEMP, KAPPA, NMIN, MAXBIN 19: READ(10,*) TMIN, TMAX, NTEMP, KAPPA, NMIN, MAXBIN
 22: CLOSE(10) 20: CLOSE(10)
 23: ALLOCATE(MEANE(NTEMP),ZSAVE(NTEMP)) 21: ALLOCATE(MEANE(NTEMP),ZSAVE(NTEMP))
 24:  22: 
 25: PRINT '(3(A,I10))','bins=',MAXBIN 23: PRINT '(3(A,I10))','bins=',MAXBIN
 26: PRINT '(A,2G15.5,A,I8)','minimum and maximum T for Cv calculation=',TMIN,TMAX,' number of Cv data points=',NTEMP 24: PRINT '(A,2G15.5,A,I8)','minimum and maximum T for Cv calculation=',TMIN,TMAX,' number of Cv data points=',NTEMP
 30:  28: 
 31: OPEN(UNIT=10,FILE='min.data',STATUS='OLD') 29: OPEN(UNIT=10,FILE='min.data',STATUS='OLD')
 32: DO J1=1,NMIN 30: DO J1=1,NMIN
 33:    READ(10,*) EMIN(J1) 31:    READ(10,*) EMIN(J1)
 34: !  PRINT '(A,I8,G20.10)','J1,emin=',J1,EMIN(J1) 32: !  PRINT '(A,I8,G20.10)','J1,emin=',J1,EMIN(J1)
 35: ENDDO 33: ENDDO
 36: CLOSE(10) 34: CLOSE(10)
 37:  35: 
 38: OPEN(UNIT=10,FILE='weights.BS',STATUS='OLD') 36: OPEN(UNIT=10,FILE='weights.BS',STATUS='OLD')
 39: NO1D(1:MAXBIN)=.TRUE. 37: NO1D(1:MAXBIN)=.TRUE.
 40: PEMIN=HUGE(1.0D0) 
 41: DELTA=0.0D0 
 42: LN1DSUM=0.0D0 
 43: DO J1=1,MAXBIN 38: DO J1=1,MAXBIN
 44:    READ(10,*,END=432) NDUMMY,DUMMY,DUMMY2,DUMMY3,DUMMY4 39:    READ(10,*,END=432) NDUMMY,DUMMY,DUMMY2,DUMMY3,DUMMY4
 45:    ENERGY(NDUMMY)=DUMMY2 40:    ENERGY(NDUMMY)=DUMMY2
 46:    LNWEIGHT(NDUMMY)=DUMMY3 41:    LNWEIGHT(NDUMMY)=DUMMY4
 47:    LN1DSUM=LN1DSUM+EXP(DUMMY3) 
 48:    NO1D(NDUMMY)=.FALSE. 42:    NO1D(NDUMMY)=.FALSE.
 49:    IF (ENERGY(NDUMMY).LT.PEMIN) THEN 
 50:       PEMIN=ENERGY(NDUMMY) 
 51:    ENDIF 
 52:    IF (NDUMMY-1.GT.0) THEN 
 53:       IF ((.NOT.NO1D(NDUMMY)).AND.(.NOT.NO1D(NDUMMY-1))) DELTA=ENERGY(NDUMMY)-ENERGY(NDUMMY-1) 
 54:    ENDIF 
 55: ENDDO 43: ENDDO
 56: 432 CONTINUE 44: 432 CONTINUE
 57: CLOSE(10) 45: CLOSE(10)
 58: LN1DSUM=LOG(LN1DSUM) 
 59:  46: 
 60: TINT=(TMAX-TMIN)/(NTEMP-1) 47: TINT=(TMAX-TMIN)/(NTEMP-1)
  48: DELTA=ENERGY(2)-ENERGY(1)
 61: ! 49: !
 62: !  Calculate Z0, Z1 and Z2 over the required T range. Omit factors of (kT/h)^kappa, 50: !  Calculate Z0, Z1 and Z2 over the required T range. Omit factors of (kT/h)^kappa,
 63: !  which cancel. 51: !  which cancel.
 64: ! 52: !
 65: OPEN(UNIT=1,FILE='Cv.out.BS',STATUS='UNKNOWN') 53: OPEN(UNIT=1,FILE='Cv.out.BS',STATUS='UNKNOWN')
 66: DO J1=1,NTEMP 54: DO J1=1,NTEMP
 67:    Z0=0.0D0 55:    Z0=0.0D0
 68:    Z1=0.0D0 56:    Z1=0.0D0
 69:    Z2=0.0D0 57:    Z2=0.0D0
 70:    TEMPERATURE=TMIN+(J1-1)*TINT 58:    TEMPERATURE=TMIN+(J1-1)*TINT
 71:    DO J2=1,MAXBIN 59:    DO J2=1,MAXBIN
 72:       IF (NO1D(J2)) CYCLE 60:       IF (NO1D(J2)) CYCLE
 73:       DUMMY=EXP(-(ENERGY(J2)-PEMIN)/TEMPERATURE+LNWEIGHT(J2)) 61:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT(J2))
 74:       Z0=Z0+DUMMY 62:       Z0=Z0+DUMMY
 75:       Z1=Z1+DUMMY*(ENERGY(J2)-PEMIN) 63:       Z1=Z1+DUMMY*(ENERGY(J2)-ENERGY(1))
 76:       Z2=Z2+DUMMY*(ENERGY(J2)-PEMIN)**2 64:       Z2=Z2+DUMMY*(ENERGY(J2)-ENERGY(1))**2
 77:    ENDDO 65:    ENDDO
 78:    MEANE(J1)=Z1/Z0 66:    MEANE(J1)=Z1/Z0
 79:    ZSAVE(J1)=Z0 67:    ZSAVE(J1)=Z0
 80:    IF (DELTA/TEMPERATURE.LT.1.0D-7) THEN 68:    IF (DELTA/TEMPERATURE.LT.1.0D-7) THEN
 81:       ONEMEXP=-DELTA/TEMPERATURE 69:       ONEMEXP=-DELTA/TEMPERATURE
 82:    ELSE 70:    ELSE
 83:       ONEMEXP= 1.0D0-EXP(DELTA/TEMPERATURE) 71:       ONEMEXP= 1.0D0-EXP(DELTA/TEMPERATURE)
 84:    ENDIF 72:    ENDIF
 85:    WRITE(1,'(7G20.10)') TEMPERATURE, Z0, Z1, Z2, & 73:    WRITE(1,'(7G20.10)') TEMPERATURE, Z0, Z1, Z2, &
 86:   &                     KAPPA/2.0D0 + & 74:   &                     KAPPA/2.0D0 + &
 89: ENDDO 77: ENDDO
 90: CLOSE(1) 78: CLOSE(1)
 91:  79: 
 92: OPEN(UNIT=1,FILE='prob.out.BS',STATUS='UNKNOWN') 80: OPEN(UNIT=1,FILE='prob.out.BS',STATUS='UNKNOWN')
 93: DO J1=1,NTEMP 81: DO J1=1,NTEMP
 94:    TEMPERATURE=TMIN+(J1-1)*TINT 82:    TEMPERATURE=TMIN+(J1-1)*TINT
 95:    PSUM=0.0D0 83:    PSUM=0.0D0
 96:    PSUM2=0.0D0 84:    PSUM2=0.0D0
 97:    DO J2=1,MAXBIN 85:    DO J2=1,MAXBIN
 98:       IF (NO1D(J2)) CYCLE 86:       IF (NO1D(J2)) CYCLE
 99:       DUMMY=EXP(-(ENERGY(J2)-PEMIN)/TEMPERATURE+LNWEIGHT(J2)) 87:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT(J2))
100:       DP=DUMMY*(ENERGY(J2)-PEMIN-MEANE(J1))/(ZSAVE(J1)*TEMPERATURE**2) 88:       DP=DUMMY*(ENERGY(J2)-ENERGY(1)-MEANE(J1))/(ZSAVE(J1)*TEMPERATURE**2)
101:       PSUM=PSUM+DP*(ENERGY(J2)-PEMIN-MEANE(J1)) 89:       PSUM=PSUM+DP*(ENERGY(J2)-ENERGY(1)-MEANE(J1))
102: !     PRINT '(A,2I6,4G20.10)','J1,J2,T,ZSAVE,DUMMY,DP=',J1,J2,TEMPERATURE,ZSAVE(J1),DUMMY,DP 90: !     PRINT '(A,2I6,4G20.10)','J1,J2,T,ZSAVE,DUMMY,DP=',J1,J2,TEMPERATURE,ZSAVE(J1),DUMMY,DP
103:       PSUM2=PSUM2+ZSAVE(J1)*(TEMPERATURE*DP)**2/DUMMY 91:       PSUM2=PSUM2+ZSAVE(J1)*(TEMPERATURE*DP)**2/DUMMY
104: !     PSUM2=PSUM2+(ZSAVE(J1)/DUMMY)*(TEMPERATURE*DP)**2 92: !     PSUM2=PSUM2+(ZSAVE(J1)/DUMMY)*(TEMPERATURE*DP)**2
105:    ENDDO 93:    ENDDO
106:    IF (DELTA/TEMPERATURE.LT.1.0D-7) THEN 94:    IF (DELTA/TEMPERATURE.LT.1.0D-7) THEN
107:       ONEMEXP=-DELTA/TEMPERATURE 95:       ONEMEXP=-DELTA/TEMPERATURE
108:    ELSE 96:    ELSE
109:       ONEMEXP= 1.0D0-EXP(DELTA/TEMPERATURE) 97:       ONEMEXP= 1.0D0-EXP(DELTA/TEMPERATURE)
110:    ENDIF 98:    ENDIF
111:    WRITE(1,'(5G20.10)') TEMPERATURE, KAPPA/2.0D0+1.0D0 - DELTA**2*EXP(DELTA/TEMPERATURE)/(ONEMEXP**2*TEMPERATURE**2)+PSUM, & 99:    WRITE(1,'(5G20.10)') TEMPERATURE, KAPPA/2.0D0+1.0D0 - DELTA**2*EXP(DELTA/TEMPERATURE)/(ONEMEXP**2*TEMPERATURE**2)+PSUM, &
126:    COLOURTHRESH=0.9D0114:    COLOURTHRESH=0.9D0
127:    READ(1,*,END=666) COLOURTHRESH115:    READ(1,*,END=666) COLOURTHRESH
128: 666 CONTINUE116: 666 CONTINUE
129:    PRINT '(2(A,G20.10))','Analysing heat capacity contributions at T=',TANAL,' colour threshold=',COLOURTHRESH117:    PRINT '(2(A,G20.10))','Analysing heat capacity contributions at T=',TANAL,' colour threshold=',COLOURTHRESH
130:    CLOSE(1)118:    CLOSE(1)
131:    TEMPERATURE=TANAL119:    TEMPERATURE=TANAL
132:    Z0=0.0D0120:    Z0=0.0D0
133:    Z1=0.0D0121:    Z1=0.0D0
134:    DO J2=1,MAXBIN122:    DO J2=1,MAXBIN
135:       IF (NO1D(J2)) CYCLE123:       IF (NO1D(J2)) CYCLE
136:       DUMMY=EXP(-(ENERGY(J2)-PEMIN)/TEMPERATURE+LNWEIGHT(J2))124:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT(J2))
137:       Z0=Z0+DUMMY125:       Z0=Z0+DUMMY
138:       Z1=Z1+DUMMY*(ENERGY(J2)-PEMIN)126:       Z1=Z1+DUMMY*(ENERGY(J2)-ENERGY(1))
139:    ENDDO127:    ENDDO
140:    MEANE(1)=Z1/Z0128:    MEANE(1)=Z1/Z0
141:    ZSAVE(1)=Z0129:    ZSAVE(1)=Z0
142:    DO J2=1,MAXBIN130:    DO J2=1,MAXBIN
143:       IF (NO1D(J2)) CYCLE131:       IF (NO1D(J2)) CYCLE
144:       DUMMY=EXP(-(ENERGY(J2)-PEMIN)/TEMPERATURE+LNWEIGHT(J2))132:       DUMMY=EXP(-(ENERGY(J2)-ENERGY(1))/TEMPERATURE+LNWEIGHT(J2))
145:       DP=DUMMY*(ENERGY(J2)-PEMIN-MEANE(1))/(ZSAVE(1)*TEMPERATURE**2)133:       DP=DUMMY*(ENERGY(J2)-ENERGY(1)-MEANE(1))/(ZSAVE(1)*TEMPERATURE**2)
146:       DUMMY=DP*(ENERGY(J2)-PEMIN-MEANE(1))134:       DUMMY=DP*(ENERGY(J2)-ENERGY(1)-MEANE(1))
147:       IF (DP.LT.0.0D0) THEN135:       IF (DP.LT.0.0D0) THEN
148:          NMINUS=NMINUS+1136:          NMINUS=NMINUS+1
149:          MINM(NMINUS)=J2137:          MINM(NMINUS)=J2
150:          PERMINM(NMINUS)=DUMMY138:          PERMINM(NMINUS)=DUMMY
151:          DMINUS=DMINUS+DUMMY139:          DMINUS=DMINUS+DUMMY
152:       ELSE140:       ELSE
153:          NPLUS=NPLUS+1141:          NPLUS=NPLUS+1
154:          MINP(NPLUS)=J2142:          MINP(NPLUS)=J2
155:          PERMINP(NPLUS)=DUMMY143:          PERMINP(NPLUS)=DUMMY
156:          DPLUS=DPLUS+DUMMY144:          DPLUS=DPLUS+DUMMY
157:       ENDIF145:       ENDIF
158: !     PRINT '(A,I6,3G20.10,2I6)','J2,ENERGY(J2)-PEMIN,MEANE,DP,NPLUS,NMINUS=',J2,ENERGY(J2)-PEMIN,MEANE(1),DP,NPLUS,NMINUS146: !     PRINT '(A,I6,3G20.10,2I6)','J2,ENERGY(J2)-ENERGY(1),MEANE,DP,NPLUS,NMINUS=',J2,ENERGY(J2)-ENERGY(1),MEANE(1),DP,NPLUS,NMINUS
159:    ENDDO147:    ENDDO
160:    DO J2=1,NMINUS148:    DO J2=1,NMINUS
161:       PERMINM(J2)=PERMINM(J2)/DMINUS149:       PERMINM(J2)=PERMINM(J2)/DMINUS
162:    ENDDO150:    ENDDO
163:    DO J2=1,NPLUS151:    DO J2=1,NPLUS
164:       PERMINP(J2)=PERMINP(J2)/DPLUS152:       PERMINP(J2)=PERMINP(J2)/DPLUS
165:    ENDDO153:    ENDDO
166:    CALL SORT(NMINUS,MAXBIN,PERMINM,MINM)154:    CALL SORT(NMINUS,MAXBIN,PERMINM,MINM)
167:    CALL SORT(NPLUS,MAXBIN,PERMINP,MINP)155:    CALL SORT(NPLUS,MAXBIN,PERMINP,MINP)
168:  
169:    DUMMY=0.0D0156:    DUMMY=0.0D0
170:    PRINT '(A)','i bins with negative dp/dT below threshold:'157:    PRINT '(A)','i bins with negative dp/dT below threshold:'
171:    DO J2=1,NMINUS158:    DO J2=1,NMINUS
172:       DUMMY=DUMMY+PERMINM(J2)159:       DUMMY=DUMMY+PERMINM(J2)
173:       PRINT '(2I6,4G20.10)',J2,MINM(J2),PERMINM(J2),DUMMY,ENERGY(MINM(J2)),ENERGY(MINM(J2))-PEMIN160:       PRINT '(2I6,4G20.10)',J2,MINM(J2),PERMINM(J2),DUMMY,ENERGY(MINM(J2)),ENERGY(MINM(J2))-ENERGY(1)
174:       IF (DUMMY.GT.COLOURTHRESH) EXIT161:       IF (DUMMY.GT.COLOURTHRESH) EXIT
175:    ENDDO162:    ENDDO
176:    DUMMY=0.0D0163:    DUMMY=0.0D0
177:    PRINT '(A)','i bins with positive dp/dT:'164:    PRINT '(A)','i bins with positive dp/dT:'
178:    DO J2=1,NPLUS165:    DO J2=1,NPLUS
179:       DUMMY=DUMMY+PERMINP(J2)166:       DUMMY=DUMMY+PERMINP(J2)
180:       PRINT '(2I6,4G20.10)',J2,MINP(J2),PERMINP(J2),DUMMY,ENERGY(MINP(J2)),ENERGY(MINP(J2))-PEMIN167:       PRINT '(2I6,4G20.10)',J2,MINP(J2),PERMINP(J2),DUMMY,ENERGY(MINP(J2)),ENERGY(MINP(J2))-ENERGY(1)
181:       IF (DUMMY.GT.COLOURTHRESH) EXIT168:       IF (DUMMY.GT.COLOURTHRESH) EXIT
182:    ENDDO169:    ENDDO
183: 170: 
184:    INQUIRE(FILE='weights.2D.BS',EXIST=DO2DT)171:    INQUIRE(FILE='weights.2D.BS',EXIST=DO2DT)
185:    IF (DO2DT) THEN172:    IF (DO2DT) THEN
186:       PRINT '(A)','Reading 2D weights from weights.2D.BS'173:       PRINT '(A)','Reading 2D weights from weights.2D.BS'
187:       OPEN(UNIT=10,FILE='Cv.BS.data',STATUS='OLD')174:       OPEN(UNIT=10,FILE='Cv.BS.data',STATUS='OLD')
188: !175: !
189: ! This itme, look for NQBINS parameter at the end for allocations176: ! This itme, look for NQBINS parameter at the end for allocations
190: !177: !
191:       READ(10,*) TMIN, TMAX, NTEMP, KAPPA, NMIN, MAXBIN, NQBINS178:       READ(10,*) TMIN, TMAX, NTEMP, KAPPA, NMIN, MAXBIN, NQBINS
192:       CLOSE(10)179:       CLOSE(10)
193:       ALLOCATE(LNW2D(NQBINS,MAXBIN),NO2D(NQBINS,MAXBIN),CQP(NQBINS),CQM(NQBINS),QENERGY(NQBINS),CHECK1D(MAXBIN))180:       ALLOCATE(LNW2D(NQBINS,MAXBIN),NO2D(NQBINS,MAXBIN),CQP(NQBINS),CQM(NQBINS),QENERGY(NQBINS))
194:       LNW2D(1:NQBINS,1:MAXBIN)=-HUGE(1.0D0)181:       LNW2D(1:NQBINS,1:MAXBIN)=-HUGE(1.0D0)
195:       NO2D(1:NQBINS,1:MAXBIN)=.TRUE.182:       NO2D(1:NQBINS,1:MAXBIN)=.TRUE.
196:       CHECK1D(1:MAXBIN)=0.0D0 
197:       OPEN(UNIT=19,FILE='weights.2D.BS',STATUS='OLD')183:       OPEN(UNIT=19,FILE='weights.2D.BS',STATUS='OLD')
 184:       NCOUNT=0
198:       DUMMY=0.0D0185:       DUMMY=0.0D0
199:       DO 186:       DO 
200:          READ(19,*,END=753) J1, J2, DUMMY2, DUMMY3, DWEIGHT187:          READ(19,*,END=753) J1, J2, DUMMY2, DUMMY, DWEIGHT
201:          QENERGY(J1)=DUMMY2188:          QENERGY(J1)=DUMMY2
202:          NO2D(J1,J2)=.FALSE.189:          NO2D(J1,J2)=.FALSE.
203:          LNW2D(J1,J2)=DWEIGHT190:          LNW2D(J1,J2)=DWEIGHT
204:          DUMMY=DUMMY+EXP(DWEIGHT)191:          DUMMY=DUMMY+EXP(DWEIGHT)
205:          CHECK1D(J2)=CHECK1D(J2)+EXP(DWEIGHT) 
206:       ENDDO192:       ENDDO
207: 753   CONTINUE193: 753   CONTINUE
208:       CLOSE(19) 
209:       LN2DSUM=LOG(DUMMY)194:       LN2DSUM=LOG(DUMMY)
210: !     PRINT '(A)','compare 1d i bin weights with sums over 2D q bins'195:       CLOSE(19)
211: !     DUMMY=0.0D0 
212: !     DO J1=1,MAXBIN 
213: !        IF (CHECK1D(J1).GT.0.0D0) CHECK1D(J1)=LOG(CHECK1D(J1)) 
214: !        PRINT '(I8,3G20.10)',J1,CHECK1D(J1),LNWEIGHT(J1),ABS(CHECK1D(J1)-LNWEIGHT(J1)) 
215: !        IF (CHECK1D(J1).GT.0.0D0) DUMMY=DUMMY+EXP(CHECK1D(J1)) 
216: !     ENDDO 
217: !     DUMMY=LOG(DUMMY) 
218: !     PRINT '(A,3G20.10)','LN2DSUM,DUMMY,LN1DSUM=',LN2DSUM,DUMMY,LN1DSUM 
219:  
220:       IF (DELTA/TEMPERATURE.LT.1.0D-7) THEN 
221:          ONEMEXP=-DELTA/TEMPERATURE 
222:       ELSE 
223:          ONEMEXP= 1.0D0-EXP(DELTA/TEMPERATURE) 
224:       ENDIF 
225:       DUMMY=KAPPA/2.0D0 + &  
226:   &                      1.0D0*(1.0D0 - DELTA**2*EXP(DELTA/TEMPERATURE)/(ONEMEXP**2*TEMPERATURE**2)) & 
227:   &                      + DPLUS+DMINUS 
228:  
229:       PRINT '(A,G20.10,A,G20.10)','Looking for total Cv at this T of ',DUMMY,' flux dependent component=',DPLUS+DMINUS 
230: 196: 
231:       DO J2=1,NQBINS197:       DO J2=1,NQBINS
232:          CQP(J2)=0.0D0198:          CQP(J2)=0.0D0
233:          CQM(J2)=0.0D0199:          CQM(J2)=0.0D0
234:          DO J1=1,NPLUS200:          DO J1=1,NPLUS
235:             IF (NO2D(J2,MINP(J1))) CYCLE201:             IF (NO2D(J2,MINP(J1))) CYCLE
236:             IF (NO1D(MINP(J1))) CYCLE202:             IF (NO1D(MINP(J1))) CYCLE
237: !           CQP(J2)=CQP(J2)+PERMINP(J1)*DPLUS*EXP(LNW2D(J2,MINP(J1))-LN2DSUM-LNWEIGHT(MINP(J1))+LN1DSUM)203:             CQP(J2)=CQP(J2)+PERMINP(J1)*DPLUS*EXP(LNW2D(J2,MINP(J1))-LN2DSUM-LNWEIGHT(MINP(J1)))
238:             CQP(J2)=CQP(J2)+PERMINP(J1)*DPLUS*EXP(LNW2D(J2,MINP(J1))-LNWEIGHT(MINP(J1))) 
239: !           PRINT '(A,3I8,6G20.10)','J2,J1,MINP,PERMINP,DPLUS,LNW2D,LN2DSUM,LNWEIGHT,LN1DSUM=', & 
240: ! &                                J2,J1,MINP(J1),PERMINP(J1),DPLUS,LNW2D(J2,MINP(J1)),LN2DSUM,LNWEIGHT(MINP(J1)),LN1DSUM 
241:          ENDDO204:          ENDDO
242:          DO J1=1,NMINUS205:          DO J1=1,NMINUS
243:             IF (NO2D(J2,MINM(J1))) CYCLE206:             IF (NO2D(J2,MINM(J1))) CYCLE
244:             IF (NO1D(MINM(J1))) CYCLE207:             IF (NO1D(MINM(J1))) CYCLE
245: !           CQM(J2)=CQM(J2)+PERMINM(J1)*DMINUS*EXP(LNW2D(J2,MINM(J1))-LN2DSUM-LNWEIGHT(MINM(J1))+LN1DSUM)208:             CQM(J2)=CQM(J2)+PERMINM(J1)*DMINUS*EXP(LNW2D(NCOUNT,MINM(J1))-LN2DSUM-LNWEIGHT(MINM(J1)))
246:             CQM(J2)=CQM(J2)+PERMINM(J1)*DMINUS*EXP(LNW2D(J2,MINM(J1))-LNWEIGHT(MINM(J1))) 
247: !           PRINT '(A,3I8,6G20.10)','J2,J1,MINM,PERMINM,DMINUS,LNW2D,LN2DSUM,LNWEIGHT,LN1DSUM=', & 
248: ! &                                J2,J1,MINM(J1),PERMINM(J1),DMINUS,LNW2D(J2,MINM(J1)),LN2DSUM,LNWEIGHT(MINM(J1)),LN1DSUM 
249:          ENDDO209:          ENDDO
250:       ENDDO210:       ENDDO
251: !211: !
252: ! Breakdown of Cv over q bins at TANAL212: ! Breakdown of Cv over q bins at TANAL
253: !213: !
254:       PRINT '(A,G20.10,A)','Breakdown of Cv over q bins at T=',TEMPERATURE,' q bin, plus, minus, sum, cumulative Cv'214:       PRINT '(A,G20.10,A)','Breakdown of Cv over q bins at T=',TEMPERATURE,' q bin, plus, minus, sum, cumulative Cv'
255:       DUMMY2=0.0D0215:       DUMMY2=0.0D0
256:       DO J2=1,NQBINS216:       DO J2=1,NQBINS
257:          DUMMY2=DUMMY2+CQP(J2)+CQM(J2)217:          DUMMY2=DUMMY2+CQP(J2)+CQM(J2)
258:          IF (CQP(J2)+CQM(J2).GT.1.0D-10) THEN 218:          PRINT '(I8,4G20.10)',J2,CQP(J2),CQM(J2),CQP(J2)+CQM(J2),DUMMY2
259:             PRINT '(I8,4G20.10)',J2,CQP(J2),CQM(J2),CQP(J2)+CQM(J2),DUMMY2 
260:          ENDIF 
261:       ENDDO219:       ENDDO
262:    220:    
263:       PRINT '(A,I8,A)','For disconnectionDPS use line TRMIN 2 ',NMIN,' min.minus.BS min.plus.BS and '221:       STOP
 222:             
 223:    ENDIF
 224: 
 225:    IF (DO2DT) THEN
 226:       PRINT '(A,I8,A)','For disconnectionDPS use line TRMIN 2 ',NMIN,' min.minus.anhshift.BS min.plus.anhshift.BS and '
264:       PRINT '(A)','CHOOSECOLOURS in the dinfo file'227:       PRINT '(A)','CHOOSECOLOURS in the dinfo file'
265:       PRINT '(A)','Assignment of minima to quench bins: minimum, pe, q bin, q energy, diff'228:       PRINT '(A)','Assignment of minima to quench bins: minimum, pe, q bin, q energy, diff'
266:  
267:       NQPLUS=0.0D0 
268:       NQMINUS=0.0D0 
269:       DQPLUS=0.0D0 
270:       DQMINUS=0.0D0 
271:       ALLOCATE(PERMINQP(NQBINS),PERMINQM(NQBINS),MINQP(NQBINS),MINQM(NQBINS)) 
272:  
273:       DO J2=1,NQBINS 
274:          IF (CQP(J2)-CQM(J2).LT.0.0D0) THEN 
275:             NQMINUS=NQMINUS+1 
276:             MINQM(NQMINUS)=J2 
277:             PERMINQM(NQMINUS)=CQM(J2)-CQP(J2) 
278:             DQMINUS=DQMINUS+CQM(J2)-CQP(J2) 
279:          ELSE 
280:             NQPLUS=NQPLUS+1 
281:             MINQP(NQPLUS)=J2 
282:             PERMINQP(NQPLUS)=CQP(J2)-CQM(J2) 
283:             DPLUS=DPLUS+CQP(J2)-CQM(J2) 
284:          ENDIF 
285:       ENDDO 
286:       DO J2=1,NQMINUS 
287:          PERMINQM(J2)=PERMINQM(J2)/MAX(DQMINUS,1.0D-200) 
288:       ENDDO 
289:       DO J2=1,NQPLUS 
290:          PERMINQP(J2)=PERMINQP(J2)/MAX(DQPLUS,1.0D-200) 
291:       ENDDO 
292:       CALL SORT(NQMINUS,NQBINS,PERMINQM,MINQM) 
293:       CALL SORT(NQPLUS,NQBINS,PERMINQP,MINQP) 
294:        
295:       DUMMY=0.0D0 
296:       PRINT '(A)','q bins with negative overall flux below threshold:' 
297:       DO J2=1,NQMINUS 
298:          DUMMY=DUMMY+PERMINQM(J2) 
299:          PRINT '(2I6,4G20.10)',J2,MINQM(J2),PERMINQM(J2),DUMMY,QENERGY(MINQM(J2)) 
300:          IF (DUMMY.GT.COLOURTHRESH) EXIT 
301:       ENDDO 
302:       DUMMY=0.0D0 
303:       PRINT '(A)','q bins with positive overall flux below threshold:' 
304:       DO J2=1,NQPLUS 
305:          DUMMY=DUMMY+PERMINQP(J2) 
306:          PRINT '(2I6,4G20.10)',J2,MINQP(J2),PERMINQP(J2),DUMMY,QENERGY(MINQP(J2)) 
307:          IF (DUMMY.GT.COLOURTHRESH) EXIT 
308:       ENDDO 
309:  
310:       DO J1=1,NMIN229:       DO J1=1,NMIN
311:          DMIN=HUGE(1.0D0)230:          DMIN=HUGE(1.0D0)
312:          DO J2=1,NQBINS231:          DO J2=1,NQBINS
313:             IF (ABS(EMIN(J1)-QENERGY(J2)).LT.DMIN) THEN232:             IF (ABS(EMIN(J1)-QENERGY(J2)).LT.DMIN) THEN
314:                QBIN(J1)=J2233:                QBIN(J1)=J2
315:                DMIN=ABS(EMIN(J1)-QENERGY(J2))234:                DMIN=ABS(EMIN(J1)-QENERGY(J2))
316:             ENDIF235:             ENDIF
317:          ENDDO236:          ENDDO
318: !        PRINT '(I8,G20.10,I8,2G20.10)',J1,EMIN(J1),QBIN(J1),QENERGY(QBIN(J1)),DMIN237:          PRINT '(I8,G20.10,I8,2G20.10)',J1,EMIN(J1),QBIN(J1),QENERGY(QBIN(J1)),DMIN
319:       ENDDO238:       ENDDO
320: 239: 
321:       OPEN(UNIT=1,FILE='min.minus.BS',STATUS='UNKNOWN')240: !     OPEN(UNIT=1,FILE='min.minus.anhshift.BS',STATUS='UNKNOWN')
322:       DUMMY=0.0D0241: !     DUMMY=0.0D0
323:       DO J2=1,NQMINUS242: !     DO J2=1,NMINUS
324:          DUMMY=DUMMY+PERMINQM(J2)243: !        DUMMY=DUMMY+PERMINM(J2)
325:          DO J1=1,NMIN244: !        DO J1=1,NMIN
326:             IF (QBIN(J1).EQ.MINQM(J2)) THEN245: !           IF (QTOI(QBIN(J1)).EQ.MINM(J2)) THEN
327:                WRITE(1,'(I6)') J1246: !              WRITE(1,'(I6)') J1
328:                WRITE(*,'(A,3I8,3G20.10)') 'J1, J2, MINQM(J2), EMIN', J1, J2, MINQM(J2), EMIN(J1)247: !              WRITE(*,'(A,3I8,3G20.10)') 'J1, J2, MINM(J2), EMIN, diff=', &
329:             ENDIF248: ! &                                   J1, J2, MINM(J2), EMIN(J1), &
330:          ENDDO249: ! &                                   ABS(EMIN(J1)-ENERGY(MINM(J2))+ANSHIFT(QBIN(J1)))
331:          IF (DUMMY.GT.COLOURTHRESH) EXIT250: !           ENDIF
332:       ENDDO251: !        ENDDO
333:       CLOSE(1)252: !        IF (DUMMY.GT.COLOURTHRESH) EXIT
 253: !     ENDDO
 254: !     CLOSE(1)
334: 255: 
335:       OPEN(UNIT=1,FILE='min.plus.BS',STATUS='UNKNOWN')256: !     OPEN(UNIT=1,FILE='min.plus.anhshift.BS',STATUS='UNKNOWN')
336:       DUMMY=0.0D0257: !     DUMMY=0.0D0
337:       DO J2=1,NQPLUS258: !     DO J2=1,NPLUS
338:          DUMMY=DUMMY+PERMINQP(J2)259: !        DUMMY=DUMMY+PERMINP(J2)
339:          DO J1=1,NMIN260: !        DO J1=1,NMIN
340:             IF (QBIN(J1).EQ.MINQP(J2)) THEN261: !           IF (QTOI(QBIN(J1)).EQ.MINP(J2)) THEN
341:                WRITE(1,'(I6)') J1262: !              WRITE(1,'(I6)') J1
342:                WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINQP(J2), EMIN=', J1, J2, MINQP(J2), EMIN(J1)263: !              WRITE(*,'(A,3I6,3G20.10)') 'J1, J2, MINP(J2), EMIN, diff=', &
343:             ENDIF264: ! &                                   J1, J2, MINP(J2), EMIN(J1), &
344:          ENDDO265: ! &                                   ABS(EMIN(J1)-ENERGY(MINP(J2))+ANSHIFT(QBIN(J1)))
345:          IF (DUMMY.GT.COLOURTHRESH) EXIT266: !           ENDIF
346:       ENDDO 267: !        ENDDO
347:       CLOSE(1)268: !        IF (DUMMY.GT.COLOURTHRESH) EXIT
 269: !     ENDDO 
 270: !     CLOSE(1)
348: 271: 
349:    ENDIF272:    ENDIF
350: ENDIF273: ENDIF
351: 274: 
352: END PROGRAM CVBS275: END PROGRAM CVBS
353: 276: 
354: SUBROUTINE SORT(N,J3,A,NA)277: SUBROUTINE SORT(N,J3,A,NA)
355: IMPLICIT NONE278: IMPLICIT NONE
356: INTEGER J1, L, N, J3, J2, NA(J3), NTEMP279: INTEGER J1, L, N, J3, J2, NA(J3), NTEMP
357: DOUBLE PRECISION TEMP, A(J3)280: DOUBLE PRECISION TEMP, A(J3)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0