program main

 implicit none

 integer              :: n,m,i,j,sp,ntmp,k,np,ierr
 real*8               :: sigma,mu,pi,cdiag,ires,r,iold,psca,x,g
 integer,allocatable  :: d(:),itmp(:)
 real*8,allocatable   :: a(:,:),b(:,:),c(:,:),tmp(:)


!1)
 write(*,*) 'Value of n'
 do
   read(*,*) n
   if (n.gt.1) exit
 enddo

 write(*,*) 'Value of m'
 do
   read(*,*) m 
   if (m.gt.1) exit
 enddo  

 write(*,*) 'Value of sigma'
 do
   read(*,*) sigma 
   if (sigma.gt.0.d0) exit
 enddo

 write(*,*) 'Value of mu'
 do
   read(*,*) mu 
   if (mu.gt.0.d0) exit
 enddo

!2)
 pi=dacos(-1.d0)

 allocate(a(n,n),b(m,m),c(n,n))

 do i=1,n
    do j=1,n
       x=dble(i+j)
       a(i,j)=g(x,sigma,mu,pi)
    enddo
 enddo   

 do i=1,m
    do j=1,m
       x=dble(i+j)
       b(i,j)=g(x,sigma,0.d0,pi)
    enddo
 enddo

!3)
 call define_c(a,b,c,n,m)

 cdiag=0.d0
 do i=1,n
    cdiag=cdiag+c(i,i)
 enddo   

 write(*,'("Sum_i C(i,i)", F12.5)') cdiag

 write(*,*) 'Value of R'
 do 
   read(*,*) r
   if (r.gt.0.d0) exit
 enddo

 np=10
 iold=1.d10
 if (cdiag.gt.0.5d0) then
    do
      call compute_int(ires,-r,r,np,sigma,mu,pi)
      if (abs(ires-iold).le.1.d-6) exit 
      np=np+5
      iold=ires
    enddo
 endif

!4)
 open(10,file='input.dat',status='old',iostat=ierr,err=200)
 allocate(tmp(n))
 do i=1,n
    read(10,*) tmp(i)
 enddo
 close(10)

 psca=0.d0
 do i=1,n
    psca=psca+tmp(i)*c(i,i)
 enddo

 allocate(d(n))
 do i=1,n
    d(i)=ceiling(b(i,i)+psca)-i 
 enddo

!5)
 if (mod(n,2).eq.0) then
    ntmp=n/2
 else
    ntmp=n/2+1
 endif   
 allocate(itmp(ntmp))

 k=0
 do i=1,n,2
    k=k+1
    itmp(k)=d(i)
 enddo

 call sort(itmp,ntmp)

 k=0
 do i=1,n,2
    k=k+1
    d(i)=itmp(k)
 enddo

!6)
 if (cdiag.gt.0.5d0) then
    open(10,file='integral.dat')
    write(10,'("Integral is", F12.5, " with points: ", I5)') ires,np 
    close(10)
 endif

 open(10,file='a.dat')
 do i=1,n
    do j=1,n
       write(10,100) i,j,a(i,j)
    enddo
 enddo
 close(10)

 open(10,file='b.dat')
 do i=1,m
    do j=1,m
       write(10,100) i,j,b(i,j)
    enddo   
 enddo
 close(10)

 open(10,file='c.dat')
 do i=1,n
    do j=1,n
       write(10,100) i,j,c(i,j)
    enddo
 enddo
 close(10)

 open(10,file='d.dat')
 do i=1,n
    write(10,110) i,d(i)
 enddo
 close(10)

100 format(I5,I5,F12.5)
110 format(I5,I5)

200 if (ierr.ne.0) then
       write(*,*) 'Error in opening the input file'
       stop
    endif

 deallocate(a,b,c)
 deallocate(d,tmp,itmp)

 stop

end program main      


real*8 function g(x,sigma,mu,pi)

 implicit none

 real*8, intent(in) :: x,sigma,mu,pi

 real*8             :: tmp

 tmp = dsqrt(2.d0*pi)*sigma
 tmp = 1.d0/tmp
 g = tmp*dexp(-(x-mu)**2/(2.d0*sigma**2))  

 return

end function g        


subroutine define_c(a,b,c,n,m)

 implicit none

 integer, intent(in)   :: n,m
 real*8,  intent(in)   :: a(n,n),b(m,m)
 real*8,  intent(out)  :: c(n,n) 

 integer               :: i,j
 real*8                :: csum

 csum=0.d0

 do i=n+1,m
    do j=1,n
       csum=csum+b(i,j)
    enddo   
 enddo   

 do i=n+1,m
    do j=n+1,m
       csum=csum+b(i,j)
    enddo
 enddo

 do i=1,n
    do j=n+1,m
       csum=csum+b(i,j)
    enddo
 enddo

 do i=1,n
    do j=1,n
       c(i,j)=a(i,j)+b(i,j)+csum
    enddo
 enddo   

 return

end subroutine define_c        

subroutine compute_int(res,a,b,n,sigma,mu,pi)

 implicit none

 integer, intent(in)  :: n
 real*8,  intent(in)  :: a,b,sigma,mu,pi
 real*8,  intent(out) :: res

 integer              :: i
 real*8               :: h,g,x

 h=(b-a)/dble(n)
 res=g(a,sigma,0.d0,pi)+g(b,sigma,0.d0,pi)
 res=0.5d0*res

 do i=1,n-1
    x = a+i*h
    res=res+g(x,sigma,0.d0,pi)
 enddo
 res=res*h

 return

end subroutine compute_int        

subroutine sort(a,n)

 implicit none

 integer, intent(in)    :: n
 integer, intent(inout) :: a(n)
 

 integer                :: i,j
 integer                :: tmp

 do i=1,n
    do j=2,n-(i-1)
       if (a(j-1).gt.a(j)) then
          tmp=a(j-1)
          a(j-1)=a(j)
          a(j)=tmp
       endif
    enddo
 enddo

 return

end subroutine sort