program main
  use sep
  implicit none

  integer :: nmx, nmy, nmz, nax, nay
  integer :: imx, imy, imz, iax, iay
  real    :: dmx, dmy, dmz, dax, day
  real    :: omx, omy, omz, oax, oay
  real, dimension(:,:,:,:), allocatable :: dati
  real, dimension(:,:)    , allocatable :: dato
  real                                  :: PI, wtx, wty, ax, ay

  PI = 3.14159265359

  call initpar()

  call from_history("n1", nax); call from_history("d1", dax); call from_history("o1", oax)
  call from_history("n2", nay); call from_history("d2", day); call from_history("o2", oay)
  call from_history("n3", nmx); call from_history("d3", dmx); call from_history("o3", omx)
  call from_history("n4", nmy); call from_history("d4", dmy); call from_history("o4", omy)
  call from_history("n5", nmz); call from_history("d5", dmz); call from_history("o5", omz)

  call to_history("n1", 1); call to_history("d1", 1.); call to_history("o1", 0.)
  call to_history("n2", 1); call to_history("d2", 1.); call to_history("o2", 0.)
   
write(0,*)nax, nay, nmx, nmy

  allocate(dati(nax,nay,nmx,nmy))
  allocate(dato(nmx,nmy))

  do imz=1, nmz
     call sep_read(dati)
     dato = 0.
     do iay=1, nay; ay = oay+(iay-1)*day; wty = cos(ay*PI/180.)
     do iax=1, nax; ax = oax+(iax-1)*dax; wtx = cos(ax*PI/180.)
        dato(:,:) = dato(:,:) + dati(iax,iay,:,:)/wty/wty/wtx/wtx
     end do
     end do
     call sep_write(dato)
  end do


end program
