program dso
use sep
use igrad1
use curvature
implicit none

logical						:: pwstk,zeronear,mco
real						:: a,b
real, allocatable					:: wgcig(:,:,:),h(:),odcig(:,:,:),lapxo(:)
integer						:: i1,i2,i3,stat1
integer                                                     :: n1,n2,n3,n4
real                                                        :: d1,d2,d3,d4,o1,o2,o3,o4

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 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",1); call to_history("d4",1.); call to_history("o4",0.)
call to_history("eout",4);call to_history("label3","")
call from_param("a",a,1.)
call from_param("pwstk",pwstk,.false.)
call from_param("zeronear",zeronear,.false.)
call from_param("mco",mco,.false.)
write(0,*) "mco ",mco
if (pwstk) call from_param("b",b,1.);write(0,*) pwstk,b
allocate(wgcig(n1,n2,n3),h(n2))
allocate(odcig(n1,n2,n3),lapxo(n2))
call sreed("in",wgcig,4*n1*n2*n3)
do i2=1,n2
   h(i2)=(i2-1)*d2+o2
   h(i2)=abs(h(i2)/((n2-1)*d2+o2))
end do

do i3=1,n3;do i1=1,n1
   wgcig(i1,:,i3)=wgcig(i1,:,i3)*h**a
end do;end do

if (mco) then
   do i3=1,n3
      do i1=1,n1
          stat1 = igrad1_lop (  .false.,  .false., wgcig(i1,:,i3), lapxo)   ! tmp = grad y
          stat1 = igrad1_lop (  .true. ,  .false.,  odcig(i1,:,i3), lapxo)   ! x = x + grad tmp
          stat1 = igrad1_lop (  .false.,  .false., wgcig(i1,:,i3), lapxo)   ! tmp = grad y
          odcig(i1,:,i3) = odcig(i1,:,i3) * (1 + lapxo**2)**(1.5) / ((1 + lapxo**2)**3 + maxval(1+lapxo**2))
      end do
   end do
   wgcig=odcig
!   do i3=1,n3
!      do i1=1,n1
!         stat1 = igrad1_lop (  .false.,  .false., odcig(i1,:,i3), lapxo)   ! tmp = grad x
!         stat1 = igrad1_lop (  .true. ,   .true., wgcig(i1,:,i3), lapxo)   ! y = y + grad tmp
!      end do
!   end do
!   do i3=1,n3
!      do i2=1,n2
!         do i1=1,n1
!            if (i2 == 1) then
!              odcig(i1,i2,i3) = (wgcig(i1,i2+1,i3)-wgcig(i1,i2,i3))/d2
!            else if (i2 == n2) then
!              odcig(i1,i2,i3) = (wgcig(i1,i2,i3)-wgcig(i1,i2-1,i3))/d2
!            else
!              odcig(i1,i2,i3) = (wgcig(i1,i2+1,i3)-2*wgcig(i1,i2,i3)+wgcig(i1,i2-1,i3))/(d2**2)
!            end if
!         end do
!      end do
!   end do
end if


do i3=1,n3;do i1=1,n1
   wgcig(i1,:,i3)=wgcig(i1,:,i3)*h**a
end do;end do

call srite("out",wgcig,4*n1*n2*n3)
deallocate(wgcig,h,lapxo,odcig)
call sep_end_prog()
end program
