implicit real*8(a-h,o-z)
dimension x(315500),y(315500),z(315500),e66(6,6),b(3,3)&
&,meshmoku(300000,8),meshkou(300000,8),meshkou2(300000,8),meshhogou(300000,8)&
&,narax(315500),naray(315500),naraz(315500)
b1=8.4d-1
b2=0.9d-2
b3=0.7d-2!補剛
b4=6.8d-2!補剛
h1=10.0d-3
h2=100.0d-3
h3=h1
h4=255.0d-3
h5=5.0d-3
h6=h1
h7=h2
h9=7.0d-3
h8=h1+h2+h3-h6-h7-h9
ell=3.5d0
p=1.000d-1!荷重
ellhaji=150.0d-3
ellhaji2=165.0d-3
ellitahan=50.0d-3
ellaida=630.0d-3
!xyz方向の要素分割数V
nx1=3
nx2=4
nx3=3
nx4=4
ny1=3
ny2=4
ny3=ny1
ny4=3
ny5=4
ny6=ny1
ny7=ny2
ny8=3
ny9=4
nx=nx1+2*nx2
nita=10
nzhaji=3!0-150
nzhaji2=3!150-315
nzitahan=4!
nzaida=4!間のぶんかつ
nza=nita*nzitahan*2+(nzhaji+nzhaji2)*2+nzaida*(nita-1)
nitimen=(nx+1)*(ny1*2+ny2+1)*2+(nx2+1)*2*(ny4-1)
ellall=(ellhaji+ellhaji2)*2.0d0+ellitahan*nita*2.0d0+ellaida*(nita-1)
!ヤング率の入力
ezz=7.694d3!z軸方向のヤング率
exx=ezz/25.0d0!x軸方向のヤング率
eyy=ezz/25.0d0!y軸方向のヤング率
ezz2=206.0d3!z軸方向のヤング率
!ポアソン比の入力
poixy=0.016d0
poixz=poixy
poiyx=poixy
poiyz=poixy
poizx=0.4d0
poizy=poizx
!せん断弾性係数の入力
gxy=ezz/15.0d0
gxz=gxy
gyz=gxy
!1つあたりの要素の長さ?
xn1=b1/nx1
xn2=b2/nx2
xn3=b3/nx3
xn4=b4/nx4
yn1=h1/ny1
yn2=h2/ny2
yn3=h3/ny3
yn4=h4/ny4
yn5=h5/ny5
yn6=h6/ny6
yn7=h7/ny7
yn8=h8/ny8
yn9=h9/ny9
zn=ell/nznaka
znhaji=ellhaji/nzhaji
znhaji2=ellhaji2/nzhaji2
znitahan=ellitahan/nzitahan
znaida=ellaida/nzaida
nitimen=(nx1+nx2*2+1)*(ny1+ny2+ny3+1)+2*(nx2+1)*(ny4-1)+(nx2+nx3+1)*(ny5)*2&
&+(nx1+2*nx2+2*nx3+1)*(ny6+ny7+ny8)+(nx1+2*nx2+2*nx3+2*nx4+2*nx5+1)*(ny9+1)
! 節点番号ごとの節点座標!i~mで始まるのがdoで使える 
!した部分の節点（軸方向はじ）
!print'("*NODE, NSET=Nall")'
do k=1,nzhaji*2+nzhaji2*2+nzitahan*2+1
do j=1,ny1*2+ny2+1
do i=1,nx+1
if (k>=1.and.k<=nzhaji+nzhaji2+nzitahan) then
nset=i+(j-1)*(nx+1)+(k-1)*nitimen
else
nset=nitimen*(nza-nzhaji-nzhaji2-nzitahan-1)&
&+i+(j-1)*(nx+1)+(k-nzhaji-nzhaji2-nzitahan)*nitimen
end if
!
if (i>=1.and.i<=nx2+1) then
x(nset)=(i-1)*xn2+b4+b3
else if (i>=nx2+2.and.i<=nx2+nx1+1) then
x(nset)=(i-nx2-1)*xn1+b2+b4+b3
else 
x(nset)=(i-nx2-nx1-1)*xn2+b2+b1+b4+b3
end if
!
if (j>=1.and.j<=ny1+1) then
y(nset)=(j-1)*yn1
else if (j>=ny1+2.and.j<=ny2+ny1+1) then
y(nset)=(j-ny1-1)*yn2+h1
else 
y(nset)=(j-ny2-ny1-1)*yn1+h2+h1
end if
!
if (k>=1.and.k<=nzhaji+1) then
z(nset)=(k-1)*znhaji
else if (k>=nzhaji+2.and.k<=nzhaji+nzhaji2+1) then
z(nset)=(k-1-nzhaji)*znhaji2+ellhaji 
else if (k>=nzhaji+nzhaji2+2.and.k<=nzhaji+nzhaji2+nzitahan) then
z(nset)=(k-1-nzhaji-nzhaji2)*znitahan+ellhaji+ellhaji2
else if (k>=nzhaji+nzhaji2+nzitahan+1.and.k<=nzhaji+nzhaji2+nzitahan*2+1) then
z(nset)=(k-1-nzhaji-nzhaji2-nzitahan)*znitahan+ellall-ellhaji-ellhaji2-ellitahan
else if (k>=nzhaji+nzhaji2+nzitahan*2+2.and.k<=nzhaji+nzhaji2*2+nzitahan*2+1) then
z(nset)=(k-1-nzhaji-nzhaji2-nzitahan*2)*znhaji2+ellall-ellhaji-ellhaji2
else 
z(nset)=(k-1-nzhaji-nzhaji2*2-nzitahan*2)*znhaji+ellall-ellhaji
end if
!print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
!上と下の節点(軸方向あいだ)
!print'("*NODE, NSET=Nall")'
do l=1,nita-1
do k=1,(nzitahan*2+nzaida)
do j=1,ny1*2+ny2+1
do i=1,nx+1
nset=i+(j-1)*(nx+1)+(k-1)*nitimen&
&+(l-1)*nitimen*(nzitahan*2+nzaida)+nitimen*(nzhaji+nzhaji2+nzitahan)
!
if (i>=1.and.i<=nx2+1) then
x(nset)=(i-1)*xn2+b4+b3
else if (i>=nx2+2.and.i<=nx2+nx1+1) then
x(nset)=(i-nx2-1)*xn1+b2+b4+b3
else 
x(nset)=(i-nx2-nx1-1)*xn2+b2+b1+b4+b3
end if
!
if (j>=1.and.j<=ny1+1) then
y(nset)=(j-1)*yn1
else if (j>=ny1+2.and.j<=ny2+ny1+1) then
y(nset)=(j-ny1-1)*yn2+h1
else 
y(nset)=(j-ny2-ny1-1)*yn1+h2+h1
end if
!
if (k>=1.and.k<=nzitahan+1) then
z(nset)=(k-1)*znitahan+ellhaji+ellhaji2+ellitahan&
&+(l-1)*(ellitahan*2.0d0+ellaida)
else if (k>=nzitahan+2.and.k<=nzitahan+nzaida+1) then
z(nset)=(k-1-nzitahan)*znaida+ellhaji+ellhaji2+ellitahan*2.0d0&
&+(l-1)*(ellitahan*2.0d0+ellaida)
else 
z(nset)=(k-1-nzitahan-nzaida)*znitahan+ellhaji+ellhaji2+ellitahan*2.0d0&
&+ellaida+(l-1)*(ellitahan*2.0d0+ellaida)
end if
!print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do
!
!120-375の間の節点(軸方向はじ)

