!clthaba.f90をまとめる
! /1sou と /8sou のdatをgnuplot用にまとめる。

implicit none
!---入力---
real    :: b1s=200.d-3/2    ! m   :短い幅/2(x)
real    :: b1l=4000.d-3/2   ! m   :長い幅/2(x)
integer :: nb1=19           !     :長さ変化の分割数
real    :: ell=2010.d-3/2   ! m   :長さ/2(z)
integer :: nx=100           !     :幅方向の分割数
integer :: nz=100           !     :長さ方向の分割数
real    :: saikax=200.d-3/2 ! m   :載荷板の幅/2(x)
real    :: saikaz=500.d-3/2 ! m   :載荷板の幅/2(z)
!---
real(8),allocatable :: v1(:,:) ! m  :1層の載荷方向(y)の節点変位
real(8),allocatable :: d1(:)   ! m  :1層の載荷方向(y)の平均変位
real(8),allocatable :: v8(:,:) ! m  :8層の載荷方向(y)の節点変位
real(8),allocatable :: d8(:)   ! m  :8層の載荷方向(y)の平均変位
integer    nitimen             !    :断面の節点数
integer    i,j,icut
real       rcut
character  ccut*20,filename1*10
real(8)    step,n1,n8,b1
allocate(v1((nx+1)*(nint(nz*saikaz/ell)+1),nb1+1),d1(nb1+1))
allocate(v8((nx+1)*(nint(nz*saikaz/ell)+1),nb1+1),d8(nb1+1))

!初期化
v1=0.d0
d1=0.d0
v8=0.d0
d8=0.d0
do j=0,nb1
   b1=b1s+(b1l-b1s)/nb1*j
   nitimen=(int(nx*saikax/b1)+1)*(int(nz*saikaz/ell)+1)
   !#--出力ファイル名
   write(filename1,'(i10)') nint(b1*1.d3)
   open(7,file="1/haba"//trim(adjustl(filename1))//"b.dat", status='old')
   open(8,file="8/haba"//trim(adjustl(filename1))//"b.dat", status='old')
   !datファイルのデータの読み飛ばし
   read(7,'()')
   read(8,'()')
   !datファイルの荷重の割合を読み込み
   read(7,*) ccut,ccut,ccut,ccut,ccut,ccut,ccut,ccut,ccut,step
   read(8,*) ccut,ccut,ccut,ccut,ccut,ccut,ccut,ccut,ccut,step
   read(7,'()')
   read(8,'()')
   !datファイルの節点変位の読み込み
   do i=1,nitimen
      read(7,*) icut,rcut,v1(i,j+1),rcut
      read(8,*) icut,rcut,v8(i,j+1),rcut
   enddo
   !変位の平均を計算
   d1(j+1)=sum(v1(:,j+1))/nitimen
   d8(j+1)=sum(v8(:,j+1))/nitimen
   close(7)
   close(8)
enddo
!gnuplot用のデータの出力
print*,'# b1/ell   (v1-v8)/v8'
do j=0,nb1
   b1=b1s+(b1l-b1s)/nb1*j
   print'(1pe13.5,1pe13.5)'&
       &,b1/ell,(d1(j+1)-d8(j+1))/d8(j+1)*1.d2
enddo
end
