program magic
*
* Implementation of the De la Lubere's algorithm
* for odd magic squares n x n  (n = 2m+1)
*
      parameter(m=5)
      parameter (n=2*m+1)
      dimension k(n,n)
*
* Initialize matrix elements to zero
*
      do 10 i = 1,n
      do 10 j =1,n
 10   k(i,j)=0
*
* Number 1 is always set in the midle position of the first row
*
      k(1,m+1)= 1
* 
* ir and ic are the coordinates of the most recently filled square
* The coordinate system is illustrated in the following few comments
*
*     (1,1) (1,2) (1,3) ... (1,n)
*     (2,1) (2,2) (2,3) ... (2,n)
*      ...
*      ...
*     (n,1) (n,2) (n,3) ... (n,n)
*
      ir = 1
      ic = m+1
*
* Position one by one the rest n**2 - 1 numbers
*
      do 1 i= 2,n**2
* irn is the new candidate row number
      irn = ir -1
      if (irn .eq. 0) irn = n
* icn is the new candidate coloumn number
      icn = ic + 1
      if (icn .eq. n+1 ) icn = 1
* check if position is free
      if (k(irn,icn).eq.0) then
* it is free. Hence update the coordinates
         ir =irn
         ic = icn
      else
* it is filled. Prepare to fill the position below the old one
        ir = ir + 1
      endif
* fill the position
         k(ir,ic) = i
1     continue
*
* Done . Print the results
*
      do 2 i = 1,n
      write(*,'(19(i4,x))')(k(i,j),j=1,n)
 2    continue
      end