program geterror c c get the error to investigate order of spatial accuracy c (hardwired specifically for this case) c dimension cd(5),dh(5) c c first, write dh vs. error to show slopes c 1st order: open(11,file='slope1.dat',form='formatted',status='unknown') dh1=1.e-3 error1=1.e-6 dh2=1. error2=1.e-3 write(11,'(''variables="dh","error"'')') write(11,'(2e17.7)') dh1,error1 write(11,'(2e17.7)') dh2,error2 c 2nd order: open(12,file='slope2.dat',form='formatted',status='unknown') dh1=1.e-3 error1=1.e-6 dh2=1. error2=1. write(12,'(''variables="dh","error"'')') write(12,'(2e17.7)') dh1,error1 write(12,'(2e17.7)') dh2,error2 c c next, read the final results from each of the runs with different dh open(13,file='n0012_65t0.res',form='formatted',status='old') open(14,file='n0012_129t0.res',form='formatted',status='old') open(15,file='n0012_257t0.res',form='formatted',status='old') open(16,file='n0012_513t0.res',form='formatted',status='old') open(17,file='n0012_1025t0.res',form='formatted',status='old') dh(1)=sqrt(1./(65.*33.)) dh(2)=sqrt(1./(129.*65.)) dh(3)=sqrt(1./(257.*129.)) dh(4)=sqrt(1./(513.*257.)) dh(5)=sqrt(1./(1025.*513.)) do n=13,17 read(n,*) read(n,*) read(n,*) read(n,*) read(n,*) num do m=1,num read(n,*) idum,dum,dum,cd(n-12),dum,dum enddo enddo c c use finest 2 results with Richardson extrapolation (assuming 2nd order) c to get an estimate for results with an infinitely small dh cdexact = 4.*cd(5)/3. - cd(4)/3. write(6,'('' extrapolated cdexact is '',e18.8)') cdexact c c finally, write out results open(19,file='cderror.dat',form='formatted',status='unknown') write(19,'(''variables="dh","error","percent"'')') do n=1,5 write(19,'(3e17.7)') dh(n),abs(cdexact-cd(n)), + abs(cdexact-cd(n))/cdexact enddo write(6,'('' output to cderror.dat'')') write(6,'('' (slope curves for reference in slope1.dat'', + '' and slope2.dat)'')') stop end