# 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_2d_mod
  use double_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
  subroutine init_tau_z_conv_2d(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(:)
    nz(1:min(3,size(n_z)))=n_z
    oz(1:min(3,size(o_z)))=o_z
    dz(1:min(3,size(d_z)))=d_z
    nt(1:min(3,size(n_t)))=n_t
    ot(1:min(3,size(o_t)))=o_t
    dt(1:min(3,size(d_t)))=d_t
    maxz=oz+dz*(nz-1)
    if (first(1:3).ne."NoT") then
      nullify(mapz,mapt,velz,velt)
      first(1:3)="NoT"
    end if
  end subroutine
  subroutine calc_sigma_2d(vel,sigma)
    real,target :: vel(:,:)
    real :: sigma(:,:)
    integer i1,i2
    if (size(vel,2).ne.size(sigma,2) .or. size(vel,1) .ne. size(sigma,1&
      &)) then
      call erexit("size of arrays must match")
    end if
    sigma(:,1)=2*(1./vel(:,2)-1./vel(:,1))/dz(2)
    sigma(:,nz(2))=2*(1./vel(:,nz(2))-1./vel(:,nz(2)-1))/dz(2)
    do i2=2,nz(2)-1
      sigma(:,i2)=(1./vel(:,i2+1)-1./vel(:,i2-1))/dz(2)
    end do
    sigma(1,:)=sigma(1,:)*dz(1)
    do i1=2,nz(1)
      sigma(i1,:)=sigma(i1-1,:)+ sigma(i1,:)*dz(1)
    end do
!sigma=1000000.
    if (associated(velz)) then
      deallocate(velz)
    end if
    allocate(velz(size(vel,1),size(vel,2)))
    velz=vel
  end subroutine
  subroutine tau_calc_sigma_2d(velt,sigmat)
    real :: velt(:,:),sigmat(:,:)
    integer i1,i2
    real, allocatable :: junk(:,:)
    allocate(junk(nt(1),nt(2)))
    junk(:,1)=(velt(:,2)-velt(:,1))/dt(2)
    junk(:,nt(2))=(velt(:,nt(2))-velt(:,nt(2)-1))/dt(2)
    do i2=2,nt(2)-1
      junk(:,i2)=.5*(velt(:,i2+1)-velt(:,i2-1))/dt(2)
    end do
    junk(1,:)=junk(1,:)*dt(1)
    do i1=2,nt(1)
      junk(i1,:)=junk(i1-1,:)+junk(i1,:)*dt(1)
    end do
    sigmat=-1./velt*junk
    deallocate(junk)
  end subroutine
  subroutine create_tau_map_2d(vel,map)
    real,target :: vel(:,:)
    real,target :: map(:,:)
    integer :: i1,i2,izero,ib
    real :: time,fzero,f
    if (size(vel,2).ne.size(map,2) .or. size(vel,1) .ne. size(map,1))&
      & then
      write(0,*) "vel",size(vel,1),size(vel,2)," map ",size(map,1),size&
        &(map,2)
      call erexit("size of arrays must match")
    end if
    if (oz(1) > 0.) then
      write(0,*) "CHECK THIS1",oz(1)
      call erexit("oz(1) must pass through zero")
    end if
    fzero=(-oz(1))/dz(1)+1.
    izero=nint(fzero)
    f=fzero-floor(fzero)
    ib=floor(fzero)
    do i2=1,nz(2)
      if (oz(1)<0.) then
        time=-2*f*dz(1)/vel(izero,i2)
        do i1=ib,1,-1
          map(i1,i2)=time
          time= time-2*dz(1)/vel(i1,i2)
        end do
      end if
      time=2*(1.-f)*dz(1)/vel(izero,i2)
      do i1=ib+1,nz(1)
        map(i1,i2)=time
        time= time+2*dz(1)/vel(i1,i2)
      end do
    end do
    if (associated(mapz)) then
      deallocate(mapz)
    end if
    if (associated(velz)) then
      deallocate(velz)
    end if
    allocate(mapz(size(map,1),size(map,2)))
    allocate(velz(size(vel,1),size(vel,2)))
    mapz=map
    velz=vel
  end subroutine
  subroutine convert_field_tau_z_2d(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
    real,optional,target :: map_in(:,:)
    real,pointer :: map(:,:)
    if (present(map_in)) then
      map=>map_in
    else
      map=>mapt
    end if
    if (size(field_t,2).ne.size(field_z,2).or. size(field_z,2).ne.size&
      &(map,2)) then
      write(0,*) size(map,2),size(field_z,2),size(field_t,2)
      call erexit("size axis 2 not equiv")
    end if
    nz1=size(map,1)
    nt1=size(field_t,1)
    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)) then
          field_z(i1,i2)=field_t(1,i2)
        else if (depth < map(size(map,1),i2)) then
          do while (map(bi,i2)< depth .and. bi< nz1)
            bi= bi+1
          end do
          bi= bi-1
