      program zgousei2 verb2
      implicit real*8(a-h,o-z)
      dimension nofn(8),b(3,3),e66(6,6),area(100,100)
!
!######## 各断面の諸定数入力 ########
!
!      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(*,*)zlkl
      zlkl=4.5d-3
!
!      print*,'下部に挿入する鋼板のy軸の長さを入れて下さい'
!      read(*,*)zmkl
      zmkl=4.d-2
!
!      print*,'上部に挿入する鋼板のx軸の長さを入れてください'
!      read(*,*)zlku
       zlku=9.d-3
!
!       print*,'上部に挿入する鋼板のy軸の長さを入れてください'
!       read(*,*)zmku
       zmku=4.d-2
!
!       print*,'鋼板は図心軸からどのくらいのところに入れますか'
!       print*,'鋼板の厚みは考慮せず入れる位置だけ指定してください'
!       read(*,*) zlt
       zlt=1.d-2
! 上部鋼板と下部鋼板の間のx軸方向の長さ
       zlt1=zlt-(zlku+zlkl)/2.d0
! 端から下部鋼板までの長さ
       zlt2=(zl-zlku-(zlkl+zlt1)*2.d0)/2.d0
!
!
!######## 要素数入力 ########
!
!       print*,'x軸方向の全要素数(偶数)を入れて下さい(鋼＋木)'
!       read(*,*) nx
!       nx=nx/2
        nx=8
!
!       print*,'y軸方向の全要素数(偶数)を入れて下さい(鋼＋木)'
!       read(*,*) ny
        ny=10
!
!       print*,'z軸方向の全要素数を入れて下さい(鋼＋木)'
!       read(*,*) nz
        nz=100
!
!       print*,'上部の鋼板だけのx軸方向の要素数(偶数)を入れてください'
!       read(*,*) nxku
       nxku=2
!       print*,'上部の鋼板だけのy軸方向の要素数を入れてください'
!       read(*,*) nyku
       nyku=4
!
!       print*,'下部の鋼板だけのx軸方向の要素数を入れてください'
!       read(*,*) nxkl
       nxkl=1
!
!       print*,'下部の鋼板だけのy軸方向の要素数を入れてください'
!       read(*,*) nxkl
       nykl=4
!
!        print*,'上部と下部の鋼板の間の木材のx方向の要素数を入れて下さい'
!        read(*,*) nxt1
        nxt1=1
!
! 端から下部鋼板までの集成材の要素数
        nxt2=(nx-nxku-(nxkl+nxt1)*2)/2
!
!
!
!######## 出力 #########
      open(7,file='zgousei2.inp')
!
      call zahyou(zl,zm,ziku,zlku,zlkl,zmku,zmkl,zlt,zlt1,zlt2,
     &nx,ny,nz,nxkl,nxku,nykl,nyku,nxt1,nxt2,nnofn,zax,zay,zaz)
!
      call youso(nx,ny,nz,nxt1,nxt2,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")')
!cccccccccccccc
!c自由端の座標c
!cccccccccccccc
!       write(7,270)
!       nnofe=(nx+1)*(ny+1)*(nz+1)-(nx+1)*(ny+1)
!       do i=1,(nx+1)*(ny+1)
!       nnofe=nnofe+1
!       write(7,*)nnofe,','
!       end do
!270    format('*NSET,NSET=LAST')
c
c
c
c
       write(7,'("*STEP")')
       write(7,'("*STATIC,SOLVER=SPOOLES")')
       write(7,'("*CLOAD")')
! 載荷
      print*,'どのように載荷しますか?'
      print*,'各節点を構成する要素数に応じて載荷:0を入れて下さい'
      print*,'各節点を構成する要素の各面積に応じて載荷:1を入れて下さい'
      read*,nsentaku
      if(nsentaku.eq.0)then
       call ekazu(nx,ny,nz,pkajyu,p2kajyu,p4kajyu)
      else if(nsentaku.eq.1)then
       call menseki(nx,ny,nz,nxt1,nxt2,zlt1,zlt2,nxkl,nxku,nykl,
     & nyku,zl,zm,ziku,zlkl,zlku,zmkl,zmku,nofnp,pakajyu,area,zakajyu)
      end if
!
       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(zl,zm,ziku,zlku,zlkl,zmku,zmkl,zlt,zlt1,zlt2,
     &nx,ny,nz,nxkl,nxku,nykl,nyku,nxt1,nxt2,nnofn,zax,zay,zaz)
      implicit real*8(a-h,o-z)
      zxkiten=-1.d0*zl/2.d0
      zykiten=1.d0*zm/2.d0
! xの座標出力
      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
       if(i.le.nxt2+1) then
        zax=zxkiten+(i-1)*zlt2/nxt2
       else if(i.le.nxt2+nxkl+1)then
        zax=zxkiten+zlt2+(i-nxt2-1)*zlkl/nxkl
       else if(i.le.nxt2+nxkl+nxt1+1)then
        zax=zxkiten+zlt2+zlkl+(i-nxt2-nxkl-1)*zlt1/nxt1
       else if(i.le.(nx+nxku)/2+1)then
        zax=zxkiten+zlt2+zlkl+zlt1+(i-nxt2-nxkl-nxt1-1)*zlku/nxku
       else if(i.le.(nx+nxku)/2+nxt1+1)then
        zax=zxkiten+(zl+zlku)/2.d0+(i-(nx+nxku)/2-1)*zlt1/nxt1
       else if(i.le.nx+1-nxt2)then
        zax=zxkiten+(zl+zlku)/2.d0+zlt1+(i-(nx+nxku)/2-nxt1-1)*zlkl/nxkl
       else
        zax=zxkiten+(zl-zlt2)+(i-nx-1+nxt2)*zlt2/nxt2
       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,F7.5,a,F7.5,a,F7.5)')nnofn,',',zax,',',zay,',',zaz
       end do
       end do
       end do
      return
      end