do m=1,2
do k=1,nzhaji*2+nzhaji2*2+nzitahan*2+1
do j=1,ny4-1
do i=1,nx2+1
if (k>=1.and.k<=nzhaji+nzhaji2+nzitahan) then
nset=i+(j-1)*(nx2+1)*2+(k-1)*nitimen+(m-1)*(nx2+1)&
&+(nx+1)*(ny1*2+ny2+1)
else
nset=nitimen*(nza-nzhaji-nzhaji2-nzitahan-1)&
&+i+(j-1)*(nx2+1)*2+(k-nzhaji-nzhaji2-nzitahan)*nitimen+(m-1)*(nx2+1)&
&+(nx+1)*(ny1*2+ny2+1)
end if
!
x(nset)=xn2*(i-1)+(m-1)*(b1+b2)+b4+b3
!
y(nset)=yn4*j+h1*2.0d0+h2
!
if (k>=1.and.k<=nzhaji+1) then
z(nset)=(k-1)*znhaji
else if (k>=nzhaji+2.and.k<=nzhaji+nzhaji2+1) then
z(nset)=(k-1-nzhaji)*znhaji2+ellhaji 
else if (k>=nzhaji+nzhaji2+2.and.k<=nzhaji+nzhaji2+nzitahan) then
z(nset)=(k-1-nzhaji-nzhaji2)*znitahan+ellhaji+ellhaji2
else if (k>=nzhaji+nzhaji2+nzitahan+1.and.k<=nzhaji+nzhaji2+nzitahan*2+1) then
z(nset)=(k-1-nzhaji-nzhaji2-nzitahan)*znitahan+ellall-ellhaji-ellhaji2-ellitahan
else if (k>=nzhaji+nzhaji2+nzitahan*2+2.and.k<=nzhaji+nzhaji2*2+nzitahan*2+1) then
z(nset)=(k-1-nzhaji-nzhaji2-nzitahan*2)*znhaji2+ellall-ellhaji-ellhaji2
else 
z(nset)=(k-1-nzhaji-nzhaji2*2-nzitahan*2)*znhaji+ellall-ellhaji
end if
!print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do
!120-375の間の節点(軸方向あいだ)

do l=1,nita-1
do m=1,2
do k=1,(nzitahan*2+nzaida)
do j=1,ny4-1
do i=1,nx2+1
nset=i+(j-1)*(nx2+1)*2+(k-1)*nitimen+(m-1)*(nx2+1)&
&+(l-1)*nitimen*(nzitahan*2+nzaida)+nitimen*(nzhaji+nzhaji2+nzitahan)&
&+(nx+1)*(ny1*2+ny2+1)
!
x(nset)=xn2*(i-1)+(m-1)*(b1+b2)+b4+b3
!
y(nset)=yn4*j+h1*2.0d0+h2
!
if (k>=1.and.k<=nzitahan+1) then
z(nset)=(k-1)*znitahan+ellhaji+ellhaji2+ellitahan&
&+(l-1)*(ellitahan*2.0d0+ellaida)
else if (k>=nzitahan+2.and.k<=nzitahan+nzaida+1) then
z(nset)=(k-1-nzitahan)*znaida+ellhaji+ellhaji2+ellitahan*2.0d0&
&+(l-1)*(ellitahan*2.0d0+ellaida)
else 
z(nset)=(k-1-nzitahan-nzaida)*znitahan+ellhaji+ellhaji2+ellitahan*2.0d0&
&+ellaida+(l-1)*(ellitahan*2.0d0+ellaida)
end if
!print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do
end do
!375-380の間の節点(軸方向はじ)
do k=1,nzhaji*2+nzhaji2*2+nzitahan*2+1
do j=1,ny5
do i=1,(nx3+nx2+1)*2
if (k>=1.and.k<=nzhaji+nzhaji2+nzitahan) then
nset=i+(j-1)*(nx3+nx2+1)*2+(k-1)*(nitimen)&
&+(nx+1)*(2*ny1+ny2+1)+(nx2+1)*(ny4-1)*2
else
nset=nitimen*(nza-nzhaji-nzhaji2-nzitahan-1)&
&+i+(j-1)*(nx3+nx2+1)*2+(k-nzhaji-nzhaji2-nzitahan)*nitimen&
&+(nx+1)*(2*ny1+ny2+1)+(nx2+1)*(ny4-1)*2
end if
!
if (i>=1.and.i<=nx3+1)then
x(nset)=b4+xn3*(i-1)
else if (i>=nx3+2.and.i<=nx3+nx2+1)then
x(nset)=b4+b3+xn2*(i-(nx3+1))
else if (i>=nx3+nx2+2.and.i<=nx3+nx2+1+nx2+1)then
x(nset)=b4+b3+b2+b1+xn2*(i-(nx3+nx2+2))
else 
x(nset)=b4+b3+b2+b1+b2+xn3*(i-(nx3+nx2+1+nx2+1))
end if

