  program Mask_diags

  use sep

  implicit none
  
  real		  		    :: a,mu,sigma,fac
  integer                           :: i1,i2,i3,n1,n2,n3
  real                              :: o1,o2,o3,d1,d2,d3
  real,dimension(:),allocatable     :: cdp,mask

  call sep_init()
  call sep_begin_prog()

!! get info from history file
  call from_history("n1",n1);  call from_history("n2",n2);  call from_history("n3",n3)
  call from_history("d1",d1);  call from_history("d2",d2);  call from_history("d3",d3)
  call from_history("o1",o1);  call from_history("o2",o2);  call from_history("o3",o3)

  call from_param("fac",fac,0.5);

  !! write info out to history file
  call to_history("n1",n1);  call to_history("n2",n2);  call to_history("n3",n3)
  call to_history("d1",d1);  call to_history("d2",d2);  call to_history("d3",d3)
  call to_history("o1",o1);  call to_history("o2",o2);  call to_history("o3",o3)

!! allocate working space
  allocate(cdp(n1*n2),mask(n1*n2))
!!
do i3=1,n3
  call sreed("in",cdp,n1*n2*4)
  mu=sum(cdp)/size(cdp)
  a=dot_product(cdp,cdp)
  sigma=sqrt(a/size(cdp)-mu**2)
  write(0,*) ' max cdp = ',maxval(cdp),' min cdp = ',minval(cdp),' mu = ',mu,' sigma = ',sigma
  do i2=1,size(cdp)
!	if (cdp(i2).lt.mu-2*sigma) then
	if (cdp(i2).lt.mu-fac*sigma) then
		mask(i2)=0.
	else
		mask(i2)=1.
	end if
  end do
  call srite("out",mask,n1*n2*4)
end do

!! clean up
  deallocate(cdp,mask)

  call sep_end_prog()
  end program Mask_diags
