module owsemb
implicit none


  integer, private        :: n1,n2,n3,nrho
  real, private           :: o1,d1,o2,d2,o3,d3,orho,drho,alfa1,alfa2

  contains
!-----------------------------------------
  integer function owsem_ini(n1_i,o1_i,d1_i,n2_i,o2_i,d2_i,n3_i,nrho_i,orho_i,drho_i,alfa1_i,alfa2_i) result(stat)
  integer        :: n1_i,n2_i,n3_i,nrho_i
  real           :: o1_i,d1_i,o2_i,d2_i,o3_i,d3_i,orho_i,drho_i,alfa1_i,alfa2_i

  n1=n1_i
  n2=n2_i
  n3=n3_i
  nrho=nrho_i
  o1=o1_i
  d1=d1_i
  o2=o2_i
  d2=d2_i
  o3=o3_i
  d3=d3_i
  orho=orho_i
  drho=drho_i
  alfa1=alfa1_i
  alfa2=alfa2_i
  write(0,*) alfa1,alfa2
  stat=0
  end function
!-----------------------------------------
  function owsemb_lop(adj,add,modl,data) result(stat)
  logical     :: adj,add
  real        :: modl(:),data(:)
  integer     :: stat

  stat=1
  if (adj) modl=0
  if(not(adj)) data=0
  call owsemb_lop2(adj,add,modl,data)
  stat=0
  end function
!----------------------------------------
  subroutine owsemb_lop2(adj,add,out,in)
  logical     :: adj,add
  real        :: out(n1,nrho,n3),in(n1,n2,n3)
  integer     :: i1,i2,i3,irho,izd,idx
  real        :: gama,sgama,cgama,zm,zd,rho,rdx,fx,gx

  do i3=1,n3
     do i2=1,n2
        gama=(i2-1)*d2+o2;sgama=sind(gama)**2;cgama=cosd(gama)**2
        do irho=1,nrho
           rho=(irho-1)*drho+orho
           do i1=1,n1
              zm=(i1-1)*d1+o1
              if (i1 < int(n1/2.)) then
                 zd=(1+(rho-1)/alfa1*sgama/(alfa1**2-sgama))*zm
              else
                 zd=(1+(rho-1)/alfa2*sgama/(alfa2**2-sgama))*zm
              end if
              rdx = (zd - o1) / d1
              idx = rdx
              izd = idx + 1
              if (izd > n1-1.or. izd < 1) cycle
              fx = rdx - idx
              gx = 1.0 - fx
              if (adj) then
                 out(i1,irho,i3)=out(i1,irho,i3)+gx*in(izd,i2,i3)+fx*in(izd+1,i2,i3)
              else
                 in(izd+1,i2,i3)  = in(izd+1,i2,i3)  + out(i1,irho,i3)*fx
                 in(izd,i2,i3)  = in(izd,i2,i3)  + out(i1,irho,i3)*gx
              end if 
           end do
        end do
     end do
  end do

  end subroutine

  end module
