module decomap_
  use util_
  implicit none

  real,     pointer, private :: lmap(:,:,:)
  real,     pointer, private :: lrho(:,:,:)

  type(myaxis), private :: la1,la2,la3,la4,la5
  type(myaxis), private :: lb1

contains
  subroutine decomap_init(a1,a2,a3,a4,a5,b1)
    type(myaxis)     :: a1,a2,a3,a4,a5,b1

    la1=a1
    la2=a2
    la3=a3
    la4=a4
    la5=a5

    lb1=b1
    nullify(lmap,lrho)
  end subroutine decomap_init

  subroutine createmap(rho,map)
    real,target :: rho(:,:,:),map(:,:,:)
    integer     :: i1,ia2,ia3,izero,ib
    real        :: time,fzero,f

    if (la1%o > 0.) call erexit("la1%o must pass through zero")

    fzero=(-la1%o)/la1%d+1.
    izero=       nint(fzero)
    f    =fzero-floor(fzero)
    ib   =      floor(fzero)

    do ia3=1,la3%n
       do ia2=1,la2%n

          if (la1%o<0.) then          
             time=-f*la1%d/rho(izero,ia2,ia3)
             do i1=ib,1,-1 
                map(i1,ia2,ia3)=time	
                time=time-la1%d/rho(i1,ia2,ia3)
             end do
          end if

          time=(1.-f)*la1%d/rho(izero,ia2,ia3)
          do i1=ib+1,la1%n
             map(i1,ia2,ia3)=time
             time=time+la1%d/rho(i1,ia2,ia3)
          end do
       end do
    end do

    allocate(lmap(la1%n,la2%n,la3%n)); lmap=map
    allocate(lrho(la1%n,la2%n,la3%n)); lrho=rho
  end subroutine createmap

  subroutine convert(fieldin,fieldou)
    real    :: fieldin(:,:,:),fieldou(:,:,:)
    integer :: ib1,ia2,ia3
    real    :: depth,fract
    integer :: tb,ta,bi

    do ia3=1,la3%n
       do ia2=1,la2%n
          bi=1
          do ib1=1,lb1%n
             depth=lb1%o+lb1%d*(ib1-1)

             if(depth <= lmap(1,ia2,ia3)) then
                fieldou(ib1,ia2,ia3)=fieldin(1,ia2,ia3)
             else if (depth < lmap(size(lmap,1),ia2,ia3)) then

                do while( lmap(bi,ia2,ia3)<depth .and. bi<la1%n)
                   bi=bi+1
                end do
                bi=bi-1
                if(bi>=la1%n) then
                   fieldou(ib1,ia2,ia3)=fieldin(la1%n,ia2,ia3)
                else
                   fract=(depth-lmap(bi,ia2,ia3))/(lmap(bi+1,ia2,ia3)-lmap(bi,ia2,ia3))
                   fieldou(ib1,ia2,ia3)=(1.-fract)*fieldin(bi,ia2,ia3)+fract*fieldin(bi+1,ia2,ia3)
                end if
             else
                fieldou(ib1,ia2,ia3)=fieldin(la1%n,ia2,ia3)
             end if
          end do
       end do
    end do

  end subroutine convert
end module decomap_




