# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!<
!tau_z_conv_mod
!
!Usage:
!call tau_z_conv_mod (arg1, arg2 arg3)
!
!Input Parameters:
!
!
!Output Parameters:
!
!
!
!Description:
!
!
!
!>
!-------------------------------------------------
!
!Author:
!
!Date Created:Wed Jan 20 9:19:10 PST 1999
!
!Purpose:
!
!
module tau_z_conv_3d_mod
  use double_mod
  use put_3d_mod
  implicit none
  integer, save, private :: nz(3),nt(3)
  real, save, private :: dz(3),oz(3),ot(3),dt(3)
  real(kind=wp), save, private :: maxz(3)
  real,pointer,private :: mapz(:,:,:),mapt(:,:,:),velz(:,:,:),velt(:,:&
    &,:)
  character(len=10) :: first
  contains
  logical function init_tau_z_conv_3d(n_z,o_z,d_z,n_t,o_t,d_t)
    real :: o_t(:),d_t(:),o_z(:),d_z(:)
    integer :: n_t(:),n_z(:)
    init_tau_z_conv_3d=.false.
    nz(1:size(n_z))=n_z
    oz(1:size(o_z))=o_z
    dz(1:size(d_z))=d_z
    nt(1:size(n_t))=n_t
    ot(1:size(o_t))=o_t
    dt(1:size(d_t))=d_t
    maxz=oz+dz*(nz-1)
    if (first(1:3).ne."NoT") then
      first(1:3)="NoT"
      nullify(mapz,mapt,velz,velt)
    end if
    init_tau_z_conv_3d=.true.
  end function
  logical function create_tau_map_3d(vel,map)
    real,target :: vel(:,:,:)
    real,target :: map(:,:,:)
    integer :: i1,i2,izero,ib,i3
    real :: time,fzero,f
    create_tau_map_3d=.false.
    if (size(vel,2).ne.size(map,2) .or. size(vel,1) .ne. size(map,1)&
      & .or. size(vel,3).ne.size(map,3)) then
      write(0,*) "vel",size(vel,1),size(vel,2)," map ",size(map,1),size&
        &(map,2)
      return
    end if
    if (oz(1) > 0.) then
      write(0,*) "oz(1) must pass through zero"
      return
    end if
    fzero=(-oz(1))/dz(1)+1.
    izero=nint(fzero)
    f=fzero-floor(fzero)
    ib=floor(fzero)
    do i3=1,nz(3)
      do i2=1,nz(2)
        if (oz(1)<0.) then
          time=-2*f*dz(1)/vel(izero,i2,i3)
          do i1=ib,1,-1
            map(i1,i2,i3)=time
            time= time-2*dz(1)/vel(i1,i2,i3)
          end do
        end if
        time=2*(1.-f)*dz(1)/vel(izero,i2,i3)
        do i1=ib+1,nz(1)
          map(i1,i2,i3)=time
          time= time+2*dz(1)/vel(i1,i2,i3)
        end do
      end do
    end do
