program OW_sqrvel

use sep

implicit none
real			:: o1,o3,d1,d3,v0,dvel,xcenter,zcenter,extentx,extentz
integer			:: i1,i3,n1,n3,nano
real,allocatable		:: velin(:,:)

call sep_init()
call sep_begin_prog()

 call from_param("nx",n1,2001); call from_param("dx",d1,10.); call from_param("ox",o1,-10000.)
 call from_param("nz",n3,201); call from_param("dz",d3,10.); call from_param("oz",o3,0.)
 call to_history("n1",n1,"slow_out"); call to_history("d1",d1,"slow_out"); call to_history("o1",o1,"slow_out")
 call to_history("n2",1,"slow_out"); call to_history("d2",1.,"slow_out"); call to_history("o2",0.,"slow_out")
 call to_history("n3",n3,"slow_out"); call to_history("d3",d3,"slow_out"); call to_history("o3",o3,"slow_out")
 call to_history("n4",1,"slow_out"); call to_history("d4",1.,"slow_out"); call to_history("o4",0.,"slow_out")
 call to_history("eout",4,"slow_out");call to_history("label3","","slow_out")

 call from_param("back_vel",v0,2000.)
 call from_param("vel_pert",dvel,0.05)
 call from_param("xcenter",xcenter,0.)
 call from_param("zcenter",zcenter,750.)
 call from_param("extentx",extentx,500.)
 call from_param("extentz",extentz,500.)
 call from_param("nano",nano,2)

allocate(velin(n1,n3));velin=v0;
do i3=1,n3
   do i1=1,n1
      if (i1>=int(((xcenter-extentx/2)-o1)/d1+1.5) .and.&
          i1<=int(((xcenter+extentx/2)-o1)/d1+1.5) .and.&
          i3>=int(((zcenter-extentz/2)-o3)/d3+1.5) .and.&
          i3<=int(((zcenter+extentz/2)-o3)/d3+1.5)) then
          velin(i1,i3)=v0*(1+dvel)
      end if
   end do
end do
if(nano == 2) then
do i3=1,n3
   do i1=1,n1
      if (i1>=int(((-xcenter-extentx/2)-o1)/d1+1.5) .and.&
          i1<=int(((-xcenter+extentx/2)-o1)/d1+1.5) .and.&
          i3>=int(((zcenter-extentz/2)-o3)/d3+1.5) .and.&
          i3<=int(((zcenter+extentz/2)-o3)/d3+1.5)) then
          velin(i1,i3)=v0*(1-dvel)
      end if
   end do
end do
end if

velin=1/velin
call srite("slow_out",velin,4*n1*n3)
deallocate(velin)
call sep_end_prog()
end program
