      program zgousei3 ver1
      implicit real*8(a-h,o-z)
      dimension area(100,100),e66(6,6),nofn(8),b(3,3)
!
!######## 各断面の諸定数入力 ########
!      print*,'x軸(断面の横の)の長さを入れてください'
!      read*,zl
      zl=6.d-2
!
!      print*,'y軸(断面の縦の)の長さを入れてください'
!      read*,zm
      zm=12.d-2
!
!      print*,'z軸方向(軸方向)の長さを入れて下さい'
!      read*,ziku
      ziku=1.08d0
!
!      print*,'上部に挿入する鋼板のx軸の長さを入れてください'
!      read*,zlku
      zlku=9.d-3
!
!      print*,'上部に挿入する鋼板のy軸の長さを入れてください'
!      read*,zmku
      zmku=4.d-2
!
!      print*,'下部に挿入する鋼板のx軸の長さを入れて下さい'
!      read*,zlkl
      zlkl=4.5d-3
!      
!      print*,'下部に挿入する鋼板のy軸の長さを入れて下さい'
!      read*,zmkl
      zmkl=4.d-2
!
!      print*,'図心軸上以外の下部鋼と、図心軸からの距離を入れて下さい'
!      read*,zlkls
      zlkls=1.5d-2
! 集成材を内側から番号付け
      zlt1=(zlku-zlkl)/2.d0
      zlt2=zlkls-(zlkl+zlku)/2.d0
      zlt3=zlkls-zlkl/2.d0
!
!######## 要素数入力 ########
!      print*,'要素数入力は、細心の注意を払って下さい'
!      print*,'まず、x軸方向の各要素数を設定します'
!
!      print*,'x軸方向の全要素数を入れて下さい(最低12以上の偶数)'
!      read*,nx
      nx=12
!ちなみに、nx=3nxkl+2nxt1+2nxt2+2nxt3
!
!      print*,'上部の鋼板のx軸方向の要素数を入れて下さい(4以上の倍数)'
!      read*,nxku
      nxku=4
!
!      print*,'下部の鋼板のx軸方向の要素数を入れて下さい(偶数)
!      read*,nxkl
!      nxkl=2
      nxkl=nxku/2
! nxkuとnxklにより、自動的に上部鋼と下部鋼の差の分の集成材要素nxt1が決定
      nxt1=(nxku-nxkl)/2
!      
!      print*,'端の下部鋼の外側の集成材のx方向の要素数を入れて下さい'
!      read*,nxt3
      nxt3=1
! 以上から、上部鋼と下部鋼の間の集成材の要素が決定
      nxt2=(nx-2*nxt1-2*nxkl-nxku)/2
!
!      print*,'次に、y軸方向の各要素数を設定します'
!      print*,'y軸方向の全要素数を入れて下さい(最低4以上の偶数)'
!      read*,ny
!ちなみに、ny=nykl+nyku+nyt(注：nytは便宜上のものでプログラム上には無い)
      ny=10
!
!      print*,'上部の鋼材のy軸方向の要素数を入れて下さい'
!      read*,nyku
      nyku=4
!
!      print*,'上部の鋼材のy軸方向の要素数を入れて下さい'
!      read*,nykl
      nykl=4
!
!      print*,'最後に、z軸方向の全要素数を入れて下さい'
!      read*,nz
      nz=100
!
!######## 出力 #########
      open(7,file='zgousei3.inp')
!
      call zahyou(nx,ny,nz,zl,zm,ziku,nxt1,nxt2,nxt3,nxkl,nxku,zlku,
     &zlkl,zlt1,zlt2,zlt3,nyku,nykl,zmku,zmkl,nnofn,zax,zay,zaz)
!
      call youso(nx,ny,nz,nxt1,nxt2,nxt3,nxkl,nykl,nxku,nyku,nofe,nofn)
!
!
ccccccccccccccccccccccccccc
c 逆行列から9成分を求める c
ccccccccccccccccccccccccccc
c
cccccccccc材料定数ccccccccccccc
c
c　ヤング率：exx,eyy,ezz
c      write(*,*)'x軸方向のヤング率を入れて下さい'
c      read(*,*) exx
       exx=5.d2