!allocate(mapz(size(vel,1),size(vel,2),size(vel,3)))
!allocate(velz(size(vel,1),size(vel,2),size(vel,3)))
    mapz=>map
    velz=>vel
    create_tau_map_3d=.true.
  end function
  logical function set_z_vel_map(vel,map)
    real, target :: vel(:,:,:),map(:,:,:)
    mapz=>map
    velz=>vel
    set_z_vel_map=.true.
  end function
  logical function set_tau_vel_map(vel,map)
    real, target :: vel(:,:,:),map(:,:,:)
    mapt=>map
    velt=>vel
    set_tau_vel_map=.true.
  end function
  logical function calc_sigmax_3d(vel,sigmax)
    real :: vel(:,:,:)
    real :: sigmax(:,:,:)
    integer i1,i2,i3
    calc_sigmax_3d=.false.
    if (size(vel,2).ne.size(sigmax,2) .or. size(vel,1) .ne. size(sigmax&
      &,1) .or. size(vel,3).ne.size(sigmax,3) ) then
      write(0,*) "size of arrays must match"
      return
    end if
    do i3=1,nz(3)
      sigmax(:,1,i3)=2*(1./vel(:,2,i3)-1./vel(:,1,i3))/dz(2)
      sigmax(:,nz(2),i3)=2*(1./vel(:,nz(2),i3)-1./vel(:,nz(2)-1,i3))/dz&
        &(2)
      do i2=2,nz(2)-1
        sigmax(:,i2,i3)=(1./vel(:,i2+1,i3)-1./vel(:,i2-1,i3))/dz(2)
      end do
      sigmax(1,:,i3)=sigmax(1,:,i3)*dz(1)
      do i1=2,nz(1)
        sigmax(i1,:,i3)=sigmax(i1-1,:,i3)+ sigmax(i1,:,i3)*dz(1)
      end do
    end do
    calc_sigmax_3d=.true.
  end function
  logical function calc_sigmay_3d(vel,sigmay)
    real,target :: vel(:,:,:)
    real :: sigmay(:,:,:)
    integer i1,i2,i3
    calc_sigmay_3d=.false.
    if (size(vel,2).ne.size(sigmay,2) .or. size(vel,1) .ne. size(sigmay&
      &,1) .or. size(vel,3).ne.size(sigmay,3) ) then
      write(0,*) "size of arrays must match"
      return
    end if
    do i2=1,nz(2)
      sigmay(:,i2,1)=2*(1./vel(:,i2,2)-1./vel(:,i2,1))/dz(3)
      sigmay(:,i2,nz(3))=2*(1./vel(:,i2,nz(3))-1./vel(:,i2,nz(3)-1))/dz&
        &(3)
      do i3=3,nz(3)-1
        sigmay(:,i2,i3)=(1./vel(:,i2,i3+1)-1./vel(:,i2,i3-1))/dz(3)
      end do
      sigmay(1,i2,:)=sigmay(1,i2,:)*dz(1)
      do i1=2,nz(1)
        sigmay(i1,i2,:)=sigmay(i1-1,i2,:)+ sigmay(i1,i2,:)*dz(1)
      end do
    end do
    calc_sigmay_3d=.true.
  end function
  logical function convert_field_tau_z_3d(field_t,field_z,map_in)
    real :: field_t(:,:,:),field_z(:,:,:)
    integer :: i1,i2,n2
    real :: d,depth,fract
    integer :: tempi,tb,ta,bi,ia,nz1,nt1,i3
    real,optional,target :: map_in(:,:,:)
    real,pointer :: map(:,:,:)
    if (present(map_in)) then
      map=>map_in
    else
      map=>mapt
    end if
    convert_field_tau_z_3d=.false.
    if (size(field_t,2).ne.size(field_z,2).or. size(field_z,2).ne.size&
      &(map,2) .or. size(field_t,3).ne.size(field_z,3)) then
      write(0,*) size(map,2),size(field_z,2),size(field_t,2)
      write(0,*) "size axis 2 not equiv"
      return
    end if
    nz1=size(map,1)
    nt1=size(field_z,1)
    do i3=1,size(field_t,3)
      do i2=1,size(field_t,2)
        bi=1
        do i1=1,size(field_z,1)
          depth=oz(1)+dz(1)*(i1-1)
          if (depth <= map(1,i2,i3)) then
            field_z(i1,i2,i3)=field_t(1,i2,i3)
          else if (depth < map(size(map,1),i2,i3)) then
            do while (map(bi,i2,i3)< depth .and. bi< nz1)
              bi= bi+1
            end do
            bi= bi-1
            if (bi>=nz(1)) then
              field_z(i1,i2,i3)=field_t(nz1,i2,i3)
            else
              fract=(depth-map(bi,i2,i3))/(map(bi+1,i2,i3)-map(bi,i2,i3&
                &))
              field_z(i1,i2,i3)=(1.-fract)*field_t(bi,i2,i3)+fract*field_t&
                &(bi+1,i2,i3)
            end if
          else
            field_z(i1,i2,i3)=field_t(nz1,i2,i3)
          end if
        end do
      end do
    end do
    convert_field_tau_z_3d=.true.
  end function
  logical function convert_field_z_tau_3d(field_z,field_t,map_in)
    real :: field_z(:,:,:),field_t(:,:,:)
    integer :: i1,i2,n2,i3
    real :: d,depth,fract
    integer :: tempi,tb,ta,bi,ia,nz1,nt1
    real,optional,target :: map_in(:,:,:)
    real,pointer :: map(:,:,:)
    convert_field_z_tau_3d=.false.
    if (present(map_in)) then
      map=>map_in
    else
      map=>mapz
    end if
    if (size(field_t,2).ne.size(field_z,2).or. size(field_z,2).ne.size&
      &(map,2) .or. size(field_t,3).ne.size(field_z,3)) then
      write(0,*) size(map,2),size(field_z,2),size(field_t,2)
      write(0,*) "SIZE 3",size(field_t,3),size(field_z,3)
      write(0,*) "size axis 2 not equiv"
      return
    end if
    nz1=size(map,1)
    nt1=size(field_t,1)
    do i3=1,size(field_z,3)
      do i2=1,size(field_z,2)
        bi=1
        do i1=1,size(field_t,1)
          depth=ot(1)+dt(1)*(i1-1)
          if (depth <= map(1,i2,i3)) then
            field_t(i1,i2,i3)=field_z(1,i2,i3)
          else if (depth < map(size(map,1),i2,i3)) then
            do while (map(bi,i2,i3)< depth .and. bi< nz1)
              bi= bi+1
            end do
            bi= bi-1
            if (bi>=nz(1)) then
              field_t(i1,i2,i3)=field_z(nz1,i2,i3)
            else
              fract=(depth-map(bi,i2,i3))/(map(bi+1,i2,i3)-map(bi,i2,i3&
                &))
              field_t(i1,i2,i3)=(1.-fract)*field_z(bi,i2,i3)+fract*field_z&
                &(bi+1,i2,i3)
            end if
          else
            field_t(i1,i2,i3)=field_z(nz1,i2,i3)
          end if
        end do
      end do
    end do
    convert_field_z_tau_3d=.true.
  end function
  logical function create_z_map_3d(vel,map)
    real,target :: vel(:,:,:),map(:,:,:)
    integer i1,i2,izero,ibr,i3
    real :: depth,f,fzero
    create_z_map_3d=.false.
    if (ot(1) >0.) then
      write(0,*) "ot(1) must pass through zero"
      return
    end if
    fzero=(-ot(1))/dt(1)+1.
    izero=nint(fzero)
    ibr=floor(fzero)
    f=fzero-floor(fzero)
