I would like to ask about the error occurred with my Fortran code. Since I'm new to fortran I can't handle this after 2 days, I also searched around but still don't know how to fix it.
PROGRAM SUBDEM
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 NTET,NPTK
INTEGER*4 NA,NC,NE,NBAND,NTMAX,IDEF
C CALL SETK08(NA,NC, A,C, PTK,NPTK, IDEF,NTET, NKMAX,NTMAX)
C CALL DOSTET(ENER,IDIME,NBAND,IDEF,NTET,XE,NE,Y,Z)
PRINT *,'Test print'
C with DIMENSION PTK(4,NKMAX),IDEF(5,NTMAX)
C NA,NC are the two k-mesh discretization parameters,
C A,C are the two parameters "a","c" of the direct lattice.
C PTK is a REAL*8 array and IDEF is an INTEGER*4 array
NA=3
NC=3
A=1.67
C=2.1
NKMAX=550
NTMAX=50.D3
CALL SETK08(NA,NC,A,C,PTK,NPTK,IDEF,NTET,NKMAX,NTMAX)
NE = 20
IDIME = 5.0
NBAND = 5
XE = 4
NE = 15
C CALL DOSTET(ENER,IDIME,NBAND,IDEF,NTET,XE,NE,Y,Z)
C ENER (REAL*8 two-dimensional array, input)
C ENER(NU,IK) is the energy of band NU, computed for the k-point IK, defined by the SETK** routine
C IDIME (INTEGER*4, input)
C First dimension of the array ENER, as defined in the calling program. IDIME must be at least equal to NBAND
C NBAND (INTEGER*4, input)
C Number of energy bands included in the summation
C IDEF (INTEGER*4 two-dimensional array,input)
C Table defining the tetrahedron corners, as obtained from the SETK** routines. The first dimension is 5.
C NTET (INTEGER*4, input)
C Number of tetrahedra filling the volume V (provided by a SETK** routine)
C XE (REAL*8 one-dimensional array, input)
C Contains the values of the energies E where the density of states and integrated density of states are to be computed.
C Dimension is at least NE.
C NE (INTEGER*4, input)
C Number of energy points where the density of states and integrated density of states are computed. Only the first NE locations of XE are used by DOSTET.
C Y (REAL*8 one-dimensional array, output)
C The NE first components of this vector contain, on return, the density of states evaluated at energy points corresponding to the NE first components of XE.
C Z (REAL*8 one-dimensional array, output)
C The NE first components of this vector contain, on return, the integrated density of states evaluated at energy points corresponding to the NE first components of XE.
END
SUBROUTINE SETK08(NA,NC,A,C,PTK,NPTK,IDEF,NTET,NKMAX,NTMAX)
C SET THE K-POINTS IN THE 1/16TH OF THE BRILLOUIN ZONE FOR A
C SIMPLE TETRAGONAL LATTICE WITH PARAMETERS A, C
C SYMMETRY IS D4H
IMPLICIT REAL*8(A-H,O-Z)
REAL*4 AVOL
DIMENSION PTK(4,NKMAX),IDEF(5,NTMAX)
EQUIVALENCE (IVOL,AVOL)
PI = 3.141592653589793238D0
IF(NA.LE.0.OR.NC.LE.0) GOTO 97
IF(A.LE.0.0D0 .OR. C.LE.0.0D0) GOTO 98
NPTK = (NA+1)*(NA+2)*(NC+1)/2
IF(NPTK.GT.NKMAX) STOP '*** <SETK08> NPTK EXCEEDS NKMAX ***'
NTET = 3*NC*NA**2
IF(NTET.GT.NTMAX) STOP '*** <SETK08> NTET EXCEEDS NTMAX ***'
C *** SET THE K-POINTS
AK=PI/A/NA
CK=PI/C/NC
WRITE(6,100) NPTK,NTET,NA*AK,NA*AK,NC*CK
W = 2.0D0/(NA*NA*NC)
NPTK=0
DO 1 I=0,NA,1
DO 1 J=0,I,1
DO 1 K=0,NC,1
C NPTK = I*(I+1)/2*NZ1 + J*NZ1 + K+1
WK = W
IF(I.EQ.0) WK = WK/2.0D0
IF(J.EQ.0) WK = WK/2.0D0
IF(J.EQ.I) WK = WK/2.0D0
IF(I.EQ.NA) WK = WK/2.0D0
IF(J.EQ.NA) WK = WK/2.0D0
IF(K.EQ.0 .OR. K.EQ.NC) WK = WK/2.0D0
NPTK=NPTK+1
PTK(1,NPTK)=I*AK
PTK(2,NPTK)=J*AK
PTK(3,NPTK)=K*CK
PTK(4,NPTK)=WK
1 CONTINUE
C *** DEFINE THE TETRAHEDRA
NZ1=NC+1
NTET=0
I7=0
I=0
4 IX=(I+1)*NZ1
J = 0
5 K=0
I7=I*IX/2+J*NZ1
6 I7=I7+1
I6=I7+IX
I2=I6+NZ1
I1=I2+1
NTET=NTET+1
IDEF(1,NTET)=I7
IDEF(2,NTET)=I6
IDEF(3,NTET)=I2
IDEF(4,NTET)=I1
I8=I7+1
I5=I6+1
NTET=NTET+1
IDEF(1,NTET)=I7
IDEF(2,NTET)=I6
IDEF(3,NTET)=I5
IDEF(4,NTET)=I1
NTET=NTET+1
IDEF(1,NTET)=I7
IDEF(2,NTET)=I8
IDEF(3,NTET)=I5
IDEF(4,NTET)=I1
IF(J.EQ.I) GOTO 7
I3=I7+NZ1
I4=I3+1
NTET=NTET+1
IDEF(1,NTET)=I7
IDEF(2,NTET)=I3
IDEF(3,NTET)=I2
IDEF(4,NTET)=I1
NTET=NTET+1
IDEF(1,NTET)=I7
IDEF(2,NTET)=I3
IDEF(3,NTET)=I4
IDEF(4,NTET)=I1
NTET=NTET+1
IDEF(1,NTET)=I7
IDEF(2,NTET)=I8
IDEF(3,NTET)=I4
IDEF(4,NTET)=I1
7 K=K+1
IF(K.LT.NC) GOTO 6
J=J+1
IF(J.LE.I) GOTO 5
I=I+1
IF(I.LT.NA) GOTO 4
AVOL=1.D0/DFLOAT(NTET)
DO 15 IT=1,NTET
15 IDEF(5,IT)=IVOL
PRINT *,NTET,NPTK
RETURN
97 WRITE(6,101)
GOTO 99
98 WRITE(6,102)
99 STOP
100 FORMAT(' SAMPLING THE 16TH PART OF A SQUARE-BASED PRISM'/
.1X,I5,' K-POINTS',I7,' TETRAHEDRA'/
.' KXMAX =',D11.4,' KYMAX =',D11.4,' KZMAX =',D11.4)
101 FORMAT(' *** <SETK08> NA OR NC IS NOT A POSITIVE INTEGER ***')
102 FORMAT(' *** <SETK08> A AND C MUST BE POSITIVE ***')
END
SUBROUTINE DOSTET(ENER,IDIME,NBAND,IDEF,NTET,XE,NE,Y,Z)
C COMPUTE A DENSITY OF STATES USING THE TETRAHEDRONS METHOD.
C XE CONTAINS THE ENERGIES, Y AND Z RETURN THE RELATED DENSITY OF
C STATES AND THE INTEGRATED DENSITY OF STATES, RESPECTIVELY.
IMPLICIT REAL*8(A-H,O-Z)
REAL*4 AVOL
DIMENSION ENER(IDIME,1),XE(1),Y(1),Z(1),IDEF(5,1),C(4)
EQUIVALENCE (IVOL,AVOL),(C(1),E1),(C(2),E2),(C(3),E3),(C(4),E4)
DATA EPS/1.0D-05/
DO 6 IX=1,NE
Y(IX)=0.D0
6 Z(IX)=0.D0
C
C LOOP OVER THE TETRAHEDRONS
DO 9 ITET=1,NTET
C
IA=IDEF(1,ITET)
IB=IDEF(2,ITET)
IC=IDEF(3,ITET)
ID=IDEF(4,ITET)
IVOL=IDEF(5,ITET)
C
C LOOP OVER THE BANDS
DO 9 NB=1,NBAND
C
C *** DEFINE E1, E2, E3, E4, AS THE CORNER ENERGIES ORDERED BY
C *** DECREASING SIZE
C(1)=ENER(NB,IA)
C(2)=ENER(NB,IB)
C(3)=ENER(NB,IC)
C(4)=ENER(NB,ID)
DO 2 I=1,4
CC=C(I)
J=I
1 J=J+1
IF(J.GT.4) GOTO 2
IF(CC.GE.C(J)) GOTO 1
C(I)=C(J)
C(J)=CC
CC=C(I)
GOTO 1
2 CONTINUE
UNITE=1.0D0
IF(E1.GT.E4) UNITE=E1-E4
E12=(E1-E2)/UNITE
E13=(E1-E3)/UNITE
E14=(E1-E4)/UNITE
E23=(E2-E3)/UNITE
E24=(E2-E4)/UNITE
E34=(E3-E4)/UNITE
FACY=3.D0*DBLE(AVOL)/UNITE
DO 9 IX=1,NE
E=XE(IX)
SURFAC=0.D0
VOLUME=1.D0
IF(E.GT.E1) GOTO 8
VOLUME=0.D0
IF(E.LT.E4) GOTO 8
EE1=(E-E1)/UNITE
IF(DABS(EE1).LT.EPS) EE1=0.D0
EE2=(E-E2)/UNITE
IF(DABS(EE2).LT.EPS) EE2=0.D0
EE3=(E-E3)/UNITE
IF(DABS(EE3).LT.EPS) EE3=0.D0
EE4=(E-E4)/UNITE
IF(DABS(EE4).LT.EPS) EE4=0.D0
IF(E.GT.E3) GOTO 5
C *** E4.LE.E.AND.E.LE.E3
IF(E4.EQ.E3) GOTO 3
SURFAC=(EE4/E34)*(EE4/E24)
VOLUME=SURFAC*EE4
GOTO 8
3 IF(E3.LT.E2) GOTO 8
IF(E2.EQ.E1) GOTO 4
SURFAC=1.D0/E12
GOTO 8
4 SURFAC=1.0D+15
VOLUME=0.5D0
GOTO 8
5 IF(E.GT.E2) GOTO 7
C *** E3.LT.E.AND.E.LE.E2
SURFAC=-(EE3*EE2/E23+EE4*EE1)/E13/E24
VOLUME=(0.5D0*EE3*(2.D0*E13*E34+E13*EE4-E34*EE1-2.D0*EE1*EE4+
+ EE3*(EE3-3.D0*EE2)/E23)/E13+E34*E34)/E24
GOTO 8
C *** E2.LT.E.AND.E.LE.E1
7 SURFAC=(EE1/E12)*(EE1/E13)
VOLUME=1.D0+SURFAC*EE1
8 Y(IX)=Y(IX)+FACY*SURFAC
Z(IX)=Z(IX)+DBLE(AVOL)*VOLUME
9 CONTINUE
RETURN
END
Seems like the code broke at line 82:
PTK(1,NPTK)=DFLOAT(I)*DK