c      write(*,*)'y軸方向のヤング率を入れて下さい'
c      read(*,*) eyy
       eyy=1.d3
c      write(*,*)'z軸方向のヤング率を入れて下さい'
c      read(*,*) ezz
       ezz=1.d4
c
c　ポアソン比：poixy,poixz,poiyz
       poixy=2.d-1
       poixz=2.5d-2
       poiyz=3.d-2
c
c　せん断弾性係数：gxy,gxz,gxy
       gxy=3.d1
       gxz=6.d2
       gyz=7.d2
c
c応力−ひずみ行列：e66
          do i=1,6
          do j=1,6
             e66(i,j)=0.d0
          end do
          end do
c  
       e66(1,1)=1.d0/exx
       e66(1,2)=-poixy/exx
       e66(1,3)=-poixz/exx
       e66(2,1)=e66(1,2)
       e66(2,2)=1.d0/eyy
       e66(2,3)=-poiyz/eyy
       e66(3,1)=e66(1,3)
       e66(3,2)=e66(2,3)
       e66(3,3)=1.d0/ezz
c
c
c
c
c 逆行列を取って、応力−ひずみ行列に変換する
       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
c
c
c      write(*,*)'xy面のせん断弾性係数を入れて下さい'
c      read(*,*) gxy
      e66(4,4)=gyz
c      write(*,*)'xz面のせん断弾性係数を入れて下さい'
c      read(*,*) gxz
      e66(5,5)=gxz
c      write(*,*)'yz面のせん断弾性係数を入れて下さい'
c      read(*,*) gyz
      e66(6,6)=gxy
c
c
! ******* 確認のため インプットデッキには書き込まれません ******
!             write(*,*) '応力-ひずみ行列：'
!            do i=1,6
!             write(*,211) (e66(i,j),j=1,6)
!            end do
!211         format(1p6d11.3)
!
cccccccccccccc
c固定端の座標c
cccccccccccccc
       write(7,210)
       do i=1,(nx+1)*(ny+1)
       write(7,*)i,','
       end do
210    format('*NSET,NSET=FIX')
c
cccccccccccccccccccccccccccccccc
cその他境界条件に必要な要素番号c
cccccccccccccccccccccccccccccccc
cccccccc固定端の中立軸の座標cccccccccc
       write(7,520)
       do i =1,nx+1
       write(7,*)ny/2*(nx+1)+i,','
       end do
520    format('*NSET,NSET=CHURITU')
c
cccccccc図心のあるy軸の節点番号cccccccccc
      write(7,'("*NSET,NSET=ZUSIN")')
      do i =1,ny+1 
      write(7,*)(i-1)*(nx+1)+(nx/2)+1,','
      end do
cccccccccc
c境界条件c
cccccccccc
       write(7,220)
       write(7,510)
220    format('*BOUNDARY')
510    format('FIX,3')
c
cccccc固定端の中立軸のyを拘束ccccccc
       write(7,521)
521    format('CHURITU,2')
c
cccccc半分に切断した面のxを拘束ccccccc
c       write(7,531)
c531    format('SETUDAN,1')
c
cccccc図心軸のxを拘束ccccccccccccc
       write(7,'("ZUSIN,1")')
c
!#############
! 材料の設定 #
!#############
!
!#########　鋼材　##########
       write(7,'("*MATERIAL,NAME=STE")') 
!デフォールトは等方性異方性ならTYPE=ANISOをつける
       write(7,'("*ELASTIC")')
!ヤング率、ポワソン比の順
       write(7,'("206000.,.3333")')
!######### 塑性を入れる場合(鋼材) #######
       write(7,'("*PLASTIC")')
       write(7,'("235.,0.")')
       write(7,'("235.,0.5")')
!
!######　直行異方性の場合(集成材)  #######
!
       write(7,231)
231    format('*MATERIAL,NAME=GLU')
       write(7,244)
