program owscale

use sep
implicit none
logical                             :: rms
real,allocatable                    :: data(:,:,:,:), a(:,:),temp(:,:)
integer                             :: i1,i2,i3,i4,i5,n1,n2,n3,n4,n5,wx,wy
real                                :: d1,d2,d3,d4,d5,o1,o2,o3,o4,o5,perc,b

call sep_init()
call sep_begin_prog()
call from_param("rms",rms,.true.)
call from_param("wx",wx,5)
call from_param("perc",perc,0.1)
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)
call to_history("eout",4);call to_history("label3","")
allocate(data(n1,n2,n3,n4),a(n1,n4),temp(n1,n4))
call sreed("in", data, 4*n1*n2*n3*n4)

  a=0
  do i4=1,n4;do i2=1,n2;do i1=1,n1
     a(i1,i4) = a(i1,i4) + sum(data(i1,i2,:,i4)**2)
  end do;end do;end do
  a = sqrt(a/n1/n2/n3)
  temp=0;
  do i4=1,n4;do i1=1,n1
     temp(i1,i4)=sum(a(max(1,i1-wx):min(n1,i1+wx),i4))/(min(n1,i1+wx)-max(1,i1-wx)+1)
  end do;end do
  
  do i4=1,n4
     b=perc*maxval(a(:,i4))
     do i1=1,n1
        if(a(i1,i4).gt.b) then
          data(i1,:,:,i4)=data(i1,:,:,i4)/a(i1,i4)
        else
          data(i1,:,:,i4)=0
        end if
     end do
  end do

call to_history("n1",n1,"scale"); call to_history("d1",d1,"scale"); call to_history("o1",o1,"scale")  
call to_history("n2",n4,"scale"); call to_history("d2",d4,"scale"); call to_history("o2",o4,"scale")  
call to_history("n3",1,"scale") ; call to_history("d3",1.,"scale") ; call to_history("o3",1.,"scale")
call to_history("esize",4,"scale")
call srite("scale",a,4*n1*n4)

call srite("out", data, 4*n1*n2*n3*n4)
deallocate(data,a,temp)
call sep_end_prog()
end program