!write(0,*) "check it",f,fzero,izero,ot(1),dt(1)
    do i3=1,nt(3)
      do i2=1,nt(2)
        if (ot(1) < 0.) then
          depth=-dt(1)*vel(izero,i2,i3)/2.*f
          do i1=ibr,1,-1
            map(i1,i2,i3)=depth
            depth= depth-dt(1)*vel(i1,i2,i3)/2.
          end do
        end if
        depth=dt(1)*vel(izero,i2,i3)/2.*(1.-f)
        do i1=ibr+1,nt(1)
          map(i1,i2,i3)=depth
          depth= depth+dt(1)*vel(i1,i2,i3)/2.
        end do
      end do
    end do
    create_z_map_3d=.true.
  end function
!OLD
  subroutine convert_ray_traj_z_tau_2d(rays)
    real :: rays(:,:),fractx,fractz,tempr,temp2r,dd,fracty
    integer :: iz,ix,n1,n2,i1,i2,iy
    real :: y
    call init_coord(nz,oz,dz)
    y=oz(2)+dz(2)*.5
    n1=size(rays,1)
    n2=size(rays,2)
    if (n1.ne.3) then
      call erexit("expecting first axis of rays to be length 3")
    end if
    do i2=1,n2
      rays(3,i2)=get_val_lin_pos(mapz,rays(3,i2),rays(1,i2),y)
    end do
  end subroutine
  subroutine convert_ray_traj_z_tau_3d(rays)
    real :: rays(:,:),fractx,fractz,tempr,temp2r,dd,fracty
    integer :: iz,ix,n1,n2,i1,i2,iy
    call init_coord(nz,oz,dz)
    n1=size(rays,1)
    n2=size(rays,2)
    if (n1.ne.3) then
      call erexit("expecting first axis of rays to be length 3")
    end if
    do i2=1,n2
      rays(3,i2)=get_val_lin_pos(mapz,rays(3,i2),rays(1,i2),rays(2,i2)&
        &)
    end do
  end subroutine
  real function convert_loc_to_depth_3d(x,y,tau)
    real :: x,tau,fx,fz,y
    integer :: ix,iz
    convert_loc_to_depth_3d=get_val_lin_pos(mapt,tau,x,y)
  end function
  logical function convert_ref_z_to_tau_2d(refs)
    complex :: refs(:,:)
    integer :: i1,i2,i3
    real :: x,y,z
    convert_ref_z_to_tau_2d=.false.
    call init_coord(nz,oz,dz)
    y=ot(2)+dt(2)*.5
    do i3=1,size(refs,2)
      do i1=1,size(refs,1)
        x=real(refs(i1,i3))
        z=aimag(refs(i1,i3))
        refs(i1,i3)=cmplx(x,get_val_lin_pos(mapz,z,x,y))
      end do
    end do
    convert_ref_z_to_tau_2d=.true.
  end function
  logical function convert_ref_tau_to_z_2d(refs,oref,dref,reject_val)
    complex :: refs(:,:)
    integer :: i1,i2,i3
    real :: x,y,z,oref(2),dref(2)
    real, optional :: reject_val
    logical :: reject
    real :: reject_it
    convert_ref_tau_to_z_2d=.false.
    if (present(reject_val)) then
      reject=.true.
      reject_it=reject_val
    else
      reject=.false.
    end if
    call init_coord(nz,oz,dz)
    y=ot(2)+dt(2)*.5
    do i3=1,size(refs,2)
      do i1=1,size(refs,1)
        x=real(refs(i1,i3))
        z=aimag(refs(i1,i3))
        if (.not.(reject .and. aimag(refs(i1,i3)).eq.reject_it)) then
          refs(i1,i3)=cmplx(x,get_val_lin_pos(mapt,z,x,y))
        end if
      end do
    end do
    convert_ref_tau_to_z_2d=.true.
  end function
  logical function convert_ref_tau_to_z_3d(refs,oref,dref,reject_val)
    real :: refs(:,:,:)
    integer :: i1,i2,i3
    real :: x,y,z,oref(2),dref(2)
    real, optional :: reject_val
    logical :: reject
    real :: reject_it
    convert_ref_tau_to_z_3d=.false.
    if (present(reject_val)) then
      reject=.true.
      reject_it=reject_val
    else
      reject=.false.
    end if
    call init_coord(nz,oz,dz)
    do i3=1,size(refs,3)
      do i2=1,size(refs,2)
        y=oref(2)+dref(2)*(i2-1)
        do i1=1,size(refs,1)
          x=oref(1)+dref(1)*(i1-1)
          if (.not.(reject .and. refs(i1,i2,i3).eq.reject_it)) then
            refs(i1,i2,i3)=get_val_lin_pos(mapt,refs(i1,i2,i3),x,y)
          end if
        end do
      end do
    end do
    convert_ref_tau_to_z_3d=.true.
  end function
  logical function convert_ref_z_to_tau_3d(refs,oref,dref,reject_val)
    real :: refs(:,:,:)
    integer :: i1,i2,i3
    real :: x,y,z,oref(2),dref(2)
    real, optional :: reject_val
    logical :: reject
    real :: reject_it
    convert_ref_z_to_tau_3d=.false.
    if (present(reject_val)) then
      reject=.true.
      reject_it=reject_val
    else
      reject=.false.
    end if
    call init_coord(nz,oz,dz)
    do i3=1,size(refs,3)
      do i2=1,size(refs,2)
        y=oref(2)+dref(2)*(i2-1)
        do i1=1,size(refs,1)
          x=oref(1)+dref(1)*(i1-1)
          if (.not.(reject .and. refs(i1,i2,i3).eq.reject_it)) then
            refs(i1,i2,i3)=get_val_lin_pos(mapz,refs(i1,i2,i3),x,y)
          end if
        end do
      end do
    end do
    convert_ref_z_to_tau_3d=.true.
  end function
  real function convert_loc_to_tau_3d(x,y,z)
    real :: x,z,fx,fz,y
    integer :: ix,iz
    call init_coord(nz,oz,dz)
    convert_loc_to_tau_3d=get_val_lin_pos(mapz,z,x,y)
  end function
  subroutine calc_vrms_3d(vint,vrms)
    real vint(:,:,:),vrms(:,:,:)
    integer :: i1,i2,i3
    real :: tot
    do i3=1,size(vint,3)
      do i2=1,size(vint,2)
        tot=0
        do i1=1,size(vint,1)
          tot= tot+vint(i1,i2,i3)**2
          vrms(i1,i2,i3)=sqrt(tot/i1)
        end do
      end do
    end do
  end subroutine
  logical function map_z_3d(mymap)
    real :: mymap(:,:,:)
    map_z_3d=.true.
    if (associated(mapz)) then
      mymap=mapz
    else
      map_z_3d=.false.
    end if
  end function
  logical function map_t_3d(mymap)
    real :: mymap(:,:,:)
    map_t_3d=.true.
    if (associated(mapt)) then
      mymap=mapt
    else
      map_t_3d=.false.
    end if
  end function
end module