244    format('*ELASTIC,TYPE=ORTHOTROPIC')
c240    format('*ELASTIC')
       write(7,243) e66(1,1),',',e66(1,2),',',e66(2,2),',',e66(1,3),','
     & ,e66(2,3),',',e66(3,3),',',e66(6,6),',',e66(5,5),','
       write(7,242) e66(4,4)
243    format(f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a,f12.4,a
     & ,f12.4,a)
242    format(f12.4,',','293.')
c########## 塑性を入れる場合(集成材) #########
!Experimental and analytical study on steel deck-gluam beam
!hybrid bridge behavior : Dr.KISS Lajos
!より集成材の降伏応力を引用
       write(7,'("*PLASTIC")')
       write(7,'("21.696,0.")')
       write(7,'("21.696,0.5")')
c
c
c
!#######################
! 要素に材料特性の付加 #
!#######################
! ここでどの要素がどんな材料か決定
      write(7,'("*SOLID SECTION,ELSET=glulam,MATERIAL=GLU")')
      write(7,'("*SOLID SECTION,ELSET=ukou,MATERIAL=STE")')
      write(7,'("*SOLID SECTION,ELSET=lkou,MATERIAL=STE")')
!
!
       write(7,'("*STEP")')
       write(7,'("*STATIC,SOLVER=SPOOLES")')
       write(7,'("*CLOAD")')
! 載荷
       call menseki(nx,ny,nz,nxt1,nxt2,zlt1,zlt2,nxkl,nxku,nykl,nyku,
     & zl,zm,ziku,zlkl,zlku,zmkl,zmku,nofnp,pakajyu,area,zakajyu,
     & nxt3,zlt3)
!
       write(7,'("*NODE PRINT,NSET=Nall")')
       write(7,'("U")')
       write(7,'("*EL PRINT,ELSET=glulam")')
       write(7,'("S")')
       write(7,'("*EL PRINT,ELSET=ukou")')
       write(7,'("S")')
       write(7,'("*EL PRINT,ELSET=lkou")')
       write(7,'("S")')
       write(7,'("*NODE FILE")')
       write(7,'("U")')
       write(7,'("*EL FILE,POSITION=AVERAGED AT NODES")')
       write(7,'("S")')
       write(7,'("*END STEP")')
!
      close(7)
      end
!
!###################### ここまでがメイン ###################
!
!
!
!
!
!
!
!
!######## 座標を出力するサブルーチン######
      subroutine zahyou(nx,ny,nz,zl,zm,ziku,nxt1,nxt2,nxt3,nxkl,nxku,
     &zlku,zlkl,zlt1,zlt2,zlt3,nyku,nykl,zmku,zmkl,nnofn,zax,zay,zaz)
      implicit real*8(a-h,o-z)
!
      zxkiten=-1.d0*zl/2.d0
      zykiten=1.d0*zm/2.d0
      write(7,'("*NODE,NSET=Nall")')
      nnofn=0
      do k=1,nz+1
      do j=1,ny+1
      do i=1,nx+1
       nnofn=nnofn+1
! xの座標出力
       if(i.le.nxt1+1)then
        zax=zxkiten+(i-1)*zlt3/nxt3
       else if(i.le.nxt3+nxkl+1)then
        zax=zxkiten+zlt3+(i-nxt3-1)*zlkl/nxkl
       else if(i.le.nxt1+nxkl+nxt2+1)then
        zax=zxkiten+zlt3+zlkl+(i-nxt1-nxkl-1)*zlt2/nxt2
       else if(i.le.(nx-nxkl)/2+1)then
        zax=zxkiten+(zl-zlku)/2.d0+(i-nxt3-nxkl-nxt2-1)*zlt1/nxt1
       else if(i.le.(nx+nxkl)/2+1)then
        zax=zxkiten+(zl-zlkl)/2.d0+(i-(nx-nxkl)/2-1)*zlkl/nxkl
       else if(i.le.(nx+nxkl)/2+nxt1+1)then
        zax=zxkiten+(zl+zlkl)/2.d0+(i-(nx+nxkl)/2-1)*zlt1/nxt1
       else if(i.le.nx-nxt3-nxkl+1)then
        zax=zxkiten+(zl+zlku)/2.d0+(i-(nx+nxkl)/2-nxt1-1)*zlt2/nxt2
       else if(i.le.nx-nxt3+1)then
        zax=zxkiten+zl-zlkl-zlt3+(i-nx+nxt3+nxkl-1)*zlkl/nxkl
       else
        zax=zxkiten+zl-zlt3+(i-nx+nxt3-1)*zlt3/nxt3
       end if
