c This program does the sum over l for the modes for given values of omega IMPLICIT NONE Integer i,ip,j,k,kk,li,lf,npts,is,n,l,LL,ltest,nom,junk,nab real*16 s,sin,dels,tlp1,xl,l12,om,omtest,stest,omi,delom real*16 S1,S2,S3,S4,S5,S2l,S3l,S4l,S5l real*16 difS2,difS3,difS4,difS5,pq,pdqqdp,dpdq real*16 Pi real*16 xi,a2,a3,a4,b2,b3,b4,A22,A33,B22,B33,A44,B44 real*16 S2acc,S3acc,S4acc,S5acc,S2Lacc,S3Lacc,S4Lacc,S5Lacc real*16 S2cutacc,S3cutacc,S4cutacc,S5cutacc real*16 pqacc,pdqqdpacc,dpdqacc Dimension lf(50),nom(50),omi(50),delom(25),S2l(200),S3l(200),S4l(200),S5l(200) Dimension pqacc(200),pdqqdpacc(200),dpdqacc(200) OPEN(12,FILE='s.dat',STATUS='OLD') OPEN(14,FILE='lz.dat',STATUS='OLD') Open(15,file="Slsumz.p1") Open(16,file="Slsumz.p2") Open(17,file="Slsumz.p3") Open(18,file='Slsum.tst') Open(19,file="Slsumz.out") Open(10,file='Slsum_acc_all.out') Open(11,file='Slsum_acc.out') c Open(9,file="Slsumz.test") Pi = 2.q0*qasin(1.q0) Print *,Pi c sin is the value of s you want to start computing the stress tensor at Read(12,*) Junk Read(12,*) sin,dels,npts READ(14,*) j Read(14,*) omi(1),delom(1),nom(1),lf(1) Print *,' omi(1) = ',omi(1),' delom(1) = ',delom(1), 'nom(1) = ',nom(1),' lf(1) = ',lf(1) Print *,' omi(1) = ',omi(1),' delom(1) = ',delom(1), ' nom(1) = ',nom(1),' lf(1) = ',lf(1) c Due to lack of disk space I have to put in the actual pqew8 files. OPEN(21,FILE='pqw8z.1',STATUS='OLD') Read(21,*) junk READ(21,*) xi,a2,b2,a3,b3,a4,b4 Read(21,*) junk READ(21,*) nab,A22,B22,A33,B33,A44,B44 Read(21,*) junk WRITE(6,*) xi,a2,b2,a3,b3,a4,b4 Write(6,*) nab,A22,B22,A33,B33,A44,B44 OPEN(41,FILE='pqacc.1',STATUS='OLD') Read(41,*) junk READ(41,*) xi,a2,b2,a3,b3 Read(41,*) junk READ(41,*) nab,A22,B22,A33,B33 Read(41,*) junk Read(41,*) junk If(j.gt.1) then Read(14,*) omi(2),delom(2),nom(2),lf(2) Print *,' omi(2) = ',omi(2),' delom(2) = ',delom(2), 'nom(2) = ',nom(2),' lf(2) = ',lf(2) OPEN(22,FILE='pqw8z.2',STATUS='OLD') Read(22,*) junk READ(22,*) xi,a2,b2,a3,b3 Read(22,*) junk READ(22,*) nab,A22,B22,A33,B33 Read(22,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(42,FILE='pqacc.2',STATUS='OLD') Read(42,*) junk READ(42,*) xi,a2,b2,a3,b3 Read(42,*) junk READ(42,*) nab,A22,B22,A33,B33 Read(42,*) junk Read(42,*) junk Endif If(j.gt.2) then Read(14,*) omi(3),delom(3),nom(3),lf(3) Print *,' omi(3) = ',omi(3),' delom(3) = ',delom(3), 'nom(3) = ',nom(3),' lf(3) = ',lf(3) OPEN(23,FILE='pqw8z.3',STATUS='OLD') Read(23,*) junk READ(23,*) xi,a2,b2,a3,b3 Read(23,*) junk READ(23,*) nab,A22,B22,A33,B33 Read(23,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(43,FILE='pqacc.3',STATUS='OLD') Read(43,*) junk READ(43,*) xi,a2,b2,a3,b3 Read(43,*) junk READ(43,*) nab,A22,B22,A33,B33 Read(43,*) junk Read(43,*) junk Endif If(j.gt.3) then Read(14,*) omi(4),delom(4),nom(4),lf(4) Print *,' omi(4) = ',omi(4),' delom(4) = ',delom(4), 'nom(4) = ',nom(4),' lf(4) = ',lf(4) OPEN(24,FILE='pqw8z.4',STATUS='OLD') Read(24,*) junk READ(24,*) xi,a2,b2,a3,b3 Read(24,*) junk READ(24,*) nab,A22,B22,A33,B33 Read(24,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(44,FILE='pqacc.4',STATUS='OLD') Read(44,*) junk READ(44,*) xi,a2,b2,a3,b3 Read(44,*) junk READ(44,*) nab,A22,B22,A33,B33 Read(44,*) junk Read(44,*) junk Endif If(j.gt.4) then Read(14,*) omi(5),delom(5),nom(5),lf(5) Print *,' omi(5) = ',omi(5),' delom(5) = ',delom(5), 'nom(5) = ',nom(5),' lf(5) = ',lf(5) OPEN(25,FILE='pqw8z.5',STATUS='OLD') Read(25,*) junk READ(25,*) xi,a2,b2,a3,b3 Read(25,*) junk READ(25,*) nab,A22,B22,A33,B33 Read(25,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(45,FILE='pqacc.5',STATUS='OLD') Read(45,*) junk READ(45,*) xi,a2,b2,a3,b3 Read(45,*) junk READ(45,*) nab,A22,B22,A33,B33 Read(45,*) junk Read(45,*) junk Endif If(j.gt.5) then Read(14,*) omi(6),delom(6),nom(6),lf(6) Print *,' omi(6) = ',omi(6),' delom(6) = ',delom(6), 'nom(6) = ',nom(6),' lf(6) = ',lf(6) OPEN(26,FILE='pqw8z.6',STATUS='OLD') Read(26,*) junk READ(26,*) xi,a2,b2,a3,b3 Read(26,*) junk READ(26,*) nab,A22,B22,A33,B33 Read(26,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(46,FILE='pqacc.6',STATUS='OLD') Read(46,*) junk READ(46,*) xi,a2,b2,a3,b3 Read(46,*) junk READ(46,*) nab,A22,B22,A33,B33 Read(46,*) junk Read(46,*) junk Endif If(j.gt.6) then Read(14,*) omi(7),delom(7),nom(7),lf(7) Print *,' omi(7) = ',omi(7),' delom(7) = ',delom(7), 'nom(7) = ',nom(7),' lf(7) = ',lf(7) OPEN(27,FILE='pqw8z.7',STATUS='OLD') Read(27,*) junk READ(27,*) xi,a2,b2,a3,b3 Read(27,*) junk READ(27,*) nab,A22,B22,A33,B33 Read(27,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(47,FILE='pqacc.7',STATUS='OLD') Read(47,*) junk READ(47,*) xi,a2,b2,a3,b3 Read(47,*) junk READ(47,*) nab,A22,B22,A33,B33 Read(47,*) junk Read(47,*) junk Endif If(j.gt.7) then Read(14,*) omi(8),delom(8),nom(8),lf(8) Print *,' omi(8) = ',omi(8),' delom(8) = ',delom(8), 'nom(8) = ',nom(8),' lf(8) = ',lf(8) OPEN(28,FILE='pqw8z.8',STATUS='OLD') Read(28,*) junk READ(28,*) xi,a2,b2,a3,b3 Read(28,*) junk READ(28,*) nab,A22,B22,A33,B33 Read(28,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(48,FILE='pqacc.8',STATUS='OLD') Read(48,*) junk READ(48,*) xi,a2,b2,a3,b3 Read(48,*) junk READ(48,*) nab,A22,B22,A33,B33 Read(48,*) junk Read(48,*) junk Endif If(j.gt.8) then Read(14,*) omi(9),delom(9),nom(9),lf(9) Print *,' omi(9) = ',omi(9),' delom(9) = ',delom(9), 'nom(9) = ',nom(9),' lf(9) = ',lf(9) OPEN(29,FILE='pqw8z.9',STATUS='OLD') Read(29,*) junk READ(29,*) xi,a2,b2,a3,b3 Read(29,*) junk READ(29,*) nab,A22,B22,A33,B33 Read(29,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(49,FILE='pqacc.9',STATUS='OLD') Read(49,*) junk READ(49,*) xi,a2,b2,a3,b3 Read(49,*) junk READ(49,*) nab,A22,B22,A33,B33 Read(49,*) junk Read(49,*) junk Endif If(j.gt.9) then Read(14,*) omi(10),delom(10),nom(10),lf(10) Print *,' omi(10) = ',omi(10),' delom(10) = ',delom(10), 'nom(10) = ',nom(10),' lf(10) = ',lf(10) OPEN(30,FILE='pqw8z.10',STATUS='OLD') Read(30,*) junk READ(30,*) xi,a2,b2,a3,b3 Read(30,*) junk READ(30,*) nab,A22,B22,A33,B33 Read(30,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(50,FILE='pqacc.10',STATUS='OLD') Read(50,*) junk READ(50,*) xi,a2,b2,a3,b3 Read(50,*) junk READ(50,*) nab,A22,B22,A33,B33 Read(50,*) junk Read(50,*) junk Endif If(j.gt.10) then Read(14,*) omi(11),delom(11),nom(11),lf(11) Print *,' omi(11) = ',omi(11),' delom(11) = ',delom(11), 'nom(11) = ',nom(11),' lf(11) = ',lf(11) OPEN(31,FILE='pqw8z.11',STATUS='OLD') Read(31,*) junk READ(31,*) xi,a2,b2,a3,b3 Read(31,*) junk READ(31,*) nab,A22,B22,A33,B33 Read(31,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(51,FILE='pqacc.11',STATUS='OLD') Read(51,*) junk READ(51,*) xi,a2,b2,a3,b3 Read(51,*) junk READ(51,*) nab,A22,B22,A33,B33 Read(51,*) junk Read(51,*) junk Endif If(j.gt.11) then Read(14,*) omi(12),delom(12),nom(12),lf(12) Print *,' omi(12) = ',omi(12),' delom(12) = ',delom(12), 'nom(12) = ',nom(12),' lf(12) = ',lf(12) OPEN(32,FILE='pqw8z.12',STATUS='OLD') Read(32,*) junk READ(32,*) xi,a2,b2,a3,b3 Read(32,*) junk READ(32,*) nab,A22,B22,A33,B33 Read(32,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(52,FILE='pqacc.12',STATUS='OLD') Read(52,*) junk READ(52,*) xi,a2,b2,a3,b3 Read(52,*) junk READ(52,*) nab,A22,B22,A33,B33 Read(52,*) junk Read(52,*) junk Endif If(j.gt.12) then Read(14,*) omi(13),delom(13),nom(13),lf(13) Print *,' omi(13) = ',omi(13),' delom(13) = ',delom(13), 'nom(13) = ',nom(13),' lf(13) = ',lf(13) OPEN(33,FILE='pqw8z.13',STATUS='OLD') Read(33,*) junk READ(33,*) xi,a2,b2,a3,b3 Read(33,*) junk READ(33,*) nab,A22,B22,A33,B33 Read(33,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(53,FILE='pqacc.13',STATUS='OLD') Read(53,*) junk READ(53,*) xi,a2,b2,a3,b3 Read(53,*) junk READ(53,*) nab,A22,B22,A33,B33 Read(53,*) junk Read(53,*) junk Endif If(j.gt.13) then Read(14,*) omi(14),delom(14),nom(14),lf(14) Print *,' omi(14) = ',omi(14),' delom(14) = ',delom(14), 'nom(14) = ',nom(14),' lf(14) = ',lf(14) OPEN(34,FILE='pqw8z.14',STATUS='OLD') Read(34,*) junk READ(34,*) xi,a2,b2,a3,b3 Read(34,*) junk READ(34,*) nab,A22,B22,A33,B33 Read(34,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(54,FILE='pqacc.14',STATUS='OLD') Read(54,*) junk READ(54,*) xi,a2,b2,a3,b3 Read(54,*) junk READ(54,*) nab,A22,B22,A33,B33 Read(54,*) junk Read(54,*) junk Endif If(j.gt.14) then Read(14,*) omi(15),delom(15),nom(15),lf(15) Print *,' omi(15) = ',omi(15),' delom(15) = ',delom(15), 'nom(15) = ',nom(15),' lf(15) = ',lf(15) OPEN(35,FILE='pqw8z.15',STATUS='OLD') Read(35,*) junk READ(35,*) xi,a2,b2,a3,b3 Read(35,*) junk READ(35,*) nab,A22,B22,A33,B33 Read(35,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(55,FILE='pqacc.15',STATUS='OLD') Read(55,*) junk READ(55,*) xi,a2,b2,a3,b3 Read(55,*) junk READ(55,*) nab,A22,B22,A33,B33 Read(55,*) junk Read(55,*) junk Endif If(j.gt.15) then Read(14,*) omi(16),delom(16),nom(16),lf(16) Print *,' omi(16) = ',omi(16),' delom(16) = ',delom(16), 'nom(16) = ',nom(16),' lf(16) = ',lf(16) OPEN(36,FILE='pqw8z.16',STATUS='OLD') Read(36,*) junk READ(36,*) xi,a2,b2,a3,b3 Read(36,*) junk READ(36,*) nab,A22,B22,A33,B33 Read(36,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(56,FILE='pqacc.16',STATUS='OLD') Read(56,*) junk READ(56,*) xi,a2,b2,a3,b3 Read(56,*) junk READ(56,*) nab,A22,B22,A33,B33 Read(56,*) junk Read(56,*) junk Endif If(j.gt.16) then Read(14,*) omi(17),delom(17),nom(17),lf(17) Print *,' omi(17) = ',omi(17),' delom(17) = ',delom(17), 'nom(17) = ',nom(17),' lf(17) = ',lf(17) OPEN(37,FILE='pqw8z.17',STATUS='OLD') Read(37,*) junk READ(37,*) xi,a2,b2,a3,b3 Read(37,*) junk READ(37,*) nab,A22,B22,A33,B33 Read(37,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(57,FILE='pqacc.17',STATUS='OLD') Read(57,*) junk READ(57,*) xi,a2,b2,a3,b3 Read(57,*) junk READ(57,*) nab,A22,B22,A33,B33 Read(57,*) junk Read(57,*) junk Endif If(j.gt.17) then Read(14,*) omi(18),delom(18),nom(18),lf(18) Print *,' omi(18) = ',omi(18),' delom(18) = ',delom(18), 'nom(18) = ',nom(18),' lf(18) = ',lf(18) OPEN(38,FILE='pqw8z.18',STATUS='OLD') Read(38,*) junk READ(38,*) xi,a2,b2,a3,b3 Read(38,*) junk READ(38,*) nab,A22,B22,A33,B33 Read(38,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(58,FILE='pqacc.18',STATUS='OLD') Read(58,*) junk READ(58,*) xi,a2,b2,a3,b3 Read(58,*) junk READ(58,*) nab,A22,B22,A33,B33 Read(58,*) junk Read(58,*) junk Endif If(j.gt.18) then Read(14,*) omi(19),delom(19),nom(19),lf(19) Print *,' omi(19) = ',omi(19),' delom(19) = ',delom(19), 'nom(19) = ',nom(19),' lf(19) = ',lf(19) OPEN(39,FILE='pqw8z.19',STATUS='OLD') Read(39,*) junk READ(39,*) xi,a2,b2,a3,b3 Read(39,*) junk READ(39,*) nab,A22,B22,A33,B33 Read(39,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(59,FILE='pqacc.19',STATUS='OLD') Read(59,*) junk READ(59,*) xi,a2,b2,a3,b3 Read(59,*) junk READ(59,*) nab,A22,B22,A33,B33 Read(59,*) junk Read(59,*) junk Endif If(j.gt.19) then Read(14,*) omi(20),delom(20),nom(20),lf(20) Print *,' omi(20) = ',omi(20),' delom(20) = ',delom(20), 'nom(20) = ',nom(20),' lf(20) = ',lf(20) OPEN(40,FILE='pqw8z.20',STATUS='OLD') Read(40,*) junk READ(40,*) xi,a2,b2,a3,b3 Read(40,*) junk READ(40,*) nab,A22,B22,A33,B33 Read(40,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(60,FILE='pqacc.20',STATUS='OLD') Read(60,*) junk READ(60,*) xi,a2,b2,a3,b3 Read(60,*) junk READ(60,*) nab,A22,B22,A33,B33 Read(60,*) junk Read(60,*) junk Endif If(j.gt.20) then Read(14,*) omi(21),delom(21),nom(21),lf(21) Print *,' omi(21) = ',omi(21),' delom(21) = ',delom(21), 'nom(21) = ',nom(21),' lf(21) = ',lf(21) OPEN(61,FILE='pqw8z.21',STATUS='OLD') Read(61,*) junk READ(61,*) xi,a2,b2,a3,b3 Read(61,*) junk READ(61,*) nab,A22,B22,A33,B33 Read(61,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(81,FILE='pqacc.21',STATUS='OLD') Read(81,*) junk READ(81,*) xi,a2,b2,a3,b3 Read(81,*) junk READ(81,*) nab,A22,B22,A33,B33 Read(81,*) junk Read(81,*) junk Endif If(j.gt.21) then Read(14,*) omi(22),delom(22),nom(22),lf(22) Print *,' omi(22) = ',omi(22),' delom(22) = ',delom(22), 'nom(22) = ',nom(22),' lf(22) = ',lf(22) OPEN(62,FILE='pqw8z.22',STATUS='OLD') Read(62,*) junk READ(62,*) xi,a2,b2,a3,b3 Read(62,*) junk READ(62,*) nab,A22,B22,A33,B33 Read(62,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(82,FILE='pqacc.22',STATUS='OLD') Read(82,*) junk READ(82,*) xi,a2,b2,a3,b3 Read(82,*) junk READ(82,*) nab,A22,B22,A33,B33 Read(82,*) junk Read(82,*) junk Endif If(j.gt.22) then Read(14,*) omi(23),delom(23),nom(23),lf(23) Print *,' omi(23) = ',omi(23),' delom(23) = ',delom(23), 'nom(23) = ',nom(23),' lf(23) = ',lf(23) OPEN(63,FILE='pqw8z.23',STATUS='OLD') Read(63,*) junk READ(63,*) xi,a2,b2,a3,b3 Read(63,*) junk READ(63,*) nab,A22,B22,A33,B33 Read(63,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(83,FILE='pqacc.23',STATUS='OLD') Read(83,*) junk READ(83,*) xi,a2,b2,a3,b3 Read(83,*) junk READ(83,*) nab,A22,B22,A33,B33 Read(83,*) junk Read(83,*) junk Endif If(j.gt.23) then Read(14,*) omi(24),delom(24),nom(24),lf(24) Print *,' omi(24) = ',omi(24),' delom(24) = ',delom(24), 'nom(24) = ',nom(24),' lf(24) = ',lf(24) OPEN(64,FILE='pqw8z.24',STATUS='OLD') Read(64,*) junk READ(64,*) xi,a2,b2,a3,b3 Read(64,*) junk READ(64,*) nab,A22,B22,A33,B33 Read(64,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(84,FILE='pqacc.24',STATUS='OLD') Read(84,*) junk READ(84,*) xi,a2,b2,a3,b3 Read(84,*) junk READ(84,*) nab,A22,B22,A33,B33 Read(84,*) junk Read(84,*) junk Endif If(j.gt.24) then Read(14,*) omi(25),delom(25),nom(25),lf(25) Print *,' omi(25) = ',omi(25),' delom(25) = ',delom(25), 'nom(25) = ',nom(25),' lf(25) = ',lf(25) OPEN(65,FILE='pqw8z.25',STATUS='OLD') Read(65,*) junk READ(65,*) xi,a2,b2,a3,b3 Read(65,*) junk READ(65,*) nab,A22,B22,A33,B33 Read(65,*) junk WRITE(6,*) xi,a2,b2,a3,b3 Write(6,*) nab,A22,B22,A33,B33 OPEN(85,FILE='pqacc.25',STATUS='OLD') Read(85,*) junk READ(85,*) xi,a2,b2,a3,b3 Read(85,*) junk READ(85,*) nab,A22,B22,A33,B33 Read(85,*) junk Read(85,*) junk Endif c c We can add to the above if necessary but need to redo the numbers above. S = sin - dELS Do is=1,NPTS S = S + dels Print *,s c Now compute a running sum of and and compute the convergence c of the sum. c c This is the do loop that runs over the various data files. c Do i=1,j If(i.le.20) then k = 20+i kk = 40+i Else k = 40+i kk = 60+i Endif Print *,'i = ',i,' k = ',k,' kk = ',kk om = omi(i)-delom(i) Do n=1,nom(i) om = om+delom(i) S2 = 0.q0 S3 = 0.q0 S4 = 0.q0 S5 = 0.q0 Do L=0,LF(i) LL = L+1 tlp1 = 2.q0*l+1.q0 xl = l l12 = xl+0.5q0 c I am going to assume the pqacc files are the same as the pq files. 19 READ(kk,*,end=999) STEST,omtest,LTEST,pqacc(LL),pdqqdpacc(LL),dpdqacc(LL) READ(k,*,end=9999) STEST,omtest,LTEST,pq,pdqqdp,dpdq IF(qabs(STEST-S).LT.1.q-14.and.qabs(om-omtest).lt.1.q-14.and.LTEST.EQ.L) THEN Continue ELSE Print *,'s = ',s,' stest = ',stest,' om = ',om,' omtest = ',omtest Print *,'L = ',L,' Ltest = ',Ltest GO TO 19 ENDIF S2l(LL) = tlp1*dpdq S3l(LL) = tlp1*l12**2*pq S4l(LL) = tlp1*pdqqdp S5l(LL) = tlp1*pq S2 = S2+S2l(LL) S3 = S3+S3l(LL) S4 = S4+S4l(LL) S5 = S5+S5l(LL) If(is.eq.1.or.is.eq.2.or.is.eq.10.or.is.eq.20.or.is.eq.30.or.is.eq.npts) then difS2 = S2l(LL)/S2 difS3 = S3l(LL)/S3 difS4 = S4l(LL)/S4 difS5 = S5l(LL)/S5 write(18,37) k,s,om,l,S2,difS2,S3,difS3,S4,difS4,S5,difS5 Endif Enddo c Compute the accuracy of the sums over L Do L=0,LF(i) LL = L+1 c If S2l is zero they are all zero. If(qabs(S2l(LL)).gt.1.q-15) then c Here we take the minimum accuracy for each nonzero contribution to each mode sum. c It is computed by taking the accuracy of the relevant product of modes and then subtracting log(difSn) c The latter takes into account the fact that some sums may not contribute much and so their accuracy c is enhanced. It also takes into account roundoff by having less digits if the contribution is greater than the total sum. difS2 = S2l(LL)/S2 difS3 = S3l(LL)/S3 difS4 = S4l(LL)/S4 difS5 = S5l(LL)/S5 S2Lacc = dpdqacc(LL) - qlog10(qabs(difS2)) S3Lacc = pqacc(LL) - qlog10(qabs(difS3)) S4Lacc = pdqqdpacc(LL) - qlog10(qabs(difS4)) S5Lacc = pqacc(LL) - qlog10(qabs(difS5)) If(L.eq.0) then S2acc = S2Lacc S3acc = S3Lacc S4acc = S4Lacc S5acc = S5Lacc Else S2acc = qmin1(S2acc,S2Lacc) S3acc = qmin1(S3acc,S3Lacc) S4acc = qmin1(S4acc,S4Lacc) S5acc = qmin1(S5acc,S5Lacc) Endif Write(10,18) s,om,L,S2Lacc,S3Lacc,S4Lacc,S5Lacc c Write(9,18) s,om,L,pqacc(LL),pdqqdpacc(LL),dpdqacc(LL) Endif Enddo c The cutoff accuracy is the one digit less than the contribution of difSn for the last term that does not have pqacc = 0. S2cutacc = -qlog10(qabs(difS2)) - 1.q0 S3cutacc = -qlog10(qabs(difS3)) - 1.q0 S4cutacc = -qlog10(qabs(difS4)) - 1.q0 S5cutacc = -qlog10(qabs(difS5)) - 1.q0 Write(10,*) '0 Accuracy of sums over L from roundoff of WKB subtraction and contribution to the sum' Write(10,29) s,om,S2acc,S3acc,S4acc,S5acc Write(10,*) '0 Cutoff accuracy of sums over L' Write(10,29) s,om,S2cutacc,S3cutacc,S4cutacc,S5cutacc S2acc = qmin1(S2acc,S2cutacc) S3acc = qmin1(S3acc,S3cutacc) S4acc = qmin1(S4acc,S4cutacc) S5acc = qmin1(S5acc,S5cutacc) Write(10,*) '0 Total accuracy of sums over L' Write(10,29) s,om,S2acc,S3acc,S4acc,S5acc Write(11,29) s,om,S2acc,S3acc,S4acc,S5acc c Now write the actual sums S1 = om**2*S5 Write(19,17) s,om,S1,S2,S3,S4,S5 If(is.eq.1) then Write(15,17) s,om,S1,S2,S3,S4,S5 else if(is.eq.5) then Write(16,17) s,om,S1,S2,S3,S4,S5 else if(is.eq.npts) then Write(17,17) s,om,S1,S2,S3,S4,S5 Endif Enddo Enddo Enddo 17 format(1x,1(f10.7,1x),6(e22.15,1x)) 18 format(1x,2(d12.5,1x),I3,1x,4(f5.1,1x)) 29 format(1x,2(d16.9,1x),4(f5.1,1x)) 37 format(1x,i3,1x,2(e12.5,1x),1x,i4,2x,9(e12.5,1x)) STOP 999 Print *,'Read to the end of pq_acc. for i = ',i,' k = ',k Print *,'s = ',s,' stest = ',stest,' om = ',om,' omtest = ',omtest Print *,' l = ',l,' ltest = ',ltest Stop 9999 Print *,'Read to the end of pqw8z. for i = ',i,' k = ',k Print *,'s = ',s,' stest = ',stest,' om = ',om,' omtest = ',omtest Print *,' l = ',l,' ltest = ',ltest Stop END