program OW_3dmask

use sep

implicit none
character                     :: label1,label2,label3,label4,label5
real			:: o1,o2,o3,o4,o5,d1,d2,d3,d4,d5,om1,om2,dm1,dm2
integer			:: i1,i2,i3,i4,i5,n1,n2,n3,n4,n5,nm1,nm2
real, allocatable	          :: mask(:,:)
complex,allocatable		:: data(:,:,:,:)

call sep_init()
call sep_begin_prog()

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 from_history("n5",n5); call from_history("d5",d5); call from_history("o5",o5)
write(0,*) n1,n2,n3,n4,n5
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("n5",n5); call to_history("d5",d5); call to_history("o5",o5)
call to_history("eout",8)
call to_history("label1","cmp_x"); call to_history("label2","cmp_y"); call to_history("label3","hx")
call to_history("label4","hy"); call to_history("label5","freq");
call from_aux("mask","n1",nm1); call from_aux("mask","d1",dm1); call from_aux("mask","o1",om1)
call from_aux("mask","n2",nm2); call from_aux("mask","d2",dm2); call from_aux("mask","o2",om2)

if (nm1.ne.n1.and.nm2.ne.n2.and.om1.ne.o1.and.om2.ne.o2.and.dm1.ne.d1.and.dm2.ne.d2) call seperr ("mask and data in different grids")
allocate(data(n1,n2,n3,n4));data=0;
allocate(mask(n1,n2));mask=0;
call sreed("mask",mask,4*n1*n2)
do i1=1,n1
   where(mask(i1,:).ne.1) mask(i1,:)=0
end do

   do i5=1,n5
      call sreed("in",data,8*n1*n2*n3*n4)
      do i4=1,n4
         do i3=1,n3
            data(:,:,i3,i4)=data(:,:,i3,i4)*mask
         end do
      end do
      call srite("out",data,8*n1*n2*n3*n4)
   end do

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