! if(i2==50) write(0,*) "tau to z",i1,bi,depth,map(max(1,min(nz(1),bi)):max(1,min(nz(1),bi+1)),i2)
          if (bi>=nt(1)) then
            field_z(i1,i2)=field_t(nt(1),i2)
          else
            fract=(depth-map(bi,i2))/(map(bi+1,i2)-map(bi,i2))
            field_z(i1,i2)=(1.-fract)*field_t(bi,i2)+fract*field_t(bi+1&
              &,i2)
          end if
        else
          field_z(i1,i2)=field_t(nt1,i2)
        end if
      end do
    end do
  end subroutine
  subroutine convert_field_z_tau_2d(field_z,field_t,map_in)
    real :: field_z(:,:),field_t(:,:)
    integer :: i1,i2,n2
    real :: d,depth,fract
    integer :: tempi,tb,ta,bi,ia,nz1,nt1
    real,optional,target :: map_in(:,:)
    real,pointer :: map(:,:)
    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)) then
      write(0,*) size(map,2),size(field_z,2),size(field_t,2)
      call erexit("size axis 2 not equiv")
    end if
    nz1=size(map,1)
    nt1=size(field_t,1)
    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)) then
          field_t(i1,i2)=field_z(1,i2)
        else if (depth < map(size(map,1),i2)) then
          do while (map(bi,i2)< depth .and. bi< nz1)
! if(i2==1) write(0,*) "this failed",depth,bi,i1,map(bi,i2)
            bi= bi+1
          end do
          bi= bi-1
! if(i2==1) write(0,*) "z  to tau",i1,bi,depth,map(max(1,min(nt(1),bi)):max(1,min(nt(1),bi+1)),i2)
          if (bi>=nz(1)) then
            field_t(i1,i2)=field_z(nz1,i2)
          else
            fract=(depth-map(bi,i2))/(map(bi+1,i2)-map(bi,i2))
            field_t(i1,i2)=(1.-fract)*field_z(bi,i2)+fract*field_z(bi+1&
              &,i2)
          end if
        else
          field_t(i1,i2)=field_z(nz1,i2)
        end if
      end do
    end do
  end subroutine
  subroutine create_z_map_2d(vel,map)
    real,target :: vel(:,:),map(:,:)
    integer i1,i2,izero,ibr
    real :: depth,f,fzero
    if (ot(1) >0.) then
      call erexit("ot(1) must pass through zero")
    end if
    fzero=(-ot(1))/dt(1)+1.
    izero=nint(fzero)
    ibr=floor(fzero)
    f=fzero-floor(fzero)
    do i2=1,nt(2)
      if (ot(1) < 0.) then
        depth=-dt(1)*vel(izero,i2)/2.*f
        do i1=ibr,1,-1
          map(i1,i2)=depth
! if(i2==1) write(0,*) "create tau",i1,ot(1)+dt(1)*(i1-1),depth
          depth= depth-dt(1)*vel(i1,i2)/2.
        end do
      end if
      depth=dt(1)*vel(izero,i2)/2.*(1.-f)
      do i1=ibr+1,nt(1)
        map(i1,i2)=depth
