!clt3ten.f90やclt5ten.f90をgnuplot用にまとめる
! /1sou と /8sou のdatをまとめて、理論値を付記する。

implicit none
!---入力---
real    :: b1=24.d-3/2    ! m   :幅/2(x)
real    :: h=240.d-3      ! m   :全層の高さ(y)
real    :: ells=240.d-3   ! m   :短い長さ(z)
real    :: elll=4800.d-3  ! m   :長い長さ(z)
integer :: nell=19        !     :長さ変化の分割数
integer :: nx=5           !     :幅方向の分割数
integer :: ny=24          !     :高さ方向の分割数
real    :: p=1.d-3        ! MN  :荷重
real    :: ezz=3.5d3      ! GPa :ヤング率
real    :: gyz=5.0d3/15   ! GPa :せん断弾性係数
!---
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,ell,ve,vt
character  ccut*20,filename1*10
real(8)    step
nitimen=(nx+1)*(ny+1)
allocate(v1(nitimen,nell+1),d1(nell+1))
allocate(v8(nitimen,nell+1),d8(nell+1))

!初期化
v1=0.d0
d1=0.d0
v8=0.d0
d8=0.d0
do j=0,nell
   ell=ells+(elll-ells)/nell*j
   !#--出力ファイル名
   write(filename1,'(i10)') nint(ell*1.d3)
   open(7,file="1sou/5ten"//trim(adjustl(filename1))//".dat", status='old')
   open(8,file="8sou/5ten"//trim(adjustl(filename1))//".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
   close(7)
   close(8)
enddo
!変位の平均の計算
do j=1,nell+1
   d1(j)=sum(v1(:,j))/size(v1(:,j))
   d8(j)=sum(v8(:,j))/size(v8(:,j))
enddo
!gnuplot用のデータの出力
print*,'# ell/h   (v1-ve)/ve*100 (v8-ve)/ve*100 (v1-vt)/vt*100 (v8-vt)/vt*100 (v1-v8)/v8*100'
do j=0,nell
   ell=ells+(elll-ells)/nell*j
   !#--3点曲げ
   !ve=p*ell**3/(1.d0/4*ezz*b1*h**3)
   !vt=p*ell**3/(1.d0/4*ezz*b1*h**3)+p*ell/(5.d0/6*gyz*b1*h)!k=5/6としている
   !#--5点曲げ
   ve=7.d0*p*ell**3/(64.d0*ezz*b1*h**3)
   vt=7.d0*p*ell**3/(64.d0*ezz*b1*h**3)+73.d0*p*ell/(256.d0*5.d0/6*gyz*b1*h)!k=5/6としている
   print'(1pe13.5,1pe13.5,1pe13.5,1pe13.5,1pe13.5,1pe13.5)'&
       &,ell/h,(d1(j+1)-ve)/ve*1.d2,(d8(j+1)-ve)/ve*1.d2,(d1(j+1)-vt)/vt*1.d2&
       &,(d8(j+1)-vt)/vt*1.d2,(d1(j+1)-d8(j+1))/d8(j+1)*1.d2
enddo
end
