 program wderiv

 use sep
 use igrad1

 implicit none

 type axis
    integer 					:: n
    real    					:: o, d
 end type axis

 real, allocatable				:: green(:,:),tandip(:,:),dxgreen(:,:),dzgreen(:,:)
 integer					:: it,ix,iz,stat
 real						:: eps,pi
 type(axis)					:: x,z,t

 call sep_init()
 call sep_begin_prog()

 call from_history("n1",z%n); call from_history("d1",z%d); call from_history("o1",z%o)
 call from_history("n2",x%n); call from_history("d2",x%d); call from_history("o2",x%o)  
 call from_history("n3",t%n); call from_history("d3",t%d); call from_history("o3",t%o)

 call from_param("eps",eps,0.0001)

 call to_history("n1",z%n); call to_history("d1",z%d); call to_history("o1",z%o)
 call to_history("n2",x%n); call to_history("d2",x%d); call to_history("o2",x%o)  
 call to_history("n3",t%n); call to_history("d3",t%d); call to_history("o3",t%o)  
 call to_history("n4",1); call to_history("d4",1); call to_history("o4",1)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 pi = atan(1.)*4
  allocate(green(z%n,x%n),dxgreen(z%n,x%n),dzgreen(z%n,x%n),tandip(z%n,x%n))

  do it=1,t%n
     if (mod(it,250).eq.0) write(0,*) " time = ",it," from ",t%n
     call sreed("in", green, 4*z%n*x%n)
     dxgreen=0.;dzgreen=0.
     do iz=1,z%n
        stat = igrad1_lop (.true.,.false.,dxgreen(iz,:),(green(iz,:)))
     end do
     do ix=1,x%n
        stat = igrad1_lop (.true.,.false.,dzgreen(:,ix),(green(:,ix)))
     end do
     do ix=1,x%n
        do iz=1,z%n
           if (dzgreen(iz,ix)==0 ) then
              if (dxgreen(iz,ix) < 0) then
                   tandip(iz,ix)=-huge(1.)
              else if (dxgreen(iz,ix)>0) then
                   tandip(iz,ix)=huge(1.)
              else
                   tandip(iz,ix)=0.
              end if
           else
              tandip(iz,ix)=dxgreen(iz,ix)*dzgreen(iz,ix)/(dzgreen(iz,ix)**2+eps)
           end if
        end do
     end do
     tandip=atand(tandip)
     call srite("out",tandip,4*z%n*x%n)
  end do
;dzgreen=0.
  do iz=1,x%n
     if (mod(it,250).eq.0) write(0,*) " time = ",it," from ",t%n
     call sreed_window("in",3,(/z%n,x%n,t%n/),(/1,x%n,t%n/),(/iz-1,0,0/),(/1,1,1/),4,green)
     dxgreen=0.;dtgreen=0.
     do it=1,t%n
        stat = igrad1_lop (.true.,.false.,dxgreen(:,it),(green(:,it)))
     end do
     do ix=1,x%n
        stat = igrad1_lop (.true.,.false.,dtgreen(ix,:),(green(ix,:)))
     end do
     do it=1,t%n
        do ix=1,x%n
           if (dxgreen(ix,it)==0 ) then
              if (dtgreen(ix,it) < 0) then
                   dtdx(ix,it)=-huge(1.)
              else if (dxgreen(iz,ix)>0) then
                   tandip(iz,ix)=huge(1.)
              else
                   tandip(iz,ix)=0.
              end if
           else
              tandip(iz,ix)=dxgreen(iz,ix)*dzgreen(iz,ix)/(dzgreen(iz,ix)**2+eps)
           end if
        end do
     end do
     tandip=atand(tandip)
  call srite_window("out",3,(/n1,n2,n3/),(/n1,n2,1/),(/0,0,i3-1/),(/1,1,1/),4,dado)
!     call srite("out",tandip,4*z%n*x%n)
  end do

  deallocate(green,dxgreen,dzgreen,tandip)

 call sep_end_prog()
end program wderiv