! if(i2==50) write(0,*) "create tau",i1,ot(1)+dt(1)*(i1-1),depth
        depth= depth+dt(1)*vel(i1,i2)/2.
      end do
    end do
    if (associated(mapt)) then
      deallocate(mapt)
    end if
    if (associated(velt)) then
      deallocate(velt)
    end if
    allocate(mapt(size(map,1),size(map,2)))
    allocate(velt(size(vel,1),size(vel,2)))
    mapt=map
    velt=vel
  end subroutine
  subroutine convert_z_tau_2d(x,z)
    real(kind=wp) :: x(:,:),z(:,:),fractx,fractz,tempr,temp2r,dd
    integer :: iz,ix,n1,n2,i1,i2
    n2=size(z,2)
    n1=size(z,1)
    do i2=1,n2
      do i1=1,n1
        fractx=(x(i1,i2)-oz(2))/dz(2) +1.
        fractz=(z(i1,i2)-oz(1))/dz(1) +1.
        ix=floor(fractx)
        fractx= fractx-ix
        iz=floor(fractz)
        fractz= fractz-iz
        if (ix>0 .and. ix < nz(2) .and. iz > 0 .and. iz < nz(1)) then
          z(i1,i2)=(1.-fractx)*(1.-fractz)*mapz(iz,ix)+(fractx)*(1.&
            &-fractz)*mapz(iz,ix+1)+(1.-fractx)*(fractz)*mapz(iz+1,ix&
            &)+(fractx)*(fractz)*mapz(iz+1,ix+1)
        else
          if (iz<=0) then
            z(i1,i2)=oz(1)
          else if (iz>=nz(1)) then
            if (ix<1) then
              z(i1,i2)=mapz(nz(1),1)+2*distance_b(maxz(2),maxz(1),x(i1&
                &,i2),z(i1,i2))/velz(nz(1),1)
            else if (ix>=nz(2)) then
              z(i1,i2)=mapz(nz(1),nz(2))+2*distance_b(maxz(2),maxz(1),x&
                &(i1,i2),z(i1,i2))/velz(nz(1),nz(2))
            else
              tempr=(1.-fractx)*velz(nz(1),ix)+fractx*velz(nz(1),ix+1)
              temp2r=(1.-fractx)*mapz(nz(1),ix)+fractx*mapz(nz(1),ix+1&
                &)
              z(i1,i2)=temp2r+2*(z(i1,i2)-maxz(1))/tempr
            end if
          else
            if (ix<1) then
              tempr=(1.-fractz)*velz(iz,1)+fractz*velz(iz+1,1)
              temp2r=(1.-fractz)*mapz(iz,1)+fractz*mapz(iz+1,1)
              dd=oz(2)-x(i1,i2)
            else
              tempr=(1.-fractz)*velz(iz,nz(2))+fractz*velz(iz+1,nz(2))
              temp2r=(1.-fractz)*mapz(iz,nz(2))+fractz*mapz(iz+1,nz(2)&
                &)
              dd=x(i1,i2)-maxz(2)
            end if
            z(i1,i2)=temp2r+2*dd/tempr
          end if
        end if
      end do
    end do
  end subroutine
  subroutine convert_ray_traj_tau_z_2d(x,z)
    real(kind=wp) :: x(:,:),z(:,:),fractx,fractz,tempr,temp2r,dd
    integer :: iz,ix,n1,n2,i1,i2
    n1=size(x,1)
    n2=size(x,2)
    do i2=1,n2
      do i1=1,n1
        z(i1,i2)=convert_loc_to_depth_2d(real(x(i1,i2)),real(z(i1,i2))&
          &)
      end do
    end do
  end subroutine
  subroutine convert_ray_traj_z_tau_2d(x,z)
    real(kind=wp) :: x(:,:),z(:,:),fractx,fractz,tempr,temp2r,dd
    integer :: iz,ix,n1,n2,i1,i2
    n1=size(x,1)
    n2=size(x,2)
    do i2=1,n2
      do i1=1,n1
! write(0,*) "converting from",i1,real(x(i1,i2)),real(z(i1,i2))
        fractx=(x(i1,i2)-oz(2))/dz(2)+1.
        fractz=(z(i1,i2)-oz(1))/dz(1)+1.
        ix=floor(fractx)
        fractx= fractx-ix
        iz=floor(fractz)
        fractz= fractz-iz
        if (ix>0 .and. ix < nz(2) .and. iz > 0 .and. iz < nz(1)) then
          z(i1,i2)=(1.-fractx)*(1.-fractz)*mapz(iz,ix)+(fractx)*(1.&
            &-fractz)*mapz(iz,ix+1)+(1.-fractx)*(fractz)*mapz(iz+1,ix&
            &)+(fractx)*(fractz)*mapz(iz+1,ix+1)
        else
          if (iz<=0) then
            if (ix<1) then
              z(i1,i2)=mapz(nz(1),1)+2*distance_b(maxz(2),real(oz(1),wp&
                &),x(i1,i2),z(i1,i2))/velz(1,1)
            else if (ix>=nz(2)) then
              z(i1,i2)=mapz(nz(1),nz(2))+2*distance_b(maxz(2),real(oz(1&
                &),wp),x(i1,i2),z(i1,i2))/velz(1,nz(2))
            else
              if (z(i1,i2) < 0.) then
                tempr=(1.-fractx)*velz(1,ix)+fractx*velz(1,ix+1)
                temp2r=(1.-fractx)*mapz(1,ix)+fractx*mapz(1,ix+1)
                z(i1,i2)=temp2r-2*(-z(i1,i2))/tempr
              else
                z(i1,i2)=(1.-fractx)*mapz(1,ix)+(fractx)*mapz(1,ix+1)
              end if
            end if
          else if (iz>=nz(1)) then
            if (ix<1) then
              z(i1,i2)=mapz(nz(1),1)+2*distance_b(maxz(2),maxz(1),x(i1&
                &,i2),z(i1,i2))/velz(nz(1),1)
            else if (ix>=nz(2)) then
              z(i1,i2)=mapz(nz(1),nz(2))+2*distance_b(maxz(2),maxz(1),x&
                &(i1,i2),z(i1,i2))/velz(nz(1),nz(2))
            else
              tempr=(1.-fractx)*velz(nz(1),ix)+fractx*velz(nz(1),ix+1)
              temp2r=(1.-fractx)*mapz(nz(1),ix)+fractx*mapz(nz(1),ix+1&
                &)
              z(i1,i2)=temp2r+2*(z(i1,i2)-maxz(1))/tempr
            end if
          else
            if (ix<1) then
              tempr=(1.-fractz)*velz(iz,1)+fractz*velz(iz+1,1)
              temp2r=(1.-fractz)*mapz(iz,1)+fractz*mapz(iz+1,1)
              dd=oz(2)-x(i1,i2)
            else
              tempr=(1.-fractz)*velz(iz,nz(2))+fractz*velz(iz+1,nz(2))
              temp2r=(1.-fractz)*mapz(iz,nz(2))+fractz*mapz(iz+1,nz(2)&
                &)
              dd=x(i1,i2)-maxz(2)
            end if
            z(i1,i2)=temp2r+2*dd/tempr
          end if
        end if