!
!######## 要素を出力するサブルーチン ########
      subroutine youso(nx,ny,nz,nxt1,nxt2,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.nxt2+nxkl+nxt1+1).and.(i.le.nxt2+nxkl+nxt1+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.nxt2+1).and.(i.le.nxt2+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-nxt2+1).and.(i.le.nx-nxt2)
     & .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 ekazu(nx,ny,nz,pkajyu,p2kajyu,p4kajyu)
      implicit real*8(a-h,o-z)
! p,2p,4pに分けて載荷するけど、そのpの求め方(単位はMN)
       print*,'載荷する荷重を入力して下さい'
       read(*,*)zkajyu
       pkajyu=zkajyu/4.d0/nx/ny
!
!****** p載荷する節点番号と荷重の方向と大きさの出力 ******
       do j =1,ny+1,ny
       do i =1,nx+1,nx
       write(7,*) (nx+1)*(ny+1)*nz+(j-1)*(nx+1)+i,',','2',',',pkajyu
       end do
       end do
!
!****** 2p載荷する節点番号と荷重の方向と大きさの出力 ******
       p2kajyu=pkajyu*2.d0
       do j =1,ny+1,ny
       do i =2,nx
       write(7,*) (nx+1)*(ny+1)*nz+(j-1)*(nx+1)+i,',','2',',',p2kajyu
       end do
       end do
!
       do j =2,ny
       do i =1,nx+1,nx
       write(7,*) (nx+1)*(ny+1)*nz+(j-1)*(nx+1)+i,',','2',',',p2kajyu
       end do
       end do
!
!****** 4p載荷する節点番号と荷重の方向と大きさの出力 ******
       p4kajyu=pkajyu*4.d0
       do j =2,ny
       do i =2,nx
       write(7,*) (nx+1)*(ny+1)*nz+(j-1)*(nx+1)+i,',','2',',',p4kajyu
       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)
      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
      do j=2,ny+1
       if(i.le.nxt2+1) then
        yoko(i)=(i-1)*zlt2/nxt2
       else if(i.le.nxt2+nxkl+1)then
        yoko(i)=zlt2+(i-nxt2-1)*zlkl/nxkl
       else if(i.le.nxt2+nxkl+nxt1+1)then
        yoko(i)=zlt2+zlkl+(i-nxt2-nxkl-1)*zlt1/nxt1
       else if(i.le.(nx+nxku)/2+1)then
        yoko(i)=zlt2+zlkl+zlt1+(i-nxt2-nxkl-nxt1-1)*zlku/nxku
       else if(i.le.(nx+nxku)/2+nxt1+1)then
        yoko(i)=(zl+zlku)/2.d0+(i-(nx+nxku)/2-1)*zlt1/nxt1
       else if(i.le.nx+1-nxt2)then
        yoko(i)=(zl+zlku)/2.d0+zlt1+(i-(nx+nxku)/2-nxt1-1)*zlkl/nxkl
       else
        yoko(i)=(zl-zlt2)+(i-nx-1+nxt2)*zlt2/nxt2
       end if
! yの座標出力
       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-zmku+(j-ny+nykl-1)*zmkl/nykl
       end if
       end do
       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
! 
