module commvar
      real*8 :: sfr(3,3),sfo(3,3),sfro(3,3),sfrr(3,3),sfoo(3,3)
      real*8 :: sfoo1(3,3),sfoo4(3,3),sfoo5(3,3),sfoo7(3,3)
      real*8 :: bfo(3,3,3),bfro(3,3,3),bfrr(3,3,3),bfoo(3,3,3,3)
      real*8 :: bfoo1(3,3,3,3),bfoo4(3,3,3,3,3),bfoo5(3,3,3,3,3),bfoo7(3,3,3,3,3)

      contains

      subroutine values
      implicit none
      integer :: i,j
      open(unit=101,file='sidat',action='read')
      read(101,*) sfr(1,1),sfo(1,1),sfrr(1,1),sfro(1,1),sfoo(1,1),sfoo1(1,1),&
                  sfoo4(1,1),sfoo5(1,1),sfoo7(1,1)
      close(101)
      open(unit=101,file='gedat',action='read')
      read(101,*) sfr(2,2),sfo(2,2),sfrr(2,2),sfro(2,2),sfoo(2,2),sfoo1(2,2),&
                  sfoo4(2,2),sfoo5(2,2),sfoo7(2,2)
      close(101)
      open(unit=101,file='cdat',action='read')
      read(101,*) sfr(3,3),sfo(3,3),sfrr(3,3),sfro(3,3),sfoo(3,3),sfoo1(3,3),&
                  sfoo4(3,3),sfoo5(3,3),sfoo7(3,3)
      close(101)
      open(unit=101,file='sicdat',action='read')
      read(101,*) sfr(1,3),sfo(1,3),sfrr(1,3),sfro(1,3),sfoo(1,3),sfoo1(1,3),&
                  sfoo4(1,3),sfoo5(1,3),sfoo7(1,3)
      close(101)
      open(unit=101,file='sicdat',action='read')
      read(101,*) sfr(3,1),sfo(3,1),sfrr(3,1),sfro(3,1),sfoo(3,1),sfoo1(3,1),&
                  sfoo4(3,1),sfoo5(3,1),sfoo7(3,1)
      close(101)
      do i=1,3
        do j=1,3
          if (i==1.and.j==3) goto 10
          if (i==3.and.j==1) goto 10
          if (i==j) goto 10
          sfr(i,j)=(sfr(i,i)+sfr(j,j))/dble(2)
          sfo(i,j)=(sfo(i,i)+sfo(j,j))/dble(2)
          sfrr(i,j)=(sfrr(i,i)+sfrr(j,j))/dble(2)
          sfro(i,j)=(sfro(i,i)+sfro(j,j))/dble(2)
          sfoo(i,j)=(sfoo(i,i)+sfoo(j,j))/dble(2)
          sfoo1(i,j)=(sfoo1(i,i)+sfoo1(j,j))/dble(2)
          sfoo4(i,j)=(sfoo4(i,i)+sfoo4(j,j))/dble(2)
          sfoo5(i,j)=(sfoo5(i,i)+sfoo5(j,j))/dble(2)
          sfoo7(i,j)=(sfoo7(i,i)+sfoo7(j,j))/dble(2)
10        continue
        end do
      end do
      end subroutine values

      end module commvar

!	MAIN PROGRAM

      program phd4
      use commvar
      implicit none
      integer:: i,j,k,l,m           
      call values
      do i=1,3
        do j=1,3
          do k=1,3
            bfo(i,j,k)=(sfo(i,j)+sfo(j,k))/dble(2)
            bfrr(i,j,k)=(sfrr(i,j)+sfrr(j,k))/dble(2)
            bfro(i,j,k)=(sfro(i,j)+sfro(j,k))/dble(2)
            do l=1,3
              bfoo(i,j,k,l)=(sfoo(i,j)+sfoo(i,k)+sfoo(i,l))/dble(3)
              bfoo1(i,j,k,l)=(sfoo1(i,j)+sfoo1(j,k)+sfoo1(k,l))/dble(3)
              do m=1,3
                bfoo4(i,j,k,l,m)=(sfoo4(i,j)+sfoo4(j,k)+sfoo4(k,l)+&
                                  sfoo4(l,m))/dble(4)
                bfoo5(i,j,k,l,m)=(sfoo5(i,j)+sfoo5(j,k)+sfoo5(k,l)+&
                                  sfoo5(l,m))/dble(4) 
                bfoo7(i,j,k,l,m)=(sfoo7(i,l)+sfoo7(j,l)+sfoo7(k,l)+&
                                  sfoo7(k,m))/dble(4)
              end do
            end do
          end do
        end do
      end do
      call writeout(1) 
      end program phd4

      subroutine writeout(a)
      use commvar
      implicit none
      integer :: a,i,j,k,l,m
      if (a==1)open(unit=22,file='begin.f90',status='replace')
      if (a==2)open(unit=22,file='param.f90',status='replace')
      if (a==3)open(unit=22,file='minim.f90',status='replace')
      write(22,*)'      module para'
      write(22,*)'      real*8::bfr(0:3,0:3),bfo(0:3,0:3,0:3)'
      write(22,*)'      real*8::bfro(0:3,0:3,0:3),bfrr(0:3,0:3,0:3)'
      write(22,*)'      real*8::bfoo(0:3,0:3,0:3,0:3)'
    write(22,*)'      real*8::bfoo1(0:3,0:3,0:3,0:3),bfoo4(0:3,0:3,0:3,0:3,0:3)'
    write(22,*)'      real*8::bfoo5(0:3,0:3,0:3,0:3,0:3),bfoo7(0:3,0:3,0:3,0:3,0:3)'
      write(22,*)'      contains'
      write(22,*)'      subroutine parameters'