y(nset)=yn5*(j-1)+h1*2.0d0+h2+h4
!
if (k>=1.and.k<=nzhaji+1) then
z(nset)=(k-1)*znhaji
else if (k>=nzhaji+2.and.k<=nzhaji+nzhaji2+1) then
z(nset)=(k-1-nzhaji)*znhaji2+ellhaji 
else if (k>=nzhaji+nzhaji2+2.and.k<=nzhaji+nzhaji2+nzitahan) then
z(nset)=(k-1-nzhaji-nzhaji2)*znitahan+ellhaji+ellhaji2
else if (k>=nzhaji+nzhaji2+nzitahan+1.and.k<=nzhaji+nzhaji2+nzitahan*2+1) then
z(nset)=(k-1-nzhaji-nzhaji2-nzitahan)*znitahan+ellall-ellhaji-ellhaji2-ellitahan
else if (k>=nzhaji+nzhaji2+nzitahan*2+2.and.k<=nzhaji+nzhaji2*2+nzitahan*2+1) then
z(nset)=(k-1-nzhaji-nzhaji2-nzitahan*2)*znhaji2+ellall-ellhaji-ellhaji2
else 
z(nset)=(k-1-nzhaji-nzhaji2*2-nzitahan*2)*znhaji+ellall-ellhaji
end if
!print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
!375-380の間の節点(軸方向あいだ)
do l=1,nita-1
do k=1,(nzitahan*2+nzaida)
do j=1,ny5
do i=1,(nx3+nx2+1)*2
nset=i+(j-1)*(nx3+nx2+1)*2+(k-1)*nitimen&
&+(l-1)*nitimen*(nzitahan*2+nzaida)+nitimen*(nzhaji+nzhaji2+nzitahan)&
&+(nx+1)*(2*ny1+ny2+1)+(nx2+1)*(ny4-1)*2
!
if (i>=1.and.i<=nx3+1)then
x(nset)=b4+xn3*(i-1)
else if (i>=nx3+2.and.i<=nx3+nx2+1)then
x(nset)=b4+b3+xn2*(i-(nx3+1))
else if (i>=nx3+nx2+2.and.i<=nx3+nx2+1+nx2+1)then
x(nset)=b4+b3+b2+b1+xn2*(i-(nx3+nx2+2))
else 
x(nset)=b4+b3+b2+b1+b2+xn3*(i-(nx3+nx2+1+nx2+1))
end if
!
y(nset)=yn5*(j-1)+h1*2.0d0+h2+h4
!
if (k>=1.and.k<=nzitahan+1) then
z(nset)=(k-1)*znitahan+ellhaji+ellhaji2+ellitahan&
&+(l-1)*(ellitahan*2.0d0+ellaida)
else if (k>=nzitahan+2.and.k<=nzitahan+nzaida+1) then
z(nset)=(k-1-nzitahan)*znaida+ellhaji+ellhaji2+ellitahan*2.0d0&
&+(l-1)*(ellitahan*2.0d0+ellaida)
else 
z(nset)=(k-1-nzitahan-nzaida)*znitahan+ellhaji+ellhaji2+ellitahan*2.0d0&
&+ellaida+(l-1)*(ellitahan*2.0d0+ellaida)
end if
!print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do
!380-493の間の節点(軸方向はじ）
do k=1,nzhaji*2+nzhaji2*2+nzitahan*2+1
do j=1,ny6+ny7+ny8
do i=1,(nx3*2+nx2*2+nx1+1)
if (k>=1.and.k<=nzhaji+nzhaji2+nzitahan) then
nset=i+(j-1)*(nx3*2+nx2*2+nx1+1)+(k-1)*(nitimen)&
&+(nx+1)*(2*ny1+ny2+1)+(nx2+1)*(ny4-1)*2+(nx3+nx2+1)*(ny5)*2
else
nset=nitimen*(nza-nzhaji-nzhaji2-nzitahan-1)&
&+i+(j-1)*(nx3*2+nx2*2+nx1+1)+(k-nzhaji-nzhaji2-nzitahan)*nitimen&
&+(nx+1)*(2*ny1+ny2+1)+(nx2+1)*(ny4-1)*2+(nx3+nx2+1)*(ny5)*2
end if
!
if (i>=1.and.i<=nx3+1)then
x(nset)=b4+xn3*(i-1)
else if (i>=nx3+2.and.i<=nx3+nx2+1)then
x(nset)=b4+b3+xn2*(i-(nx3+1))
else if (i>=nx3+nx2+2.and.i<=nx3+nx2+nx1+1)then
x(nset)=b4+b3+b2+xn1*(i-(nx3+nx2+1))
else if (i>=nx3+nx2+nx1+2.and.i<=nx3+nx2+nx1+nx2+1)then
x(nset)=b4+b3+b2+b1+xn2*(i-(nx3+nx2+nx1+1))
else 
x(nset)=b4+b3+b2+b1+b2+xn3*(i-(nx3+nx2+nx1+1+nx2))
end if
!
if (j>=1.and.j<=ny6+1)then
y(nset)=yn6*(j-1)+h1*2.0d0+h2+h4+h5
else if (j>=ny6+2.and.j<=ny6+ny7+1)then
y(nset)=yn7*(j-ny6-1)+h1*2+h2+h4+h5+h6
else 
y(nset)=yn8*(j-ny6-ny7-1)+h1*2+h2+h4+h5+h6+h7
end if 
!
if (k>=1.and.k<=nzhaji+1) then
z(nset)=(k-1)*znhaji
else if (k>=nzhaji+2.and.k<=nzhaji+nzhaji2+1) then
z(nset)=(k-1-nzhaji)*znhaji2+ellhaji 
else if (k>=nzhaji+nzhaji2+2.and.k<=nzhaji+nzhaji2+nzitahan) then
z(nset)=(k-1-nzhaji-nzhaji2)*znitahan+ellhaji+ellhaji2
else if (k>=nzhaji+nzhaji2+nzitahan+1.and.k<=nzhaji+nzhaji2+nzitahan*2+1) then
z(nset)=(k-1-nzhaji-nzhaji2-nzitahan)*znitahan+ellall-ellhaji-ellhaji2-ellitahan
else if (k>=nzhaji+nzhaji2+nzitahan*2+2.and.k<=nzhaji+nzhaji2*2+nzitahan*2+1) then
z(nset)=(k-1-nzhaji-nzhaji2-nzitahan*2)*znhaji2+ellall-ellhaji-ellhaji2
else 
z(nset)=(k-1-nzhaji-nzhaji2*2-nzitahan*2)*znhaji+ellall-ellhaji
end if
!print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
!380-493の間の節点(軸方向あいだ)
do l=1,nita-1
do k=1,(nzitahan*2+nzaida)
do j=1,ny6+ny7+ny8
do i=1,(nx3*2+nx2*2+nx1+1)
nset=i+(j-1)*(nx3*2+nx2*2+nx1+1)+(k-1)*(nitimen)&
&+(nx+1)*(2*ny1+ny2+1)+(nx2+1)*(ny4-1)*2+(nx3+nx2+1)*(ny5)*2&
&+(l-1)*nitimen*(nzitahan*2+nzaida)+nitimen*(nzhaji+nzhaji2+nzitahan)!

if (i>=1.and.i<=nx3+1)then
x(nset)=b4+xn3*(i-1)
else if (i>=nx3+2.and.i<=nx3+nx2+1)then
x(nset)=b4+b3+xn2*(i-(nx3+1))
else if (i>=nx3+nx2+2.and.i<=nx3+nx2+nx1+1)then
x(nset)=b4+b3+b2+xn1*(i-(nx3+nx2+1))
else if (i>=nx3+nx2+nx1+2.and.i<=nx3+nx2+nx1+nx2+1)then
x(nset)=b4+b3+b2+b1+xn2*(i-(nx3+nx2+nx1+1))
else 
x(nset)=b4+b3+b2+b1+b2+xn3*(i-(nx3+nx2+nx1+1+nx2))
end if
!
if (j>=1.and.j<=ny6+1)then
y(nset)=yn6*(j-1)+h1*2.0d0+h2+h4+h5
else if (j>=ny6+2.and.j<=ny6+ny7+1)then
y(nset)=yn7*(j-ny6-1)+h1*2.0d0+h2+h4+h5+h6
else 
y(nset)=yn8*(j-ny6-ny7-1)+h1*2.0d0+h2+h4+h5+h6+h7
end if
if (k>=1.and.k<=nzitahan+1) then
z(nset)=(k-1)*znitahan+ellhaji+ellhaji2+ellitahan&
&+(l-1)*(ellitahan*2.0d0+ellaida)
else if (k>=nzitahan+2.and.k<=nzitahan+nzaida+1) then
z(nset)=(k-1-nzitahan)*znaida+ellhaji+ellhaji2+ellitahan*2.0d0&
&+(l-1)*(ellitahan*2.0d0+ellaida)
else 
z(nset)=(k-1-nzitahan-nzaida)*znitahan+ellhaji+ellhaji2+ellitahan*2.0d0&
&+ellaida+(l-1)*(ellitahan*2.0d0+ellaida)
end if
!
!print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do
!493-500の間の節点(軸方向はじ）
do k=1,nzhaji*2+nzhaji2*2+nzitahan*2+1
do j=1,ny9+1
do i=1,(nx4*2+nx3*2+nx2*2+nx1+1)
if (k>=1.and.k<=nzhaji+nzhaji2+nzitahan) then
nset=i+(j-1)*(nx4*2+nx3*2+nx2*2+nx1+1)+(k-1)*(nitimen)&
&+(nx+1)*(2*ny1+ny2+1)+(nx2+1)*(ny4-1)*2+(nx3+nx2+1)*(ny5)*2&
&+(nx3*2+nx+1)*(ny6+ny7+ny8)
else
nset=nitimen*(nza-nzhaji-nzhaji2-nzitahan-1)&
&+i+(j-1)*(nx4*2+nx3*2+nx2*2+nx1+1)+(k-nzhaji-nzhaji2-nzitahan)*nitimen&
&+(nx+1)*(2*ny1+ny2+1)+(nx2+1)*(ny4-1)*2+(nx3+nx2+1)*(ny5)*2&
&+(nx3*2+nx+1)*(ny6+ny7+ny8)
end if
!

if (i>=1.and.i<=nx4+1)then
x(nset)=xn4*(i-1)
else if (i>=nx4+2.and.i<=nx4+nx3+1)then
x(nset)=b4+xn3*(i-(nx4+1))
else if (i>=nx4+nx3+2.and.i<=nx4+nx3+nx2+1)then
x(nset)=b4+b3+xn2*(i-(nx4+nx3+1))
else if (i>=nx4+nx3+nx2+2.and.i<=nx4+nx3+nx2+nx1+1)then
x(nset)=b4+b3+b2+xn1*(i-(nx4+nx3+nx2+1))
else if (i>=nx4+nx3+nx2+nx1+2.and.i<=nx4+nx3+nx2+nx1+nx2+1)then
x(nset)=b4+b3+b2+b1+xn2*(i-(nx4+nx3+nx2+nx1+1))
else if (i>=nx4+nx3+nx2+nx1+nx2+2.and.i<=nx4+nx3+nx2+nx1+nx2+nx3+1)then
x(nset)=b4+b3+b2+b1+b2+xn3*(i-(nx4+nx3+nx2+nx1+nx2+1))
else 
x(nset)=b4+b3+b2+b1+b2+b3+xn4*(i-(nx4+nx3+nx2+nx1+nx2+nx3+1))
end if
!
y(nset)=yn9*(j-1)+h1*2+h2+h4+h5+h6+h7+h8
!
if (k>=1.and.k<=nzhaji+1) then
z(nset)=(k-1)*znhaji
else if (k>=nzhaji+2.and.k<=nzhaji+nzhaji2+1) then
z(nset)=(k-1-nzhaji)*znhaji2+ellhaji 
else if (k>=nzhaji+nzhaji2+2.and.k<=nzhaji+nzhaji2+nzitahan) then
z(nset)=(k-1-nzhaji-nzhaji2)*znitahan+ellhaji+ellhaji2
else if (k>=nzhaji+nzhaji2+nzitahan+1.and.k<=nzhaji+nzhaji2+nzitahan*2+1) then
z(nset)=(k-1-nzhaji-nzhaji2-nzitahan)*znitahan+ellall-ellhaji-ellhaji2-ellitahan
else if (k>=nzhaji+nzhaji2+nzitahan*2+2.and.k<=nzhaji+nzhaji2*2+nzitahan*2+1) then
z(nset)=(k-1-nzhaji-nzhaji2-nzitahan*2)*znhaji2+ellall-ellhaji-ellhaji2
else 
z(nset)=(k-1-nzhaji-nzhaji2*2-nzitahan*2)*znhaji+ellall-ellhaji
end if
!print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
!493-500の間の節点(軸方向あいだ)
do l=1,nita-1
do k=1,(nzitahan*2+nzaida)
do j=1,ny9+1
do i=1,(nx4*2+nx3*2+nx2*2+nx1+1)
nset=i+(j-1)*(nx4*2+nx3*2+nx2*2+nx1+1)+(k-1)*(nitimen)&
&+(nx+1)*(2*ny1+ny2+1)+(nx2+1)*(ny4-1)*2+(nx3+nx2+1)*(ny5)*2&
&+(nx3*2+nx+1)*(ny6+ny7+ny8)&
&+(l-1)*nitimen*(nzitahan*2+nzaida)+nitimen*(nzhaji+nzhaji2+nzitahan)!

if (i>=1.and.i<=nx4+1)then
x(nset)=xn4*(i-1)
else if (i>=nx4+2.and.i<=nx4+nx3+1)then
x(nset)=b4+xn3*(i-(nx4+1))
else if (i>=nx4+nx3+2.and.i<=nx4+nx3+nx2+1)then
x(nset)=b4+b3+xn2*(i-(nx4+nx3+1))
else if (i>=nx4+nx3+nx2+2.and.i<=nx4+nx3+nx2+nx1+1)then
x(nset)=b4+b3+b2+xn1*(i-(nx4+nx3+nx2+1))
else if (i>=nx4+nx3+nx2+nx1+2.and.i<=nx4+nx3+nx2+nx1+nx2+1)then
x(nset)=b4+b3+b2+b1+xn2*(i-(nx4+nx3+nx2+nx1+1))
else if (i>=nx4+nx3+nx2+nx1+nx2+2.and.i<=nx4+nx3+nx2+nx1+nx2+nx3+1)then
x(nset)=b4+b3+b2+b1+b2+xn3*(i-(nx4+nx3+nx2+nx1+nx2+1))
else 
x(nset)=b4+b3+b2+b1+b2+b3+xn4*(i-(nx4+nx3+nx2+nx1+nx2+nx3+1))
end if
!
y(nset)=yn9*(j-1)+h1*2+h2+h4+h5+h6+h7+h8
if (k>=1.and.k<=nzitahan+1) then
z(nset)=(k-1)*znitahan+ellhaji+ellhaji2+ellitahan&
&+(l-1)*(ellitahan*2.0d0+ellaida)
else if (k>=nzitahan+2.and.k<=nzitahan+nzaida+1) then
z(nset)=(k-1-nzitahan)*znaida+ellhaji+ellhaji2+ellitahan*2.0d0&
&+(l-1)*(ellitahan*2.0d0+ellaida)
else 
z(nset)=(k-1-nzitahan-nzaida)*znitahan+ellhaji+ellhaji2+ellitahan*2.0d0&
&+ellaida+(l-1)*(ellitahan*2.0d0+ellaida)
end if

!print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
!&nset,',',x(nset),',',y(nset),',',z(nset)
end do
end do
end do
end do 
!do i=1,nitimen*(nza+1)
!i5~6が整数,aが文字,f9,6が小数.1peとかはEの＊乗
!接点番号と座標の出力
nsettensuu=(nza+1)*nitimen

!print'("*NODE, NSET=Nall")'
!do i=1,nsettensuu
!print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
!&i,',',x(i),',',y(i),',',z(i)
!end do

!ここから
!要素番号とその要素の8接点の接点番号
!下部分の木材
nyou=0
do k=0,nza-1
do j=1,ny1*2+ny2
do i=1,nx1
nyou=nyou+1
meshmoku(nyou,1)=nx2+i+(j-1)*(2*nx2+nx1+1)&
&+k*nitimen
meshmoku(nyou,2)=meshmoku(nyou,1)+nx+1
meshmoku(nyou,3)=meshmoku(nyou,2)+nitimen
meshmoku(nyou,4)=meshmoku(nyou,1)+nitimen
meshmoku(nyou,5)=meshmoku(nyou,1)+1
meshmoku(nyou,6)=meshmoku(nyou,2)+1
meshmoku(nyou,7)=meshmoku(nyou,3)+1
meshmoku(nyou,8)=meshmoku(nyou,4)+1
end do
end do
end do
!!上の段の木材要素
do k=0,nza-1
do j=1,ny6+ny7+ny8-1
do i=1,nx1
nyou=nyou+1
meshmoku(nyou,1)=(nx+1)*(2*ny1+ny2+1)+(nx2+1)*(ny4-1)*2+(nx3+nx2+1)*(ny5)*2&
&+nx3+nx2+i+(j-1)*(nx+nx3*2+1)&
&+k*nitimen
meshmoku(nyou,2)=meshmoku(nyou,1)+(nx+nx3*2+1)
meshmoku(nyou,3)=meshmoku(nyou,2)+nitimen
meshmoku(nyou,4)=meshmoku(nyou,1)+nitimen
meshmoku(nyou,5)=meshmoku(nyou,1)+1
meshmoku(nyou,6)=meshmoku(nyou,2)+1
meshmoku(nyou,7)=meshmoku(nyou,3)+1
meshmoku(nyou,8)=meshmoku(nyou,4)+1
end do
end do
end do
do k=0,nza-1
do i=1,nx1
nyou=nyou+1
meshmoku(nyou,1)=(nx+1)*(2*ny1+ny2+1)+(nx2+1)*(ny4-1)*2+(nx3+nx2+1)*(ny5)*2&
&+nx3+nx2+i+(ny6+ny7+ny8-1)*(nx+nx3*2+1)&
&+k*nitimen
meshmoku(nyou,2)=(nx+1)*(2*ny1+ny2+1)+(nx2+1)*(ny4-1)*2+(nx3+nx2+1)*(ny5)*2&
&+(ny6+ny7+ny8)*(nx+nx3*2+1)+i+nx4+nx3+nx2&
&+k*nitimen
meshmoku(nyou,3)=meshmoku(nyou,2)+nitimen
meshmoku(nyou,4)=meshmoku(nyou,1)+nitimen
meshmoku(nyou,5)=meshmoku(nyou,1)+1
meshmoku(nyou,6)=meshmoku(nyou,2)+1
meshmoku(nyou,7)=meshmoku(nyou,3)+1
meshmoku(nyou,8)=meshmoku(nyou,4)+1
end do
end do
do k=0,nza-1
do j=1,ny9
do i=1,nx1
nyou=nyou+1
meshmoku(nyou,1)=(nx+1)*(2*ny1+ny2+1)+(nx2+1)*(ny4-1)*2+(nx3+nx2+1)*(ny5)*2&
&+(ny6+ny7+ny8)*(nx+nx3*2+1)+i+nx4+nx3+nx2+(j-1)*(nx4*2+nx3*2+nx+1)&
&+k*nitimen
meshmoku(nyou,2)=meshmoku(nyou,1)+nx4*2+nx3*2+nx+1
meshmoku(nyou,3)=meshmoku(nyou,2)+nitimen
meshmoku(nyou,4)=meshmoku(nyou,1)+nitimen
meshmoku(nyou,5)=meshmoku(nyou,1)+1
meshmoku(nyou,6)=meshmoku(nyou,2)+1
meshmoku(nyou,7)=meshmoku(nyou,3)+1
meshmoku(nyou,8)=meshmoku(nyou,4)+1
end do
end do
end do
!下部分の鋼板
nyou=0
do l=0,1
do k=0,nza-1 
do j=1,ny1*2+ny2
do i=1,nx2
nyou=nyou+1
meshkou(nyou,1)=i+(j-1)*(nx+1)&
&+k*nitimen&
&+l*(nx2+nx1)
meshkou(nyou,2)=meshkou(nyou,1)+nx+1
meshkou(nyou,3)=meshkou(nyou,2)+nitimen
meshkou(nyou,4)=meshkou(nyou,1)+nitimen
meshkou(nyou,5)=meshkou(nyou,1)+1
meshkou(nyou,6)=meshkou(nyou,2)+1
meshkou(nyou,7)=meshkou(nyou,3)+1
meshkou(nyou,8)=meshkou(nyou,4)+1
end do
end do 
end do
end do
!真ん中の鋼板(最初の一段目)
do m=0,1
do k=0,nza-1
do i=1,nx2
nyou=nyou+1
meshkou(nyou,1)=i+(nx+1)*(ny1*2+ny2)&
&+k*nitimen&
&+m*(nx2+nx1)
meshkou(nyou,2)=(nx+1)*(ny1*2+ny2+1)+i&
&+k*nitimen&
&+m*(nx2+1)
meshkou(nyou,3)=meshkou(nyou,2)+nitimen
meshkou(nyou,4)=meshkou(nyou,1)+nitimen
meshkou(nyou,5)=meshkou(nyou,1)+1
meshkou(nyou,6)=meshkou(nyou,2)+1
meshkou(nyou,7)=meshkou(nyou,3)+1
meshkou(nyou,8)=meshkou(nyou,4)+1
end do
end do 
end do
!真ん中の鋼板(2段目から最後の前まで)
do m=0,1
do k=0,nza-1
do j=1,ny4-2!ny4>2
do i=1,nx2
nyou=nyou+1
meshkou(nyou,1)=(nx+1)*(ny1*2+ny2+1)+i&
&+k*nitimen&
&+m*(nx2+1)+(j-1)*(nx2+1)*2
meshkou(nyou,2)=meshkou(nyou,1)+(nx2+1)*2
meshkou(nyou,3)=meshkou(nyou,2)+nitimen
meshkou(nyou,4)=meshkou(nyou,1)+nitimen
meshkou(nyou,5)=meshkou(nyou,1)+1
meshkou(nyou,6)=meshkou(nyou,2)+1
meshkou(nyou,7)=meshkou(nyou,3)+1
meshkou(nyou,8)=meshkou(nyou,4)+1
end do
end do 
end do
end do
!
!真ん中の鋼板(最後の段)
do m=0,1
do k=0,nza-1
do i=1,nx2
nyou=nyou+1
meshkou(nyou,1)=(nx+1)*(ny1*2+ny2+1)+(nx2+1)*(ny4-2)*2+i&
&+k*nitimen&
&+m*(nx2+1)
meshkou(nyou,2)=(nx+1)*(ny1*2+ny2+1)+(nx2+1)*(ny4-1)*2+nx3+i&
&+k*nitimen&
&+m*(nx2+1)
meshkou(nyou,3)=meshkou(nyou,2)+nitimen
meshkou(nyou,4)=meshkou(nyou,1)+nitimen
meshkou(nyou,5)=meshkou(nyou,1)+1
meshkou(nyou,6)=meshkou(nyou,2)+1
meshkou(nyou,7)=meshkou(nyou,3)+1
meshkou(nyou,8)=meshkou(nyou,4)+1
end do 
end do
end do
!下から3番目の鋼材
do m=0,1
do k=0,nza-1
do j=1,ny5
do i=1,nx2
nyou=nyou+1
if (j>=1.and.j<=ny5-1) then
meshkou(nyou,1)=(nx+1)*(ny1*2+ny2+1)+(nx2+1)*(ny4-1)*2&
&+nx3+i+(j-1)*(nx3+nx2+1)*2&
&+k*nitimen+m*(nx2+1)
meshkou(nyou,2)=meshkou(nyou,1)+(nx3+nx2+1)*2
meshkou(nyou,3)=meshkou(nyou,2)+nitimen
meshkou(nyou,4)=meshkou(nyou,1)+nitimen
meshkou(nyou,5)=meshkou(nyou,1)+1
meshkou(nyou,6)=meshkou(nyou,2)+1
meshkou(nyou,7)=meshkou(nyou,3)+1
meshkou(nyou,8)=meshkou(nyou,4)+1
else 
meshkou(nyou,1)=(nx+1)*(ny1*2+ny2+1)+(nx2+1)*(ny4-1)*2&
&+(nx3+nx2+1)*2*(ny5-1)+nx3+i+k*nitimen+m*(nx2+1)
meshkou(nyou,2)=(nx+1)*(ny1*2+ny2+1)+(nx2+1)*(ny4-1)*2&
&+(nx3+nx2+1)*2*(ny5)+nx3+i+k*nitimen+m*(nx2+nx1)
meshkou(nyou,3)=meshkou(nyou,2)+nitimen
meshkou(nyou,4)=meshkou(nyou,1)+nitimen
meshkou(nyou,5)=meshkou(nyou,1)+1
meshkou(nyou,6)=meshkou(nyou,2)+1
meshkou(nyou,7)=meshkou(nyou,3)+1
meshkou(nyou,8)=meshkou(nyou,4)+1
end if
end do
end do
end do
end do
!ここからmeshkou2 h6~h8の最後の段の1つまえまで
nyou=0
do m=0,1
do k=0,nza-1
do j=1,ny6+ny7+ny8-1
do i=1,nx2
nyou=nyou+1
meshkou2(nyou,1)=(nx+1)*(ny1*2+ny2+1)+2*(nx2+1)*(ny4-1)&
&+(nx3+nx2+1)*2*(ny5)+nx3+i+(j-1)*(nx+nx3*2+1)+k*nitimen&
&+m*(nx2+nx1)
meshkou2(nyou,2)=meshkou2(nyou,1)+nx+nx3*2+1
meshkou2(nyou,3)=meshkou2(nyou,2)+nitimen
meshkou2(nyou,4)=meshkou2(nyou,1)+nitimen
meshkou2(nyou,5)=meshkou2(nyou,1)+1
meshkou2(nyou,6)=meshkou2(nyou,2)+1
meshkou2(nyou,7)=meshkou2(nyou,3)+1
meshkou2(nyou,8)=meshkou2(nyou,4)+1
end do
end do
end do
end do
!h6~h8の最後の段
do m=0,1
do k=0,nza-1
do i=1,nx2
nyou=nyou+1
meshkou2(nyou,1)=(nx+1)*(ny1*2+ny2+1)+2*(nx2+1)*(ny4-1)&
&+(nx3+nx2+1)*2*(ny5)+nx3+i+(nx+nx3*2+1)*(ny6+ny7+ny8-1)+k*nitimen&
&+m*(nx2+nx1)
meshkou2(nyou,2)=(nx+1)*(ny1*2+ny2+1)+2*(nx2+1)*(ny4-1)&
&+(nx3+nx2+1)*2*(ny5)+(nx+nx3*2+1)*(ny6+ny7+ny8)+nx4+nx3+i+k*nitimen&
&+m*(nx2+nx1)
meshkou2(nyou,3)=meshkou2(nyou,2)+nitimen
meshkou2(nyou,4)=meshkou2(nyou,1)+nitimen
meshkou2(nyou,5)=meshkou2(nyou,1)+1
meshkou2(nyou,6)=meshkou2(nyou,2)+1
meshkou2(nyou,7)=meshkou2(nyou,3)+1
meshkou2(nyou,8)=meshkou2(nyou,4)+1
end do
end do
end do
!h9のだん全部
do m=0,1
do k=0,nza-1
do j=1,ny9
do i=1,nx2
nyou=nyou+1
meshkou2(nyou,1)=(nx+1)*(ny1*2+ny2+1)+2*(nx2+1)*(ny4-1)&
&+(nx3+nx2+1)*2*(ny5)+(nx+nx3*2+1)*(ny6+ny7+ny8)+nx4+nx3+i+k*nitimen&
&+m*(nx2+nx1)+(j-1)*(nx+nx3*2+nx4*2+1)
meshkou2(nyou,2)=meshkou2(nyou,1)+nx+nx3*2+nx4*2+1
meshkou2(nyou,3)=meshkou2(nyou,2)+nitimen
meshkou2(nyou,4)=meshkou2(nyou,1)+nitimen
meshkou2(nyou,5)=meshkou2(nyou,1)+1
meshkou2(nyou,6)=meshkou2(nyou,2)+1
meshkou2(nyou,7)=meshkou2(nyou,3)+1
meshkou2(nyou,8)=meshkou2(nyou,4)+1
end do
end do
end do
end do
!ここまでmeshkou2
!補剛材鋼板の隣（内部）h5の段の最後の段の1つ前の段まで
nyou=0
do m=0,1
do k=0,nza-1
do j=1,ny5-1
do i=1,nx3
nyou=nyou+1
meshhogou(nyou,1)=(nx+1)*(ny1*2+ny2+1)+2*(nx2+1)*(ny4-1)&
&+i+(j-1)*(nx3+nx2+1)*2+k*nitimen+m*(nx3+nx2*2+1)
meshhogou(nyou,2)=meshhogou(nyou,1)+(nx3+nx2+1)*2
meshhogou(nyou,3)=meshhogou(nyou,2)+nitimen
meshhogou(nyou,4)=meshhogou(nyou,1)+nitimen
meshhogou(nyou,5)=meshhogou(nyou,1)+1
meshhogou(nyou,6)=meshhogou(nyou,2)+1
meshhogou(nyou,7)=meshhogou(nyou,3)+1
meshhogou(nyou,8)=meshhogou(nyou,4)+1
end do
end do
end do
end do
do m=0,1
do k=0,nza-1
do i=1,nx3
nyou=nyou+1
meshhogou(nyou,1)=(nx+1)*(ny1*2+ny2+1)+2*(nx2+1)*(ny4-1)&
&+i+(nx3+nx2+1)*2*(ny5-1)+k*nitimen+m*(nx3+nx2*2+1)
meshhogou(nyou,2)=(nx+1)*(ny1*2+ny2+1)+2*(nx2+1)*(ny4-1)&
&+i+(nx3+nx2+1)*2*(ny5)+k*nitimen+m*(nx3+nx2+nx1+nx2)
meshhogou(nyou,3)=meshhogou(nyou,2)+nitimen
meshhogou(nyou,4)=meshhogou(nyou,1)+nitimen
meshhogou(nyou,5)=meshhogou(nyou,1)+1
meshhogou(nyou,6)=meshhogou(nyou,2)+1
meshhogou(nyou,7)=meshhogou(nyou,3)+1
meshhogou(nyou,8)=meshhogou(nyou,4)+1
end do
end do
end do
!h6~h8のだん最後の1つ前まで
do m=0,1
do k=0,nza-1
do j=1,ny6+ny7+ny8-1
do i=1,nx3
nyou=nyou+1
meshhogou(nyou,1)=(nx+1)*(ny1*2+ny2+1)+2*(nx2+1)*(ny4-1)&
&+(nx3+nx2+1)*2*(ny5)+i+(j-1)*(nx+nx3*2+1)+k*nitimen+m*(nx3+nx2+nx1+nx2)
meshhogou(nyou,2)=meshhogou(nyou,1)+nx+nx3*2+1
meshhogou(nyou,3)=meshhogou(nyou,2)+nitimen
meshhogou(nyou,4)=meshhogou(nyou,1)+nitimen
meshhogou(nyou,5)=meshhogou(nyou,1)+1
meshhogou(nyou,6)=meshhogou(nyou,2)+1
meshhogou(nyou,7)=meshhogou(nyou,3)+1
meshhogou(nyou,8)=meshhogou(nyou,4)+1
end do
end do
end do
end do
do m=0,1
do k=0,nza-1
do i=1,nx3
nyou=nyou+1
meshhogou(nyou,1)=(nx+1)*(ny1*2+ny2+1)+2*(nx2+1)*(ny4-1)&
&+(nx3+nx2+1)*2*(ny5)+i+(nx+nx3*2+1)*(ny6+ny7+ny8-1)+k*nitimen+m*(nx3+nx2+nx1+nx2)
meshhogou(nyou,2)=(nx+1)*(ny1*2+ny2+1)+2*(nx2+1)*(ny4-1)&
&+(nx3+nx2+1)*2*(ny5)+nx4+i+(nx+nx3*2+1)*(ny6+ny7+ny8)+k*nitimen+m*(nx3+nx2+nx1+nx2)
meshhogou(nyou,3)=meshhogou(nyou,2)+nitimen
meshhogou(nyou,4)=meshhogou(nyou,1)+nitimen
meshhogou(nyou,5)=meshhogou(nyou,1)+1
meshhogou(nyou,6)=meshhogou(nyou,2)+1
meshhogou(nyou,7)=meshhogou(nyou,3)+1
meshhogou(nyou,8)=meshhogou(nyou,4)+1
end do
end do
end do
!h9のだん全部
do m=0,1
do k=0,nza-1
do j=1,ny9
do i=1,nx3
nyou=nyou+1
meshhogou(nyou,1)=(nx+1)*(ny1*2+ny2+1)+2*(nx2+1)*(ny4-1)&
&+(nx3+nx2+1)*2*(ny5)+(nx+nx3*2+1)*(ny6+ny7+ny8)&
&+nx4+i+(j-1)*(nx4*2+nx3*2+nx+1)+k*nitimen+m*(nx3+nx2+nx1+nx2)
meshhogou(nyou,2)=meshhogou(nyou,1)+nx4*2+nx3*2+nx+1
meshhogou(nyou,3)=meshhogou(nyou,2)+nitimen
meshhogou(nyou,4)=meshhogou(nyou,1)+nitimen
meshhogou(nyou,5)=meshhogou(nyou,1)+1
meshhogou(nyou,6)=meshhogou(nyou,2)+1
meshhogou(nyou,7)=meshhogou(nyou,3)+1
meshhogou(nyou,8)=meshhogou(nyou,4)+1
end do
end do
end do
end do
!補剛材突き出し部(nx4の部分）
do m=0,1
do k=0,nza-1
do j=1,ny9
do i=1,nx4
nyou=nyou+1
meshhogou(nyou,1)=(nx+1)*(ny1*2+ny2+1)+2*(nx2+1)*(ny4-1)&
&+(nx2+nx3+1)*2*ny5+(nx+nx3*2+1)*(ny6+ny7+ny8)&
&+i+(j-1)*(nx4*2+nx+nx3*2+1)+k*nitimen+m*(nx4+nx3*2+nx)
meshhogou(nyou,2)=meshhogou(nyou,1)+nx+nx3*2+nx4*2+1
meshhogou(nyou,3)=meshhogou(nyou,2)+nitimen
meshhogou(nyou,4)=meshhogou(nyou,1)+nitimen
meshhogou(nyou,5)=meshhogou(nyou,1)+1
meshhogou(nyou,6)=meshhogou(nyou,2)+1
meshhogou(nyou,7)=meshhogou(nyou,3)+1
meshhogou(nyou,8)=meshhogou(nyou,4)+1
end do
end do
end do
end do
!
!節点増やすための作業
!naraに-1を入れる
nsetten=nitimen*(nza+1)
do i=1,nsetten
narax(i)=-1
naray(i)=-1
naraz(i)=-1
end do
!増やさなくていい節点を設定a!onsi3tenza.f90からコピペ
do j=1,nita-1
do i=1,nsetten
if&
&(z(i)>=ellhaji+ellhaji2+ellitahan+(j-1)*(ellitahan*2.0d0+ellaida)-0.000001d0&
&.and.z(i)<=ellhaji+ellhaji2+ellitahan+ellitahan*2.0d0+(j-1)*(ellitahan*2.0d0+ellaida)&
&+0.000001d0) then
naraz(i)=-2
end if
end do
end do
!改変
do i=1,nsetten
if (y(i)>=h1-0.000001d0.and.y(i)<=h1+h2+0.0000001d0) then
naray(i)=-2
else if (y(i)>=h6+h5+h4+h2+h1*2.0d0-0.000001d0&
&.and.y(i)<=h7+h6+h5+h4+h2+h1*2.0d0+0.0000001d0) then
naray(i)=-3
end if
!追加
if (x(i)>=b4+b3+b2-0.000001d0.and.x(i)<=b4+b3+b2+0.000001d0) then
narax(i)=-2
else if (x(i)>=b4+b3+b2+b1-0.000001d0.and.x(i)<=b4+b3+b2+b1+0.000001d0) then
narax(i)=-2
else if (x(i)>=b4+b3-0.000001d0.and.x(i)<=b4+b3+0.000001d0) then
narax(i)=-3
else if (x(i)>=b4+b3+b2+b1+b2-0.000001d0&
&.and.x(i)<=b4+b3+b2+b1+b2+0.000001d0) then
narax(i)=-3
end if
end do
!
do i=1,nsetten
if (naraz(i)==-2.and.naray(i)==-2.and.narax(i)==-2) then
narax(i)=-4
else if (naraz(i)==-2.and.naray(i)==-3.and.narax(i)==-2) then
narax(i)=-4
else if (naraz(i)==-2.and.naray(i)==-3.and.narax(i)==-3) then
narax(i)=-4
end if
end do
!b4+b3追加
!木材と鋼板の間の節点を増やす
iset=nsetten
do i=1,nx1*((ny1*2+ny2)+(ny6+ny7+ny8+ny9))*nza
if (narax(meshmoku(i,1))==-2) then
iset=iset+1
narax(meshmoku(i,1))=iset
x(iset)=b2+b4+b3
y(iset)=y(meshmoku(i,1))
z(iset)=z(meshmoku(i,1))
meshmoku(i,1)=iset
end if
if (narax(meshmoku(i,2))==-2) then
iset=iset+1
narax(meshmoku(i,2))=iset
x(iset)=b2+b4+b3
y(iset)=y(meshmoku(i,2))
z(iset)=z(meshmoku(i,2))
meshmoku(i,2)=iset
end if
if (narax(meshmoku(i,2))==-2) then
iset=iset+1
narax(meshmoku(i,3))=iset
x(iset)=b2+b4+b3
y(iset)=y(meshmoku(i,3))
z(iset)=z(meshmoku(i,3))
meshmoku(i,3)=iset
end if
if (narax(meshmoku(i,4))==-2) then
iset=iset+1
narax(meshmoku(i,4))=iset
x(iset)=b2+b4+b3
y(iset)=y(meshmoku(i,4))
z(iset)=z(meshmoku(i,4))
meshmoku(i,4)=iset
end if
if (narax(meshmoku(i,5))==-2) then
iset=iset+1
narax(meshmoku(i,5))=iset
x(iset)=b2+b1+b4+b3
y(iset)=y(meshmoku(i,5))
z(iset)=z(meshmoku(i,5))
meshmoku(i,5)=iset
end if
if (narax(meshmoku(i,6))==-2) then
iset=iset+1
narax(meshmoku(i,6))=iset
x(iset)=b2+b1+b4+b3
y(iset)=y(meshmoku(i,6))
z(iset)=z(meshmoku(i,6))
meshmoku(i,6)=iset
end if
if (narax(meshmoku(i,7))==-2) then
iset=iset+1
narax(meshmoku(i,7))=iset
x(iset)=b2+b1+b4+b3
y(iset)=y(meshmoku(i,7))
z(iset)=z(meshmoku(i,7))
meshmoku(i,7)=iset
end if
if (narax(meshmoku(i,8))==-2) then
iset=iset+1
narax(meshmoku(i,8))=iset
x(iset)=b2+b1+b4+b3
y(iset)=y(meshmoku(i,8))
z(iset)=z(meshmoku(i,8))
meshmoku(i,8)=iset
end if
end do
!b4+b3追加
!鋼板と補剛材の間の節点を増やす

do i=1,2*(nx4*ny9+nx3*(ny5+ny6+ny7+ny8+ny9))*nza
if (narax(meshhogou(i,1))==-3) then
iset=iset+1
narax(meshhogou(i,1))=iset
x(iset)=x(meshhogou(i,1))
y(iset)=y(meshhogou(i,1))
z(iset)=z(meshhogou(i,1))
meshhogou(i,1)=iset
end if
if (narax(meshhogou(i,2))==-3) then
iset=iset+1
narax(meshhogou(i,2))=iset
x(iset)=x(meshhogou(i,2))
y(iset)=y(meshhogou(i,2))
z(iset)=z(meshhogou(i,2))
meshhogou(i,2)=iset
end if
if (narax(meshhogou(i,2))==-3) then
iset=iset+1
narax(meshhogou(i,3))=iset
x(iset)=x(meshhogou(i,3))
y(iset)=y(meshhogou(i,3))
z(iset)=z(meshhogou(i,3))
meshhogou(i,3)=iset
end if
if (narax(meshhogou(i,4))==-3) then
iset=iset+1
narax(meshhogou(i,4))=iset
x(iset)=x(meshhogou(i,4))
y(iset)=y(meshhogou(i,4))
z(iset)=z(meshhogou(i,4))
meshhogou(i,4)=iset
end if
if (narax(meshhogou(i,5))==-3) then
iset=iset+1
narax(meshhogou(i,5))=iset
x(iset)=x(meshhogou(i,5))
y(iset)=y(meshhogou(i,5))
z(iset)=z(meshhogou(i,5))
meshhogou(i,5)=iset
end if
if (narax(meshhogou(i,6))==-3) then
iset=iset+1
narax(meshhogou(i,6))=iset
x(iset)=x(meshhogou(i,6))
y(iset)=y(meshhogou(i,6))
z(iset)=z(meshhogou(i,6))
meshhogou(i,6)=iset
end if
if (narax(meshhogou(i,7))==-3) then
iset=iset+1
narax(meshhogou(i,7))=iset
x(iset)=x(meshhogou(i,7))
y(iset)=y(meshhogou(i,7))
z(iset)=z(meshhogou(i,7))
meshhogou(i,7)=iset
end if
if (narax(meshhogou(i,8))==-3) then
iset=iset+1
narax(meshhogou(i,8))=iset
x(iset)=x(meshhogou(i,8))
y(iset)=y(meshhogou(i,8))
z(iset)=z(meshhogou(i,8))
meshhogou(i,8)=iset
end if
end do

print'("*NODE, NSET=Nall")'
do i=1,iset
print'(i6,a,1pe13.5,a,1pe13.5,a,1pe13.5)',&
&i,',',x(i),',',y(i),',',z(i)
end do


!設定
print'("*ELEMENT, TYPE=C3D8, ELSET=Ekou")'
do i=1,2*nx2*(ny1+ny2+ny3+ny4+ny5)*nza
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&i,',',&
&meshkou(i,1),',',&
&meshkou(i,2),',',&
&meshkou(i,3),',',&
&meshkou(i,4),',',&
&meshkou(i,5),',',&
&meshkou(i,6),',',&
&meshkou(i,7),',',&
&meshkou(i,8)
end do
nyousosuukou=2*nx2*(ny1+ny2+ny3+ny4+ny5)*nza
print'("*ELEMENT, TYPE=C3D8, ELSET=Ekou2")'
do i=1,2*nx2*(ny6+ny7+ny8+ny9)*nza
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&i+nyousosuukou,',',&
&meshkou2(i,1),',',&
&meshkou2(i,2),',',&
&meshkou2(i,3),',',&
&meshkou2(i,4),',',&
&meshkou2(i,5),',',&
&meshkou2(i,6),',',&
&meshkou2(i,7),',',&
&meshkou2(i,8)
end do
nyousosuukou2=2*nx2*(ny6+ny7+ny8+ny9)*nza
print'("*ELEMENT, TYPE=C3D8, ELSET=Ehogou")'
do i=1,2*(nx4*ny9+nx3*(ny5+ny6+ny7+ny8+ny9))*nza
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&i+nyousosuukou+nyousosuukou2,',',&
&meshhogou(i,1),',',&
&meshhogou(i,2),',',&
&meshhogou(i,3),',',&
&meshhogou(i,4),',',&
&meshhogou(i,5),',',&
&meshhogou(i,6),',',&
&meshhogou(i,7),',',&
&meshhogou(i,8)
end do
nyousosuuhogou=2*(nx4*ny9+nx3*(ny5+ny6+ny7+ny8+ny9))*nza
print'("*ELEMENT, TYPE=C3D8, ELSET=Emoku")'
do i=1,nx1*(ny1+ny2+ny3+ny6+ny7+ny8+ny9)*nza
print'(i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6,a,i6)',&
&i+nyousosuukou+nyousosuuhogou+nyousosuukou2,',',&
&meshmoku(i,1),',',&
&meshmoku(i,2),',',&
&meshmoku(i,3),',',&
&meshmoku(i,4),',',&
&meshmoku(i,5),',',&
&meshmoku(i,6),',',&
&meshmoku(i,7),',',&
&meshmoku(i,8)
end do

!載荷接点、拘束接点の設定!少し変えた
print'("*NSET,NSET=kousoku")'

do i=1,nitimen*(nza+1)
if (z(i)>=ellhaji-0.000001d0.and.z(i)<=ellhaji+0.00001d0&
&.and.y(i)<=0.000001d0.and.y(i)>=-0.00001d0) then
nkousoku=i
print'(i6,a)',nkousoku,','
end if
end do
print'("*NSET,NSET=kousoku2")'
do i=1,nitimen*(nza+1)
if (z(i)>=ellall-ellhaji-0.000001d0.and.z(i)<=ellall-ellhaji+0.00001d0&
&.and.y(i)<=0.000001d0.and.y(i)>=-0.00001d0) then
nkousoku=i
print'(i6,a)',nkousoku,','
end if
end do
print'("*NSET,NSET=saika")'
do i=1,nitimen*(nza+1)
if (z(i)>=ellall/2.0d0-0.00001d0.and.z(i)<=ellall/2.0d0+0.00001d0&
&.and.y(i)>=h1+h2+h3+h4+h5+h6+h7+h8+h9-0.000001d0&
&.and.y(i)<=h1+h2+h3+h4+h5+h6+h7+h8+h9+0.00001d0&
&.and.x(i)>=b4+b3+b2+xn1-0.00001d0&
&.and.x(i)<=b4+b3+b2+b1-xn1+0.00001d0) then
nsaika=i
print'(i6,a)',nsaika,','
end if
end do
! 載荷点の下の節点番号コピペ
!print'("*NSET,NSET=saikasita")'
print'("*NSET,NSET=saikasita")'
do i=1,iset
if (z(i)<=ellall/2.0d0+0.000001d0.and.z(i)>=ellall/2.0d0-0.000001d0&
&.and.y(i)<=0.000001d0.and.y(i)>=-0.000001d0) then
print'(i6,a)',i,','
end if 
end do
print'("*NSET,NSET=saikamen")'
do i=1,iset
if (z(i)<=ellall/2.0d0+0.000001d0.and.z(i)>=ellall/2.0d0-0.000001d0) then
print'(i6,a)',i,','
end if 
end do
! 境界条件
print'("*BOUNDARY")'
print'("kousoku,1,3")'
print'("kousoku2,1,2")'
! 木材の設定
print'("*MATERIAL,NAME=ELmoku")'
print'("*ELASTIC,TYPE=ORTHOROPIC")'
call ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
& gxz,gxy,poiyx,poizx,poizy)
print'(f12.4,a,f12.4,a,f12.4,a,&
&f12.4,a,f12.4,a,f12.4,a,f12.4,&
& a ,f12.4,a)',&
&e66(1,1),',',e66(1,2),',',e66(2,2),',',e66(1,3),','& 
&,e66(2,3),',',e66(3,3),',',e66(6,6),',',e66(5,5),','
print'(f12.4,a,a)',e66(4,4),',','293'
print'("*MATERIAL,NAME=ELkou")'
print'("*ELASTIC")'
print'("206.d3, 0.3")'
print'("*SOLID SECTION,ELSET=Emoku,MATERIAL=ELmoku")'
print'("*SOLID SECTION,ELSET=Ekou,MATERIAL=ELkou")'
print'("*SOLID SECTION,ELSET=Ekou2,MATERIAL=ELkou")'
print'("*SOLID SECTION,ELSET=Ehogou,MATERIAL=ELkou")'
! 載荷条件
print'("*STEP")'
print'("*STATIC,SOLVER=SPOOLES")'
print'("*BUCKLE")'
print'("3,0.01")'
print'("*CLOAD")'
print'(A,a,i6,a,1pe13.5)','saika',',',2,',',p/(nx+nx3*2+nx4*2+1)
! datのオプション
print'("*NODE PRINT,NSET=Nall")'
!print'("*NODE PRINT,NSET=saika")'
print'("U")'
print'("*NODE FILE")'
print'("U")'
print'("*END STEP")'
!


end
!######## 直交異方性の定数を出すプログラム ########
subroutine ortho(e66,exx,eyy,ezz,poixy,poixz,poiyz,gyz,&
&gxz,gxy,poiyx,poizx,poizy)
implicit real*8(a-h,o-z)
dimension e66(6,6),b(3,3)
!応力－ひずみ行列：e66
do i=1,6
do j=1,6
e66(i,j)=0.d0
end do
end do
!  
e66(1,1)=1.d0/exx
e66(1,2)=-poixy/exx
e66(1,3)=-poixz/exx
e66(2,1)=-poiyx/eyy
e66(2,2)=1.d0/eyy
e66(2,3)=-poiyz/eyy
e66(3,1)=-poizx/ezz
e66(3,2)=-poizy/ezz
e66(3,3)=1.d0/ezz
!
!****** 逆行列を取って、応力－ひずみ行列に変換する******
det=e66(1,1)*e66(2,2)*e66(3,3)+e66(2,1)*e66(3,2)*e66(1,3)&
&+e66(3,1)*e66(1,2)*e66(2,3)-e66(2,1)*e66(1,2)*e66(3,3)&
&-e66(3,1)*e66(2,2)*e66(1,3)-e66(1,1)*e66(3,2)*e66(2,3)
b(1,1)=(-1.d0)**(1+1)*(e66(2,2)*e66(3,3)-e66(3,2)*e66(2,3))
b(1,2)=(-1.d0)**(1+2)*(e66(2,1)*e66(3,3)-e66(3,1)*e66(2,3))
b(1,3)=(-1.d0)**(1+3)*(e66(2,1)*e66(3,2)-e66(3,1)*e66(2,2))
b(2,1)=(-1.d0)**(2+1)*(e66(1,2)*e66(3,3)-e66(3,2)*e66(1,3))
b(2,2)=(-1.d0)**(2+2)*(e66(1,1)*e66(3,3)-e66(3,1)*e66(1,3))
b(2,3)=(-1.d0)**(2+3)*(e66(1,1)*e66(3,2)-e66(3,1)*e66(1,2))
b(3,1)=(-1.d0)**(3+1)*(e66(1,2)*e66(2,3)-e66(2,2)*e66(1,3))
b(3,2)=(-1.d0)**(3+2)*(e66(1,1)*e66(2,3)-e66(2,1)*e66(1,3))
b(3,3)=(-1.d0)**(3+3)*(e66(1,1)*e66(2,2)-e66(2,1)*e66(1,2))
do i=1,3
do j=1,3
e66(i,j)=b(j,i)/det
end do
end do
!
!****** せん断補正係数を入力 ******
!print*,'xy面のせん断弾性係数を入れて下さい'
!read*, gxy
e66(4,4)=gyz
!print*,'xz面のせん断弾性係数を入れて下さい'
!read*, gxz
e66(5,5)=gxz
!print*,'yz面のせん断弾性係数を入れて下さい'
!read*, gyz
e66(6,6)=gxy
!
!
return
end