! yの座標出力
       if(j.le.nyku+1) then
        zay=zykiten+(j-1)*zmku/nyku
       else if(j.le.ny-nykl)then
        zay=zykiten+zmku+(j-nyku-1)*(zm-zmku-zmkl)/(ny-nyku-nykl)
       else 
        zay=zykiten+zm-zmkl+(j-ny+nykl-1)*zmkl/nykl
       end if
! zの座標出力
        zaz=(k-1)*ziku/nz
!
      write(7,'(i5,a,F10.5,a,F10.5,a,F10.5)')nnofn,',',zax,',',zay,',',
     &zaz
      end do
      end do
      end do
!
      return
      end
!
!######## 要素を出力するサブルーチン ########
      subroutine youso(nx,ny,nz,nxt1,nxt2,nxt3,nxkl,nykl,nxku,nyku,
     &nofe,nofn)
      implicit real*8(a-h,o-z)
      dimension nofn(8)
!cccccccccccc 鋼材と木材の節点番号出力 cccccccccccccccc
      do k=1,nz
      do j=1,ny
      do i=1,nx
!****** 上部鋼出力 ******
      if((i.ge.nxt3+nxkl+nxt2+1).and.(i.le.nxt3+nxkl+nxt2+nxku)
     &.and.(j.le.nyku)) then
        nofe=(k-1)*nx*ny+nx*(j-1)+i
        nofn(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
        nofn(2)=nofn(1)+1
        nofn(3)=nofn(1)+(nx+1)
        nofn(4)=nofn(3)+1
        nofn(5)=nofn(1)+(nx+1)*(ny+1)
        nofn(6)=nofn(2)+(nx+1)*(ny+1)
        nofn(7)=nofn(3)+(nx+1)*(ny+1)
        nofn(8)=nofn(4)+(nx+1)*(ny+1)
        write(7,'("*ELEMENT,TYPE=C3D8,ELSET=ukou")')
        write(7,*)nofe,',',nofn(1),',',nofn(3),',',nofn(7),
     & ',',nofn(5),',',nofn(2),',',nofn(4),',',nofn(8),',',nofn(6)
!
!****** 下部鋼出力 ******
      else if((i.ge.nxt3+1).and.(i.le.nxt3+nxkl)
     & .and.(j.ge.ny-nykl+1).and.(j.le.ny))then
        nofe=(k-1)*nx*ny+nx*(j-1)+i
        nofn(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
        nofn(2)=nofn(1)+1
        nofn(3)=nofn(1)+(nx+1)
        nofn(4)=nofn(3)+1
        nofn(5)=nofn(1)+(nx+1)*(ny+1)
        nofn(6)=nofn(2)+(nx+1)*(ny+1)
        nofn(7)=nofn(3)+(nx+1)*(ny+1)
        nofn(8)=nofn(4)+(nx+1)*(ny+1)
        write(7,'("*ELEMENT,TYPE=C3D8,ELSET=lkou")')
        write(7,*)nofe,',',nofn(1),',',nofn(3),',',nofn(7),
     & ',',nofn(5),',',nofn(2),',',nofn(4),',',nofn(8),',',nofn(6)
!
      else if((i.ge.(nx-nxkl)/2+1).and.(i.le.(nx+nxkl)/2)
     & .and.(j.ge.ny-nykl+1).and.(j.le.ny))then
        nofe=(k-1)*nx*ny+nx*(j-1)+i
        nofn(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
        nofn(2)=nofn(1)+1
        nofn(3)=nofn(1)+(nx+1)
        nofn(4)=nofn(3)+1
        nofn(5)=nofn(1)+(nx+1)*(ny+1)
        nofn(6)=nofn(2)+(nx+1)*(ny+1)
        nofn(7)=nofn(3)+(nx+1)*(ny+1)
        nofn(8)=nofn(4)+(nx+1)*(ny+1)
        write(7,'("*ELEMENT,TYPE=C3D8,ELSET=lkou")')
        write(7,*)nofe,',',nofn(1),',',nofn(3),',',nofn(7),
     & ',',nofn(5),',',nofn(2),',',nofn(4),',',nofn(8),',',nofn(6)
!
      else if((i.ge.nx-nxkl-nxt3+1).and.(i.le.nx-nxt3)
     & .and.(j.ge.ny-nykl+1).and.(j.le.ny))then
        nofe=(k-1)*nx*ny+nx*(j-1)+i
        nofn(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
        nofn(2)=nofn(1)+1
        nofn(3)=nofn(1)+(nx+1)
        nofn(4)=nofn(3)+1
        nofn(5)=nofn(1)+(nx+1)*(ny+1)
        nofn(6)=nofn(2)+(nx+1)*(ny+1)
        nofn(7)=nofn(3)+(nx+1)*(ny+1)
        nofn(8)=nofn(4)+(nx+1)*(ny+1)
        write(7,'("*ELEMENT,TYPE=C3D8,ELSET=lkou")')
        write(7,*)nofe,',',nofn(1),',',nofn(3),',',nofn(7),
     & ',',nofn(5),',',nofn(2),',',nofn(4),',',nofn(8),',',nofn(6)
!
!****** 集成材出力 ******
       else
        nofe=(k-1)*nx*ny+nx*(j-1)+i
        nofn(1)=(nx+1)*(ny+1)*(k-1)+(nx+1)*(j-1)+i
        nofn(2)=nofn(1)+1
        nofn(3)=nofn(1)+(nx+1)
        nofn(4)=nofn(3)+1
        nofn(5)=nofn(1)+(nx+1)*(ny+1)
        nofn(6)=nofn(2)+(nx+1)*(ny+1)
        nofn(7)=nofn(3)+(nx+1)*(ny+1)
        nofn(8)=nofn(4)+(nx+1)*(ny+1)
        write(7,'("*ELEMENT,TYPE=C3D8,ELSET=glulam")')
        write(7,*)nofe,',',nofn(1),',',nofn(3),',',nofn(7),
     & ',',nofn(5),',',nofn(2),',',nofn(4),',',nofn(8),',',nofn(6)
       end if
       end do
       end do
       end do
!
      return
      end
!
!######## 各要素の面積に応じて載荷するサブルーチン ########
      subroutine menseki(nx,ny,nz,nxt1,nxt2,zlt1,zlt2,nxkl,nxku,nykl,
     & nyku,zl,zm,ziku,zlkl,zlku,zmkl,zmku,nofnp,pakajyu,area,zakajyu,
     & nxt3,zlt3)
      implicit real*8(a-h,o-z)
      dimension area(100,100),yoko(100),tate(100),ax(100),ay(100)
!******* 各要素の面積を計算 ******
! xの座標出力
      do i=2,nx+1
       if(i.le.nxt1+1)then
        yoko(i)=(i-1)*zlt3/nxt3
       else if(i.le.nxt3+nxkl+1)then
        yoko(i)=zlt3+(i-nxt3-1)*zlkl/nxkl
       else if(i.le.nxt1+nxkl+nxt2+1)then
        yoko(i)=zlt3+zlkl+(i-nxt1-nxkl-1)*zlt2/nxt2
       else if(i.le.(nx-nxkl)/2+1)then
        yoko(i)=(zl-zlku)/2.d0+(i-nxt3-nxkl-nxt2-1)*zlt1/nxt1
       else if(i.le.(nx+nxkl)/2+1)then
        yoko(i)=(zl-zlkl)/2.d0+(i-(nx-nxkl)/2-1)*zlkl/nxkl
       else if(i.le.(nx+nxkl)/2+nxt1+1)then
        yoko(i)=(zl+zlkl)/2.d0+(i-(nx+nxkl)/2-1)*zlt1/nxt1
       else if(i.le.nx-nxt3-nxkl+1)then
        yoko(i)=(zl+zlku)/2.d0+(i-(nx+nxkl)/2-nxt1-1)*zlt2/nxt2
       else if(i.le.nx-nxt3+1)then
        yoko(i)=zl-zlkl-zlt3+(i-nx+nxt3+nxkl-1)*zlkl/nxkl
       else
        yoko(i)=zl-zlt3+(i-nx+nxt3-1)*zlt3/nxt3
       end if
      end do
! yの座標出力
      do j=2,ny+1
       if(j.le.nyku+1) then
        tate(j)=(j-1)*zmku/nyku
       else if(j.le.ny-nykl)then
        tate(j)=zmku+(j-nyku-1)*(zm-zmku-zmkl)/(ny-nyku-nykl)
       else 
        tate(j)=zm-zmkl+(j-ny+nykl-1)*zmkl/nykl
       end if
      end do
! 各要素の横の長さ
       yoko(1)=0
       do i=1,nx
       ax(i)=yoko(i+1)-yoko(i)
       end do
! 各要素の縦の長さ
       tate(1)=0
       do j=1,ny
       ay(j)=tate(j+1)-tate(j)
       end do
! 各要素の面積
      do j=1,ny
      do i=1,nx
      area(j,i)=ay(j)*ax(i)
      end do
      end do
!****** 自由端の各節点に載荷 ******
      write(*,'("載荷する荷重を入力して下さい")')
      read*,zkajyu
!       zkajyu=100.d0
       zakajyu=zkajyu/zl/zm*-1.d0
!++++++ 4要素と接する節点 ++++++
! 
      do j=1,ny
      do i=1,nx
       pakajyu=(area(j-1,i-1)+area(j-1,i)+area(j,i-1)+area(j,i))/4.d0
     &*zakajyu
       if((i.gt.1).and.(i.lt.nx+1).and.(j.gt.1).and.(j.lt.ny+1))then
       nofnp=(nx+1)*(ny+1)*nz+(j-1)*(nx+1)+i
       write(7,*)nofnp,',','2',',',pakajyu
      end if
      end do
      end do
!
!++++++ 2要素と接する節点 ++++++
! 上
      do i =1,nx
       pakajyu=(area(1,i-1)+area(1,i))/4.d0*zakajyu
       if(i.gt.1)then
       nofnp=(nx+1)*(ny+1)*nz+i
       write(7,*)nofnp,',','2',',',pakajyu
       end if
      end do
! 下
      do i =1,nx
       pakajyu=(area(ny,i-1)+area(ny,i))/4.d0*zakajyu
       if(i.gt.1)then
       nofnp=(nx+1)*(ny+1)*nz+ny*(nx+1)+i
       write(7,*)nofnp,',','2',',',pakajyu
       end if
      end do
! 右
      do j =1,ny
       pakajyu=(area(j-1,nx)+area(j,nx))/4.d0*zakajyu
       if(j.gt.1)then
       nofnp=(nx+1)*(ny+1)*nz+(j-1)*(nx+1)+1
      write(7,*)nofnp,',','2',',',pakajyu
      end if
      end do
! 左
      do j =1,ny
       pakajyu=(area(j-1,1)+area(j,1))/4.d0*zakajyu
       if(j.gt.1)then
       nofnp=(nx+1)*(ny+1)*nz+(j-1)*(nx+1)+nx+1
      write(7,*)nofnp,',','2',',',pakajyu
       end if
      end do
!++++++ どことも接しない節点 ++++++
!
      write(7,*)(nx+1)*(ny+1)*nz+1,',','2',',',area(1,1)/4.d0*zakajyu
      write(7,*)(nx+1)*(ny+1)*nz+nx+1,',','2',',',
     &area(1,nx)/4.d0*zakajyu
      write(7,*)(nx+1)*(ny+1)*nz+ny*(nx+1)+1,',','2',',',
     &area(ny,1)/4.d0*zakajyu
      write(7,*)(nx+1)*(ny+1)*nz+ny*(nx+1)+nx+1,',','2',',',
     &area(ny,nx)/4.d0*zakajyu
!
      return
      end
! 