! write(0,*) "converting to",x(i1,i2),z(i1,i2)
      end do
    end do
  end subroutine
  real function convert_loc_to_depth_2d(x,tau)
    real :: x,tau,fx,fz
    integer :: ix,iz
    fx=min(nt(2)-1.,max(1.,(x-ot(2))/dt(2)+1.))
    fz=min(nt(1)-1.,max(1.,(tau-ot(1))/dt(1)+1.))
    ix=fx
    iz=fz
    fx= fx-ix
    fz= fz-iz
    if (ix<1 .or. ix>=nt(2) .or. iz <1 .or. iz>=nt(1)) then
      write(0,*) "check",ot,dt
      write(0,*) "in convert i2",x,tau,ix,iz
      call erexit("point outside map bounds")
    end if
    convert_loc_to_depth_2d=(1.-fx)*(1.-fz)*mapt(iz,ix)+ (fx)*(1.-fz)&
      &*mapt(iz,ix+1)+(1.-fx)*(fz)*mapt(iz+1,ix)+ (fx)*(fz)*mapt(iz+1,ix+1)
  end function
  subroutine convert_ref_tau_to_z_2d(refs)
    complex :: refs(:,:)
    integer :: i1,i2
    real :: x,z
    do i2=1,size(refs,2)
      do i1=1,size(refs,1)
        x=real(refs(i1,i2))
        z=aimag(refs(i1,i2))
        refs(i1,i2)=cmplx(x,convert_loc_to_depth_2d(x,z))
      end do
    end do
  end subroutine
  subroutine convert_ref_z_to_tau_2d(refs)
    complex :: refs(:,:)
    integer :: i1,i2
    real :: x,z
    do i2=1,size(refs,2)
      do i1=1,size(refs,1)
        x=real(refs(i1,i2))
        z=aimag(refs(i1,i2))
        refs(i1,i2)=cmplx(x,convert_loc_to_tau_2d(x,z))
      end do
    end do
  end subroutine
  real function convert_loc_to_tau_2d(x,z)
    real :: x,z,fx,fz
    integer :: ix,iz
    fx=min(nz(2)-1.,max(1.,(x-oz(2))/dz(2)+1.))
    fz=min(nz(1)-1.,max(1.,(z-oz(1))/dz(1)+1.))
    ix=fx
    iz=fz
    fx= fx-ix
    fz= fz-iz
    if (ix<1 .or. ix>=nz(2) .or. iz <1 .or. iz>=nz(1)) then
      write(0,*) "check",oz,dz
      call erexit("point  outside map bounds")
    end if
    convert_loc_to_tau_2d=(1.-fx)*(1.-fz)*mapz(iz,ix)+ (fx)*(1.-fz)*mapz&
      &(iz,ix+1)+(1.-fx)*(fz)*mapz(iz+1,ix)+ (fx)*(fz)*mapz(iz+1,ix+1)
  end function
  subroutine calc_vrms_2d(vint,vrms)
    real vint(:,:),vrms(:,:)
    integer :: i1,i2
    real :: tot
    do i2=1,size(vint,2)
      tot=0
      do i1=1,size(vint,1)
        tot= tot+vint(i1,i2)**2
        vrms(i1,i2)=sqrt(tot/i1)
      end do
    end do
  end subroutine
  logical function map_z_2d(mymap)
    real :: mymap(:,:)
    map_z_2d=.true.
    if (associated(mapz)) then
      mymap=mapz
    else
      map_z_2d=.false.
    end if
  end function
  logical function map_t_2d(mymap)
    real :: mymap(:,:)
    map_t_2d=.true.
    if (associated(mapt)) then
      mymap=mapt
    else
      map_t_2d=.false.
    end if
  end function
end module