1     format(a10,i1,a1,i1,a2,f13.8)
      do i=0,3
        do j=0,3
      if(i*j/=0)write(22,1)'      bfr(',i,',',j,')=',sfr(i,j)  
      if(i*j==0) write(22,1)'      bfr(',i,',',j,')=',dble(0)  
        end do
      end do
2     format(a10,i1,a1,i1,a1,i1,a2,f10.8)
      do i=0,3
         do j=0,3
           do k=0,3
      if(i*j*k/=0)write(22,2)'      bfo(',i,',',j,',',k,')=',bfo(i,j,k)
      if(i*j*k==0)write(22,2)'      bfo(',i,',',j,',',k,')=',dble(0)
           end do
         end do
      end do
3     format(a11,i1,a1,i1,a1,i1,a2,f10.8)
      do i=0,3
         do j=0,3
           do k=0,3
      if(i*j*k/=0)write(22,3)'      bfro(',i,',',j,',',k,')=',bfro(i,j,k)
      if(i*j*k==0)write(22,3)'      bfro(',i,',',j,',',k,')=',dble(0)
           end do
         end do
      end do
      do i=0,3
         do j=0,3
           do k=0,3
      if(i*j*k/=0)write(22,3)'      bfrr(',i,',',j,',',k,')=',bfrr(i,j,k)
      if(i*j*k==0)write(22,3)'      bfrr(',i,',',j,',',k,')=',dble(0)
           end do
         end do
      end do
4     format(a11,i1,a1,i1,a1,i1,a1,i1,a2,f10.8)
      do i=0,3
         do j=0,3
           do k=0,3
             do l=0,3
      if(i*j*k*l/=0)write(22,4)'      bfoo(',i,',',j,',',k,',',l,')=',&
                                      bfoo(i,j,k,l)
      if(i*j*k*l==0)write(22,4)'      bfoo(',i,',',j,',',k,',',l,')=',&
                                      dble(0)
             end do 
           end do
         end do
      end do
5     format(a12,i1,a1,i1,a1,i1,a1,i1,a2,f10.8)
      do i=0,3
         do j=0,3
           do k=0,3
             do l=0,3
      if(i*j*k*l/=0)write(22,5)'      bfoo1(',i,',',j,',',k,',',l,')=',&
                                      bfoo1(i,j,k,l)
      if(i*j*k*l==0)write(22,5)'      bfoo1(',i,',',j,',',k,',',l,')=',&
                                      dble(0)
             end do
           end do
         end do
      end do
6     format(a12,i1,a1,i1,a1,i1,a1,i1,a1,i1,a2,f10.8)
      do i=0,3
         do j=0,3
           do k=0,3
             do l=0,3
               do m=0,3
      if(i*j*k*l*m/=0)&
      write(22,6)'      bfoo4(',i,',',j,',',k,',',l,',',m,')=',bfoo4(i,j,k,l,m)
      if(i*j*k*l*m==0)&
      write(22,6)'      bfoo4(',i,',',j,',',k,',',l,',',m,')=',dble(0)
               end do
             end do 
           end do
         end do
      end do
      do i=0,3
         do j=0,3
           do k=0,3
             do l=0,3
               do m=0,3
      if(i*j*k*l*m/=0)&
      write(22,6)'      bfoo5(',i,',',j,',',k,',',l,',',m,')=',bfoo5(i,j,k,l,m)
      if(i*j*k*l*m==0)&
      write(22,6)'      bfoo5(',i,',',j,',',k,',',l,',',m,')=',dble(0)
               end do
             end do 
           end do
         end do
      end do
      do i=0,3
         do j=0,3
           do k=0,3
             do l=0,3
               do m=0,3
      if(i*j*k*l*m/=0)&
      write(22,6)'      bfoo7(',i,',',j,',',k,',',l,',',m,')=',bfoo7(i,j,k,l,m)
      if(i*j*k*l*m==0)&
      write(22,6)'      bfoo7(',i,',',j,',',k,',',l,',',m,')=',dble(0)
               end do
             end do 
           end do
         end do
      end do
      write(22,*)'      end subroutine parameters'
      write(22,*)'      end module para'
      close(22)
      return
      end subroutine writeout