program OW_t2d
! time-to-depth/depth-to-time conversion
  use sep
  use OW_type

  implicit none

    logical                   :: adj
    real,allocatable 	:: modl(:,:,:,:,:),data(:,:,:,:,:),vel(:,:,:),vmed(:,:,:)
    integer		:: vtype,minvx,minvy,i1,i2,i3,nta,ix,ihx,iy,ihy,iz,idt,idt1
    real			:: a,zm,tau,rt,fx,gx,b
    type(mod_grid)            :: sog
    type(velo)                :: v

  call sep_init()
  call sep_begin_prog()

  call from_param("adj",adj,.false.)

  call from_history("n1",sog%xm%n);  call from_history("d1",sog%xm%d);    call from_history("o1",sog%xm%o)
  call from_history("n2",sog%ym%n);  call from_history("d2",sog%ym%d);    call from_history("o2",sog%ym%o)
  call from_history("n3",sog%shx%n); call from_history("d3",sog%shx%d);   call from_history("o3",sog%shx%o)
  call from_history("n4",sog%shy%n); call from_history("d4",sog%shy%d);   call from_history("o4",sog%shy%o)
  if(adj) then
     call from_history("n5",sog%t%n);   call from_history("d5",sog%t%d);     call from_history("o5",sog%t%o)  
  else
     call from_history("n5",sog%z%n);   call from_history("d5",sog%z%d);     call from_history("o5",sog%z%o)  
  end if
  call from_aux("vel","n1",v%xm%n);  call from_aux("vel","d1",v%xm%d);    call from_aux("vel","o1",v%xm%o)
  call from_aux("vel","n2",v%ym%n);  call from_aux("vel","d2",v%ym%d);    call from_aux("vel","o2",v%ym%o)
  call from_aux("vel","n3",v%z%n);   call from_aux("vel","d3",v%z%d);     call from_aux("vel","o3",v%z%o)

  call from_param("vtype",vtype,1)
  allocate(vel(v%xm%n,v%ym%n,v%z%n))
  call sreed("vel",vel,4*v%xm%n*v%ym%n*v%z%n)
  allocate(vmed(sog%xm%n,sog%ym%n,v%z%n))
  minvx=(sog%xm%o-v%xm%o)/v%xm%d + 1.5
  minvy=(sog%ym%o-v%ym%o)/v%ym%d + 1.5

  if (vtype.eq.1) then
     do i3=1,v%z%n;do i2=1,sog%ym%n;do i1=1,sog%xm%n;
        vmed(i1,i2,i3)=sum(vel(minvx+i1,minvy+i2,1:i3))/i3
     end do;end do;end do
  else
     vmed=vel(minvx:sog%xm%n,minvy:sog%ym%n,:)
  end if
  deallocate(vel)

  a=-99999.; b=99999.
  do i3=1,sog%xm%n;do i2=1,sog%ym%n
     if (vmed(i3,i2,v%z%n) > a) a = vmed(i3,i2,v%z%n)
     if (vmed(i3,i2,1)     < b) b = vmed(i3,i2,1)
  end do;end do

  if (adj) then
     call from_param("dz",sog%z%d,10.);
     nta=((sog%t%n-1)*sog%t%d+sog%t%o)*a/2/sog%z%d
     call from_param("nz",sog%z%n,nta)
  else
     call from_param("dt",sog%t%d,0.008);
     nta=2*((sog%z%n-1)*sog%z%d+sog%z%o)/a/sog%t%d
     call from_param("nt",sog%t%n,nta)
     call from_param("ot",sog%t%o,0.)
  end if

  if(adj) then
     call to_history("n5",sog%z%n,"out"); call to_history("d5",sog%z%d,"out"); call to_history("o5",sog%z%o,"out")  
  else
     call to_history("n5",sog%t%n,"out"); call to_history("d5",sog%t%d,"out"); call to_history("o5",sog%t%o,"out")  
  end if
  call to_history("n4",sog%shy%n,"out"); call to_history("d4",sog%shy%d,"out");call to_history("o4",sog%shy%o,"out")
  call to_history("n2",sog%ym%n,"out");  call to_history("d2",sog%ym%d,"out"); call to_history("o2",sog%ym%o,"out")
  call to_history("n3",sog%shx%n,"out"); call to_history("d3",sog%shx%d,"out");call to_history("o3",sog%shx%o,"out")
  call to_history("n1",sog%xm%n,"out");  call to_history("d1",sog%xm%d,"out"); call to_history("o1",sog%xm%o,"out")

  allocate(data(sog%xm%n,sog%ym%n,sog%shx%n,sog%shy%n,sog%t%n),modl(sog%xm%n,sog%ym%n,sog%shx%n,sog%shy%n,sog%z%n))
  if (adj) then
     call sreed("in",data,4*sog%t%n*sog%ym%n*sog%shx%n*sog%shy%n*sog%xm%n)
  else
     call sreed("in",modl,4*sog%z%n*sog%ym%n*sog%shx%n*sog%shy%n*sog%xm%n)
  end if
  modl=0.;data=0.
  do ix=1, sog%xm%n
     do ihy=1,sog%shy%n
        do ihx=1,sog%shx%n
           do iy=1,sog%ym%n
              do iz=1,sog%z%n
                 zm  = sog%z%o + (iz-1) * sog%z%d
                 tau = 2*zm/vmed(ix,iy,iz)
                 rt  = (tau - sog%t%o) / sog%t%d
                 idt = rt
                 fx = rt - idt
                 gx = 1.0 - fx
                 idt1= idt + 1
                 if ( idt1 < sog%t%n .and. 1 < idt1) then
                    if (adj) then
                       modl(ix,iy,ihx,ihy,iz) = modl(ix,iy,ihx,ihy,iz) + gx*data(ix,iy,ihx,ihy,idt) +fx*data(ix,iy,ihx,ihy,idt1)
                    else
                       data(ix,iy,ihx,ihy,idt1) = data(ix,iy,ihx,ihy,idt1) + modl(ix,iy,ihx,ihy,iz)*fx
                       data(ix,iy,ihx,ihy,idt)  = data(ix,iy,ihx,ihy,idt)  + modl(ix,iy,ihx,ihy,iz)*gx
                    end if 
                 end if
              end do
           end do
        end do
     end do
  end do
  if (adj) then
     call srite("out",modl,4*sog%z%n*sog%ym%n*sog%shx%n*sog%shy%n*sog%xm%n)
  else
     call srite("out",data,4*sog%t%n*sog%ym%n*sog%shx%n*sog%shy%n*sog%xm%n)
  end if
  deallocate(data,modl,vmed)
  call sep_end_prog()
end program
