module ztconv_

  use double_mod
  implicit none
  character(len=10):: first
  real(kind=wp), save, private :: maxz(3)

  integer,       save, private :: nz(3),nt(3)
  real,          save, private :: dz(3),dt(3)
  real,          save, private :: oz(3),ot(3)

  real,        pointer,private :: mapz(:,:),mapt(:,:)
  real,        pointer,private :: velz(:,:),velt(:,:)

contains

  subroutine init_ztconv(n_z,o_z,d_z,n_t,o_t,d_t)
    real    :: o_t(:),o_z(:)
    real    :: d_t(:),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 init_ztconv

  subroutine createmap_z(vel,map)
    real,target :: vel(:,:),map(:,:)
    integer :: i1,i2,izero,ib
    real    :: time,fzero,f

    if (oz(1) > 0.) call erexit("oz(1) must pass through zero")

    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)) deallocate(mapz)
    if (associated(velz)) deallocate(velz)
    allocate(mapz(size(map,1),size(map,2))); mapz=map
    allocate(velz(size(vel,1),size(vel,2))); velz=vel

  end subroutine createmap_z

  subroutine createmap_t(vel,map)
    real,target :: vel(:,:),map(:,:)
    integer :: i1,i2,izero,ib
    real    :: depth,fzero,f

    if (ot(1) >0.) call erexit("ot(1) must pass through zero")

    fzero=(-ot(1))/dt(1)+1.
    izero=        nint(fzero)	
    f     =fzero-floor(fzero)
    ib    =      floor(fzero)

    do i2=1,nt(2)
      if (ot(1) < 0.) then
        depth=-dt(1)*vel(izero,i2)/2.*f
        do i1=ib,1,-1 
          map(i1,i2)=depth
          depth=depth-dt(1)*vel(i1,i2)/2.
        end do
      end if
      depth=dt(1)*vel(izero,i2)/2.*(1.-f)

      do i1=ib+1,nt(1)
        map(i1,i2)=depth
        depth=depth+dt(1)*vel(i1,i2)/2.
      end do
    end do 

    if (associated(mapt)) deallocate(mapt)
    if (associated(velt)) deallocate(velt)
    allocate(mapt(size(map,1),size(map,2))); mapt=map
    allocate(velt(size(vel,1),size(vel,2))); velt=vel

  end subroutine createmap_t

  subroutine convert_t2z(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 (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 convert_t2z

  subroutine convert_z2t(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)
            bi=				 bi+1
          end do 
          bi=			bi-1
          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 convert_z2t

  subroutine zzmap(vel,map)
    real,target :: vel(:,:),map(:,:)
    integer :: i1,i2,izero,ib
    real    :: time,fzero,f

    if (oz(1) > 0.) call erexit("oz(1) must pass through zero")

    fzero=(-oz(1))/dz(1)*4+1.
    izero=       nint(fzero)
    f    =fzero-floor(fzero)
    ib   =      floor(fzero)

    do i2=1,nz(2)
      if (oz(1)<0.) then
        time=-f*dz(1)*4/vel(izero,i2)
        do i1=ib,1,-1 
          map(i1,i2)=time	
          time=time-dz(1)*4/vel(i1,i2)
        end do
      end if
      time=(1.-f)*dz(1)*4/vel(izero,i2)
      do i1=ib+1,nz(1)
        map(i1,i2)=time
        time=time+dz(1)*4/vel(i1,i2)
      end do
    end do 

    if (associated(mapz)) deallocate(mapz)
    if (associated(velz)) deallocate(velz)
    allocate(mapz(size(map,1),size(map,2))); mapz=map
    allocate(velz(size(vel,1),size(vel,2))); velz=vel

  end subroutine zzmap

end module ztconv_




