program owscale

use sep
implicit none

real,allocatable                    :: data(:,:,:,:)
real, allocatable                   :: a(:,:)
integer                             :: i1,i2,i3,i4,i5,i6,n1,n2,n3,n4,n5,n6,node
real                                :: d1,d2,d3,d4,d5,d6,o1,o2,o3,o4,o5,o6
integer, external                   :: omp_get_thread_num,omp_get_num_threads

call sep_init()
call sep_begin_prog()
call from_param("node",node,14)
call from_history("n1",n1); call from_history("d1",d1); call from_history("o1",o1)
call from_history("n2",n2); call from_history("d2",d2); call from_history("o2",o2)
call from_history("n3",n3); call from_history("d3",d3); call from_history("o3",o3)
call from_history("n4",n4); call from_history("d4",d4); call from_history("o4",o4)

call to_history("n1",n1); call to_history("d1",d1); call to_history("o1",o1)
call to_history("n2",n2); call to_history("d2",d2); call to_history("o2",o2)
call to_history("n3",n3); call to_history("d3",d3); call to_history("o3",o3)
call to_history("n4",n4); call to_history("d4",d4); call to_history("o4",o4)
allocate(data(n1,n2,n3,n4),a(n2,n4))
call sreed("w",a, 4*n4*n2)

        data=0.
        call sreed("in", data, 4*n1*n2*n3*n4)
        !$OMP parallel do private (i1,i2,i3,i4)  schedule(dynamic,1) num_threads(node) ordered
        do i4=1,n4
           do i2=1,n2
              if(a(i2,i4).ne.0)  data(:,i2,:,i4)=data(:,i2,:,i4)/a(i2,i4)
           end do
        end do
       !$OMP END PARALLEL DO
       call srite("out", data, 4*n1*n2*n3*n4)

 deallocate(data,a)
 call sep_end_prog()
 end program
