mod_trac.t

Go to the documentation of this file.

module mod_trac
  use mod_global_parameters
  use mod_mhd
  use mod_functions_bfield, only: mag
 
  implicit none
  private
 
  ! common
  integer :: numFL,numLP
  double precision :: dL,Tmax,trac_delta,T_bott
  double precision, allocatable :: xFi(:,:)
  ! TRACB
  integer :: numxT^D
  double precision :: dxT^D 
  double precision :: xTmin(ndim),xTmax(ndim)
  double precision, allocatable :: xT(:^D&,:)
  ! TRACL mask
  integer, allocatable :: trac_grid(:),ground_grid(:)
  integer :: ngrid_trac,ngrid_ground
  logical, allocatable :: trac_pe(:)
 
  public :: initialize_trac_after_settree
 
contains
 
  subroutine init_trac_line(mask)
    logical, intent(in) :: mask 
    integer :: refine_factor,ix^D,ix(ndim),j,iFL,numL(ndim),finegrid
    double precision :: lengthFL
    double precision :: xprobmin(ndim),xprobmax(ndim),domain_nx(ndim)
    integer :: memxFi
 
    trac_delta=0.25d0
    !----------------- init magnetic field -------------------!
    refine_factor=2**(refine_max_level-1)
    ^d&xprobmin(^d)=xprobmin^d\
    ^d&xprobmax(^d)=xprobmax^d\
    ^d&domain_nx(^d)=domain_nx^d\
    dl=(xprobmax(ndim)-xprobmin(ndim))/(domain_nx(ndim)*refine_factor)
    ! max length of a field line
    if (.not. mask) phys_trac_mask=xprobmax^nd
    {^iftwod
      lengthfl=(phys_trac_mask-xprobmin(2))*3.d0
    }
    {^ifthreed
      lengthfl=(phys_trac_mask-xprobmin(3))*3.d0
    }
    
    numlp=floor(lengthfl/dl)
    numl=1
    numfl=1
    do j=1,ndim-1
      ! number of field lines, every 4 finest grid, every direction
      finegrid=mhd_trac_finegrid
      numl(j)=floor((xprobmax(j)-xprobmin(j))/dl/finegrid)
      numfl=numfl*numl(j)
    end do
    allocate(xfi(numfl,ndim))
    xfi(:,ndim)=xprobmin(ndim)+dl/50.d0
    {do ix^db=1,numl(^db)\}
      ^d&ix(^d)=ix^d\
      ifl=0
      do j=ndim-1,1,-1
        ifl=ifl+(ix(j)-(ndim-1-j))*(numl(j))**(ndim-1-j)
      end do
      xfi(ifl,1:ndim-1)=xprobmin(1:ndim-1)+finegrid*dl*ix(1:ndim-1)-finegrid*dl/2.d0
    {end do\}
    {^iftwod
    if(mype .eq. 0) write(*,*) 'NOTE:  2D TRAC method take the y-dir == grav-dir'
    }
    {^ifthreed
    if(mype .eq. 0) write(*,*) 'NOTE:  3D TRAC method take the z-dir == grav-dir'
    }
 
    memxfi=floor(8*numfl*numlp*ndim/1e6)
    if (mype==0) write(*,*) 'Memory requirement for each processor in TRAC:'
    if (mype==0) write(*,*)  memxfi,' MB'
 
    allocate(trac_pe(0:npe-1),trac_grid(max_blocks),ground_grid(max_blocks))
 
  end subroutine init_trac_line
 
  subroutine init_trac_block(mask)
    logical, intent(in) :: mask 
    integer :: refine_factor,finegrid,iFL,j
    integer :: ix^D,ixT^L
    integer :: numL(ndim),ix(ndim)
    double precision :: lengthFL
    double precision :: ration,a0
    double precision :: xprobmin(ndim),xprobmax(ndim),dxT(ndim)
 
    refine_factor=2**(refine_max_level-1)
    ^d&xprobmin(^d)=xprobmin^d\
    ^d&xprobmax(^d)=xprobmax^d\
    ^d&dxt^d=(xprobmax^d-xprobmin^d)/(domain_nx^d*refine_factor/block_nx^d)\
    ^d&dxt(ndim)=dxt^d\
    finegrid=mhd_trac_finegrid
    {^ifoned
      dl=dxt1/finegrid
    }
    {^nooned
      dl=min(dxt^d)/finegrid
    }
    ! table for interpolation
    ^d&xtmin(^d)=xprobmin^d\
    ^d&xtmax(^d)=xprobmax^d\
    if(mask) xtmax(ndim)=phys_trac_mask
    ! max length of a field line
    if(mask) then
      lengthfl=maxval(xtmax-xprobmin)*3.d0
    else
      lengthfl=maxval(xprobmax-xprobmin)*3.d0
    end if
    numlp=floor(lengthfl/dl)
    ^d&numxt^d=ceiling((xtmax(^d)-xtmin(^d)-smalldouble)/dxt^d)\
    allocate(xt(numxt^d,ndim))
    ixtmin^d=1\
    ixtmax^d=numxt^d\
    {do j=1,numxt^d
      xt(j^d%ixT^s,^d)=(j-0.5d0)*dxt^d+xtmin(^d)
    end do\}
    if(mask) xtmax(ndim)=maxval(xt(:^d&,ndim))+half*dxt(ndim)
    numl=1
    numfl=1
    do j=1,ndim-1
      ! number of field lines, every 4 finest grid, every direction
      numl(j)=floor((xprobmax(j)-xprobmin(j))/dl)
      numfl=numfl*numl(j)
    end do
    allocate(xfi(numfl,ndim))
    xfi(:,ndim)=xprobmin(ndim)+dl/50.d0
    {do ix^db=1,numl(^db)\}
      ^d&ix(^d)=ix^d\
      ifl=0
      do j=ndim-1,1,-1
        ifl=ifl+(ix(j)-(ndim-1-j))*(numl(j))**(ndim-1-j)
      end do
      xfi(ifl,1:ndim-1)=xprobmin(1:ndim-1)+dl*ix(1:ndim-1)-dl/2.d0
    {end do\}
    {^iftwod
    if(mype .eq. 0) write(*,*) 'NOTE:  2D TRAC method take the y-dir == grav-dir'
    }
    {^ifthreed
    if(mype .eq. 0) write(*,*) 'NOTE:  3D TRAC method take the z-dir == grav-dir'
    }
  end subroutine init_trac_block
 
  subroutine trac_simple(tco_global,trac_alfa,T_peak)
    double precision, intent(in) :: tco_global, trac_alfa,T_peak
    integer :: iigrid, igrid
 
    do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
      {^ifoned
      ps(igrid)%special_values(1)=tco_global
      }
      if(ps(igrid)%special_values(1)<trac_alfa*ps(igrid)%special_values(2)) then
        ps(igrid)%special_values(1)=trac_alfa*ps(igrid)%special_values(2)
      end if
      if(ps(igrid)%special_values(1) .lt. t_bott) then
        ps(igrid)%special_values(1)=t_bott
      else if(ps(igrid)%special_values(1) .gt. 0.2d0*t_peak) then
        ps(igrid)%special_values(1)=0.2d0*t_peak
      end if
      ps(igrid)%wextra(ixg^t,iw_tcoff)=ps(igrid)%special_values(1)
      !> special values(2) to save old tcutoff
      ps(igrid)%special_values(2)=ps(igrid)%special_values(1)
    end do
  end subroutine trac_simple
 
  subroutine ltrac(T_peak)
    double precision, intent(in) :: T_peak
    integer :: iigrid, igrid
    integer :: ixO^L,trac_tcoff
 
    ixo^l=ixm^ll^ladd1;
    trac_tcoff=iw_tcoff
    do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
      where(ps(igrid)%wextra(ixo^s,trac_tcoff) .lt. t_bott)
        ps(igrid)%wextra(ixo^s,trac_tcoff)=t_bott
      else where(ps(igrid)%wextra(ixo^s,trac_tcoff) .gt. 0.2d0*t_peak)
        ps(igrid)%wextra(ixo^s,trac_tcoff)=0.2d0*t_peak
      end where
    end do
  end subroutine ltrac
 
  subroutine tracl(mask,T_peak)
    logical, intent(in) :: mask
    double precision, intent(in) :: T_peak
    integer :: ix^D,j
    integer :: iigrid, igrid
    double precision :: xF(numFL,numLP,ndim)
    integer :: numR(numFL),ix^L
    double precision :: Tlcoff(numFL)
    integer :: ipel(numFL,numLP),igridl(numFL,numLP)
    logical :: forwardl(numFL)
 
    tmax=t_peak
    xf=zero
    
    do j=1,ndim
      xf(:,1,j)=xfi(:,j)
    end do
 
    call mpi_barrier(icomm,ierrmpi)
    ! record pe and grid for mask region
    call update_pegrid()
    ! get temperature for the calculation of Te gradient, 
    ! which will be stored in wextra(ixI^S,Tcoff_) temporarily
    call get_te_grid()
    ! find out direction for tracing B field
    call get_btracing_dir(ipel,igridl,forwardl)
    ! trace Bfield and get Tcoff for each field line
    call get_tcoff_line(xf,numr,tlcoff,ipel,igridl,forwardl,mask)
    ! init vairable Tcoff_ and Tweight_
    call init_trac_tcoff()
    ! get cell Tcoff via interpolation
    call interp_tcoff(xf,ipel,igridl,numr,tlcoff)
    ! convert primitive data back to conserved data
    call mpi_barrier(icomm,ierrmpi)
 
 
  end subroutine tracl
 
  subroutine tracb(mask,T_peak)
    logical, intent(in) :: mask
    double precision, intent(in) :: T_peak
    integer :: peArr(numxT^D),gdArr(numxT^D),numR(numFL)
    double precision :: Tcoff(numxT^D),Tcmax(numxT^D),Bdir(numxT^D,ndim)
    double precision :: xF(numFL,numLP,ndim),Tcoff_line(numFL)
    integer :: xpe(numFL,numLP,2**ndim)
    integer :: xgd(numFL,numLP,2**ndim)
 
    tmax=t_peak
    tcoff=zero
    tcmax=zero
    bdir=zero
    pearr=-1
    gdarr=-1
    call block_estable(mask,tcoff,tcmax,bdir,pearr,gdarr)
    xf=zero
    numr=0
    tcoff_line=zero
    call block_trace_mfl(mask,tcoff,tcoff_line,tcmax,bdir,pearr,gdarr,xf,numr,xpe,xgd)
    call block_interp_grid(mask,xf,numr,xpe,xgd,tcoff_line)
  end subroutine tracb
 
  subroutine block_estable(mask,Tcoff,Tcmax,Bdir,peArr,gdArr)
    logical :: mask
    double precision :: Tcoff(numxT^D),Tcoff_recv(numxT^D)
    double precision :: Tcmax(numxT^D),Tcmax_recv(numxT^D)
    double precision :: Bdir(numxT^D,ndim),Bdir_recv(numxT^D,ndim)
    integer :: peArr(numxT^D),peArr_recv(numxT^D)
    integer :: gdArr(numxT^D),gdArr_recv(numxT^D)
    integer :: xc^L,xd^L,ix^D
    integer :: iigrid,igrid,numxT,intab
    double precision :: xb^L
 
    tcoff_recv=zero
    tcmax_recv=zero
    bdir_recv=zero
    pearr_recv=-1
    gdarr_recv=-1
    !> combine table from different processors 
    xcmin^d=nghostcells+1\
    xcmax^d=block_nx^d+nghostcells\
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      ps(igrid)%wextra(:^d&,tweight_)=zero
      ps(igrid)%wextra(:^d&,tcoff_)=zero
      ^d&xbmin^d=rnode(rpxmin^d_,igrid)-xtmin(^d)\
      ^d&xbmax^d=rnode(rpxmax^d_,igrid)-xtmin(^d)\
      xdmin^d=nint(xbmin^d/dxt^d)+1\
      xdmax^d=ceiling((xbmax^d-smalldouble)/dxt^d)\
      {do ix^d=xdmin^d,xdmax^d \}
        intab=0
        {if (ix^d .le. numxt^d) intab=intab+1 \}
        if(intab .eq. ndim) then
          !> in principle, no overlap will happen here
          tcoff(ix^d)=max(tcoff(ix^d),ps(igrid)%special_values(1))
          tcmax(ix^d)=ps(igrid)%special_values(2)
          !> abs(Bdir) <= 1, so that Bdir+2 should always be positive
          bdir(ix^d,1:ndim)=ps(igrid)%special_values(3:3+ndim-1)+2.d0
          pearr(ix^d)=mype
          gdarr(ix^d)=igrid
        end if
      {end do\}
    end do 
    call mpi_barrier(icomm,ierrmpi)
    numxt={numxt^d*}
    call mpi_allreduce(pearr,pearr_recv,numxt,mpi_integer,&
                       mpi_max,icomm,ierrmpi)
    call mpi_allreduce(gdarr,gdarr_recv,numxt,mpi_integer,&
                       mpi_max,icomm,ierrmpi)
    call mpi_allreduce(tcoff,tcoff_recv,numxt,mpi_double_precision,&
                       mpi_max,icomm,ierrmpi)
    call mpi_allreduce(bdir,bdir_recv,numxt*ndim,mpi_double_precision,&
                       mpi_max,icomm,ierrmpi)
    if(.not. mask) then
      call mpi_allreduce(tcmax,tcmax_recv,numxt,mpi_double_precision,&
                         mpi_max,icomm,ierrmpi)
    end if
    pearr=pearr_recv
    gdarr=gdarr_recv
    tcoff=tcoff_recv
    bdir=bdir_recv-2.d0
    if(.not. mask) tcmax=tcmax_recv
  end subroutine block_estable
 
  subroutine block_trace_mfl(mask,Tcoff,Tcoff_line,Tcmax,Bdir,peArr,gdArr,xF,numR,xpe,xgd)
    integer :: i,j,k,k^D,ix_next^D
    logical :: mask,flag,first
    double precision :: Tcoff(numxT^D),Tcoff_line(numFL)
    double precision :: Tcmax(numxT^D),Tcmax_line(numFL)
    double precision :: xF(numFL,numLP,ndim)
    integer :: ix_mod(ndim,2),numR(numFL)
    double precision :: alfa_mod(ndim,2)
    double precision :: nowpoint(ndim),nowgridc(ndim)
    double precision :: Bdir(numxT^D,ndim)
    double precision :: init_dir,now_dir1(ndim),now_dir2(ndim)
    integer :: peArr(numxT^D),xpe(numFL,numLP,2**ndim)
    integer :: gdArr(numxT^D),xgd(numFL,numLP,2**ndim)
 
    do i=1,numfl
      flag=.true.
      ^d&k^d=ceiling((xfi(i,^d)-xtmin(^d)-smalldouble)/dxt^d)\
      tcoff_line(i)=tcoff(k^d)
      if(.not. mask) tcmax_line(i)=tcmax(k^d)
      ix_next^d=k^d\
      j=1
      xf(i,j,:)=xfi(i,:)
      do while(flag)
        nowpoint(:)=xf(i,j,:)
        nowgridc(:)=xt(ix_next^d,:)
        first=.true.
        if(j .eq. 1) then 
          call rk_bdir(nowgridc,nowpoint,ix_next^d,now_dir1,bdir,&
                       ix_mod,first,init_dir)
        else
          call rk_bdir(nowgridc,nowpoint,ix_next^d,now_dir1,bdir,&
                       ix_mod,first)
        end if
        {^iftwod
          xgd(i,j,1)=gdarr(ix_mod(1,1),ix_mod(2,1))
          xgd(i,j,2)=gdarr(ix_mod(1,2),ix_mod(2,1))
          xgd(i,j,3)=gdarr(ix_mod(1,1),ix_mod(2,2))
          xgd(i,j,4)=gdarr(ix_mod(1,2),ix_mod(2,2))
          xpe(i,j,1)=pearr(ix_mod(1,1),ix_mod(2,1))
          xpe(i,j,2)=pearr(ix_mod(1,2),ix_mod(2,1))
          xpe(i,j,3)=pearr(ix_mod(1,1),ix_mod(2,2))
          xpe(i,j,4)=pearr(ix_mod(1,2),ix_mod(2,2))
        }
        {^ifthreed
          xgd(i,j,1)=gdarr(ix_mod(1,1),ix_mod(2,1),ix_mod(3,1))
          xgd(i,j,2)=gdarr(ix_mod(1,2),ix_mod(2,1),ix_mod(3,1))
          xgd(i,j,3)=gdarr(ix_mod(1,1),ix_mod(2,2),ix_mod(3,1))
          xgd(i,j,4)=gdarr(ix_mod(1,2),ix_mod(2,2),ix_mod(3,1))
          xgd(i,j,5)=gdarr(ix_mod(1,1),ix_mod(2,1),ix_mod(3,2))
          xgd(i,j,6)=gdarr(ix_mod(1,2),ix_mod(2,1),ix_mod(3,2))
          xgd(i,j,7)=gdarr(ix_mod(1,1),ix_mod(2,2),ix_mod(3,2))
          xgd(i,j,8)=gdarr(ix_mod(1,2),ix_mod(2,2),ix_mod(3,2))
          xpe(i,j,1)=pearr(ix_mod(1,1),ix_mod(2,1),ix_mod(3,1))
          xpe(i,j,2)=pearr(ix_mod(1,2),ix_mod(2,1),ix_mod(3,1))
          xpe(i,j,3)=pearr(ix_mod(1,1),ix_mod(2,2),ix_mod(3,1))
          xpe(i,j,4)=pearr(ix_mod(1,2),ix_mod(2,2),ix_mod(3,1))
          xpe(i,j,5)=pearr(ix_mod(1,1),ix_mod(2,1),ix_mod(3,2))
          xpe(i,j,6)=pearr(ix_mod(1,2),ix_mod(2,1),ix_mod(3,2))
          xpe(i,j,7)=pearr(ix_mod(1,1),ix_mod(2,2),ix_mod(3,2))
          xpe(i,j,8)=pearr(ix_mod(1,2),ix_mod(2,2),ix_mod(3,2))
        }
        nowpoint(:)=nowpoint(:)+init_dir*now_dir1*dl
        {if(nowpoint(^d) .gt. xtmax(^d) .or. nowpoint(^d) .lt. xtmin(^d)) then
          flag=.false.
        end if\}
        if(mask .and. nowpoint(ndim) .gt. phys_trac_mask) then
          flag=.false.
        end if
        if(flag) then
          first=.false.
          ^d&ix_next^d=ceiling((nowpoint(^d)-xtmin(^d)-smalldouble)/dxt^d)\
          nowgridc(:)=xt(ix_next^d,:)
          call rk_bdir(nowgridc,nowpoint,ix_next^d,now_dir2,bdir,&
                       ix_mod,first)
          xf(i,j+1,:)=xf(i,j,:)+init_dir*dl*half*(now_dir1+now_dir2)
          {if(xf(i,j+1,^d) .gt. xtmax(^d) .or. xf(i,j+1,^d) .lt. xtmin(^d)) then
            flag=.false.
          end if\}
          if(mask .and. xf(i,j+1,ndim) .gt. phys_trac_mask) then
            flag=.false.
          end if
          if(flag) then
            ^d&ix_next^d=ceiling((xf(i,j+1,^d)-xtmin(^d)-smalldouble)/dxt^d)\
            j=j+1
            tcoff_line(i)=max(tcoff_line(i),tcoff(ix_next^d))
            if(.not.mask) tcmax_line(i)=max(tcmax_line(i),tcmax(ix_next^d))
          end if
        end if
      end do
      numr(i)=j
      if(mask) then
        if(tcoff_line(i) .gt. tmax*0.2d0) then
          tcoff_line(i)=tmax*0.2d0
        end if
      else
        if(tcoff_line(i) .gt. tcmax_line(i)*0.2d0) then
          tcoff_line(i)=tcmax_line(i)*0.2d0
        end if
      end if
    end do
  end subroutine block_trace_mfl
 
  subroutine rk_bdir(nowgridc,nowpoint,ix_next^D,now_dir,Bdir,ix_mod,first,init_dir)
    double precision :: nowpoint(ndim),nowgridc(ndim)
    integer :: ix_mod(ndim,2)
    double precision :: alfa_mod(ndim,2)
    integer :: ix_next^D,k^D
    double precision :: now_dir(ndim)
    double precision :: Bdir(numxT^D,ndim)
    logical :: first
    double precision, optional :: init_dir
 
    {if(nowpoint(^d) .gt. xtmin(^d)+half*dxt^d .and. nowpoint(^d) .lt. xtmax(^d)-half*dxt^d) then
      if(nowpoint(^d) .le. nowgridc(^d)) then
        ix_mod(^d,1)=ix_next^d-1
        ix_mod(^d,2)=ix_next^d
        alfa_mod(^d,1)=abs(nowgridc(^d)-nowpoint(^d))/dxt^d
        alfa_mod(^d,2)=one-alfa_mod(^d,1)
      else
        ix_mod(^d,1)=ix_next^d
        ix_mod(^d,2)=ix_next^d+1
        alfa_mod(^d,2)=abs(nowgridc(^d)-nowpoint(^d))/dxt^d
        alfa_mod(^d,1)=one-alfa_mod(^d,2)
      end if
    else
      ix_mod(^d,:)=ix_next^d
      alfa_mod(^d,:)=half
    end if\}
    now_dir=zero
    {^iftwod
    do k1=1,2
    do k2=1,2
      now_dir=now_dir + bdir(ix_mod(1,k1),ix_mod(2,k2),:)*alfa_mod(1,k1)*alfa_mod(2,k2)
    end do
    end do
    } 
    {^ifthreed
    do k1=1,2
    do k2=1,2
    do k3=1,2
      now_dir=now_dir + bdir(ix_mod(1,k1),ix_mod(2,k2),ix_mod(3,k3),:)&
             *alfa_mod(1,k1)*alfa_mod(2,k2)*alfa_mod(3,k3)
    end do
    end do
    end do
    }
    if(present(init_dir)) then
      init_dir=sign(one,now_dir(ndim))
    end if
  end subroutine rk_bdir
 
  subroutine block_interp_grid(mask,xF,numR,xpe,xgd,Tcoff_line)
    logical :: mask
    double precision :: xF(numFL,numLP,ndim)
    integer :: numR(numFL)
    integer :: xpe(numFL,numLP,2**ndim)
    integer :: xgd(numFL,numLP,2**ndim)
    double precision :: Tcoff_line(numFL)
    double precision :: weightIndex,weight,ds
    integer :: i,j,k,igrid,iigrid,ixO^L,ixc^L,ixc^D
    double precision :: dxMax^D,dxb^D
 
    !> interpolate lines into grid cells
    weightindex=2.d0
    dxmax^d=2.d0*dl;
    ixo^l=ixm^ll;
    do i=1,numfl
      do j=1,numr(i) 
        do k=1,2**ndim
          if(mype .eq. xpe(i,j,k)) then
            igrid=xgd(i,j,k)
            if(igrid .le. igrids(igridstail)) then
              ^d&dxb^d=rnode(rpdx^d_,igrid)\
              ^d&ixcmin^d=floor((xf(i,j,^d)-dxmax^d-ps(igrid)%x(ixomin^dd,^d))/dxb^d)+ixomin^d\
              ^d&ixcmax^d=floor((xf(i,j,^d)+dxmax^d-ps(igrid)%x(ixomin^dd,^d))/dxb^d)+ixomin^d\
              {if (ixcmin^d<ixomin^d) ixcmin^d=ixomin^d\}
              {if (ixcmax^d>ixomax^d) ixcmax^d=ixomax^d\}
              {do ixc^d=ixcmin^d,ixcmax^d\}
                ds=0.d0
                {ds=ds+(xf(i,j,^d)-ps(igrid)%x(ixc^dd,^d))**2\}
                ds=sqrt(ds)
                if(ds .le. 0.099d0*dl) then
                  weight=(1/(0.099d0*dl))**weightindex
                else
                  weight=(1/ds)**weightindex
                endif
                ps(igrid)%wextra(ixc^d,tweight_)=ps(igrid)%wextra(ixc^d,tweight_)+weight
                ps(igrid)%wextra(ixc^d,tcoff_)=ps(igrid)%wextra(ixc^d,tcoff_)+weight*tcoff_line(i)
              {enddo\}
            else
              call mpistop("we need to check here 366Line in mod_trac.t")
            end if
          end if
        end do
      end do
    end do
    ! finish interpolation
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      where (ps(igrid)%wextra(ixo^s,tweight_)>0.d0)
        ps(igrid)%wextra(ixo^s,tcoff_)=ps(igrid)%wextra(ixo^s,tcoff_)/ps(igrid)%wextra(ixo^s,tweight_)
      elsewhere
        ps(igrid)%wextra(ixo^s,tcoff_)=0.2d0*tmax
      endwhere
    enddo
  end subroutine block_interp_grid
 
   ! TODO remove, not used? 
!  subroutine get_Tmax_grid(x,w,ixI^L,ixO^L,Tmax_grid)
!    integer, intent(in) :: ixI^L,ixO^L
!    double precision, intent(in) :: x(ixI^S,1:ndim)
!    double precision, intent(out) :: w(ixI^S,1:nw)
!    double precision :: Tmax_grid
!    double precision :: pth(ixI^S),Te(ixI^S)
!
!    call mhd_get_pthermal(w,x,ixI^L,ixO^L,pth)
!    Te(ixO^S)=pth(ixO^S)/w(ixO^S,rho_)
!    Tmax_grid=maxval(Te(ixO^S))
!  end subroutine get_Tmax_grid  
 
  subroutine init_trac_tcoff()
    integer :: ixI^L,ixO^L,igrid,iigrid
 
    ixi^l=ixg^ll;
    ixo^l=ixm^ll;
 
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      ps(igrid)%wextra(ixi^s,tcoff_)=0.d0
      ps(igrid)%wextra(ixi^s,tweight_)=0.d0
    end do
 
  end subroutine init_trac_tcoff
 
  subroutine update_pegrid()
 
    ngrid_trac=0
    ngrid_ground=0
    trac_pe(:)=.false.
    trac_grid(:)=0
    ground_grid(:)=0
 
    ! traverse gridtable to information of grid inside mask region
    call traverse_gridtable()
 
  end subroutine update_pegrid
 
  subroutine traverse_gridtable()
    use mod_global_parameters
 
    double precision :: dxb^D,xb^L
    integer :: iigrid,igrid,j
    logical, allocatable :: trac_pe_recv(:)
    double precision :: hcmax_bt
 
    allocate(trac_pe_recv(npe))
 
    hcmax_bt=xprobmin^nd+(xprobmax^nd-xprobmin^nd)/(domain_nx^nd*2**(refine_max_level-1))
 
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      xbmin^nd=rnode(rpxmin^nd_,igrid)
      if (xbmin^nd<phys_trac_mask) then
        trac_pe(mype)=.true.
        ngrid_trac=ngrid_trac+1
        trac_grid(ngrid_trac)=igrid
        if (xbmin^nd<hcmax_bt) then
          ngrid_ground=ngrid_ground+1
          ground_grid(ngrid_ground)=igrid
        endif
      endif
    enddo
 
    call mpi_allreduce(trac_pe,trac_pe_recv,npe,mpi_logical,mpi_lor,icomm,ierrmpi)
    trac_pe=trac_pe_recv
 
    deallocate(trac_pe_recv)
 
  end subroutine traverse_gridtable
 
  subroutine get_te_grid()
    integer :: ixI^L,ixO^L,igrid,iigrid,j
 
    ixi^l=ixg^ll;
    ixo^l=ixm^ll;
 
    do j=1,ngrid_trac
      igrid=trac_grid(j)
      
      call mhd_get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixi^l,ps(igrid)%wextra(ixi^s,tcoff_))
      !TODO move  outside loop for optimziation? 
      if(has_equi_rho0) then
        ps(igrid)%wextra(ixi^s,tcoff_)=ps(igrid)%wextra(ixi^s,tcoff_)/&
          (ps(igrid)%w(ixi^s,rho_) + ps(igrid)%equi_vars(ixi^s,equi_rho0_,0))
      else
        ps(igrid)%wextra(ixi^s,tcoff_)=ps(igrid)%wextra(ixi^s,tcoff_)/ps(igrid)%w(ixi^s,rho_)
      endif
    enddo
 
  end subroutine get_te_grid
 
  subroutine get_btracing_dir(ipel,igridl,forwardl)
    integer :: ipel(numFL,numLP),igridl(numFL,numLP)
    logical :: forwardl(numFL)
 
    integer :: igrid,ixO^L,iFL,j,ix^D,idir,ixb^D,ixbb^D
    double precision :: xb^L,dxb^D,xd^D,factor,Bh
    integer :: numL(ndim),ixmin(ndim),ixmax(ndim),ix(ndim)
    logical :: forwardRC(numFL)
 
    ipel=-1
    igridl=0
    forwardl=.true.
 
    ^d&ixomin^d=ixmlo^d;
    ^d&ixomax^d=ixmhi^d;
 
    do j=1,ngrid_ground
      igrid=ground_grid(j)
      ^d&dxb^d=rnode(rpdx^d_,igrid);
      ^d&xbmin^d=rnode(rpxmin^d_,igrid);
      ^d&xbmax^d=rnode(rpxmax^d_,igrid);
      ^d&ixmin(^d)=floor((xbmin^d-xprobmin^d)/(mhd_trac_finegrid*dl))+1;
      ^d&ixmax(^d)=floor((xbmax^d-xprobmin^d)/(mhd_trac_finegrid*dl));
      ^d&numl(^d)=floor((xprobmax^d-xprobmin^d)/(mhd_trac_finegrid*dl));
      ixmin(^nd)=1
      ixmax(^nd)=1
      numl(^nd)=1
 
      {do ix^db=ixmin(^db),ixmax(^db)\}
        ^d&ix(^d)=ix^d;
        ifl=0
        do idir=ndim-1,1,-1
          ifl=ifl+(ix(idir)-(ndim-1-idir))*(numl(idir))**(ndim-1-idir)
        enddo
        ipel(ifl,1)=mype
        igridl(ifl,1)=igrid
 
        ^d&ixb^d=floor((xfi(ifl,^d)-ps(igrid)%x(ixomin^dd,^d))/dxb^d)+ixomin^d;
        ^d&xd^d=(xfi(ifl,^d)-ps(igrid)%x(ixb^dd,^d))/dxb^d;
        bh=0.d0
        {do ixbb^d=0,1\}
          factor={abs(1-ix^d-xd^d)*}
          if (b0field) then
            bh=bh+factor*(ps(igrid)%w(ixb^d+ixbb^d,mag(^nd))+ps(igrid)%B0(ixb^d+ixbb^d,^nd,0))
          else
            bh=bh+factor*ps(igrid)%w(ixb^d+ixbb^d,mag(^nd))
          endif
        {enddo\}
        if (bh>0) then
          forwardl(ifl)=.true.
        else
          forwardl(ifl)=.false.
        endif
      {enddo\}
    enddo
 
    call mpi_allreduce(forwardl,forwardrc,numfl,mpi_logical,&
                         mpi_land,icomm,ierrmpi)
    forwardl=forwardrc
 
  end subroutine get_btracing_dir
 
  subroutine get_tcoff_line(xFL,numR,TcoffFL,ipeFL,igridFL,forwardFL,mask)
    use mod_trace_field
 
    double precision :: xFL(numFL,numLP,ndim)
    integer :: numR(numFL)
    double precision :: TcoffFL(numFL),TmaxFL(numFL)
    integer :: ipeFL(numFL,numLP),igridFL(numFL,numLP)
    logical :: forwardFL(numFL)
    logical, intent(in) :: mask
 
    integer :: nwP,nwL,iFL,iLP
    double precision :: wPm(numFL,numLP,2),wLm(numFL,1+2)
    character(len=std_len) :: ftype,tcondi
 
    nwp=2
    nwl=2
    ftype='Bfield'
    tcondi='TRAC'    
    call trace_field_multi(xfl,wpm,wlm,dl,numfl,numlp,nwp,nwl,forwardfl,ftype,tcondi)
    do ifl=1,numfl
      numr(ifl)=int(wlm(ifl,1))
      tcofffl(ifl)=wlm(ifl,2)
      tmaxfl(ifl)=wlm(ifl,3)
      if(mask) then
        if(tcofffl(ifl)>0.2d0*tmax) tcofffl(ifl)=0.2d0*tmax
      else
        tmaxfl(ifl)=wlm(ifl,3)
        if(tcofffl(ifl)>0.2d0*tmaxfl(ifl)) tcofffl(ifl)=0.2d0*tmaxfl(ifl)
      end if
 
      if(tcofffl(ifl)<t_bott) tcofffl(ifl)=t_bott
    enddo
 
    do ifl=1,numfl
      if (numr(ifl)>0) then
        do ilp=1,numr(ifl)
          ipefl(ifl,ilp)=int(wpm(ifl,ilp,1))
          igridfl(ifl,ilp)=int(wpm(ifl,ilp,2))
        enddo
      endif
    enddo
 
 
  end subroutine get_tcoff_line
 
  subroutine interp_tcoff(xF,ipel,igridl,numR,Tlcoff)
    double precision :: xF(numFL,numLP,ndim)
    integer :: numR(numFL),ipel(numFL,numLP),igridl(numFL,numLP)
    double precision :: Tlcoff(numFL)
 
    integer :: iFL,iLP,ixO^L,ixI^L,ixc^L,ixb^L,ixc^D
    integer :: igrid,j,ipmin,ipmax,igrid_nb
    double precision :: dxb^D,dxMax^D,xb^L,Tcnow
    double precision :: xFnow(ndim)
    integer :: weightIndex,idn^D,ixmax^ND
    double precision :: ds,weight
 
    weightindex=2
 
    ^d&iximin^d=ixglo^d;
    ^d&iximax^d=ixghi^d;
    ^d&ixomin^d=ixmlo^d;
    ^d&ixomax^d=ixmhi^d;
 
    do ifl=1,numfl
      ilp=1
      tcnow=tlcoff(ifl)
      do while (ilp<=numr(ifl))      
        ! find points in the same grid
        ipmin=ilp
        do while (ipel(ifl,ipmin)/=mype .and. ipmin<=numr(ifl))
          ipmin=ipmin+1
        enddo
        igrid=igridl(ifl,ipmin)
        ipmax=ipmin
        do while (ipel(ifl,ipmax)==mype .and. igridl(ifl,ipmax+1)==igrid .and. ipmax<numr(ifl))
          ipmax=ipmax+1
        enddo
 
        ! parameters for the grid
        ^d&dxb^d=rnode(rpdx^d_,igrid);
        ^d&dxmax^d=4*dxb^d;
 
        do ilp=ipmin,ipmax
          xfnow(:)=xf(ifl,ilp,:)
          ^d&ixbmin^d=floor((xfnow(^d)-dxmax^d-ps(igrid)%x(ixomin^dd,^d))/dxb^d)+ixomin^d;
          ^d&ixbmax^d=floor((xfnow(^d)+dxmax^d-ps(igrid)%x(ixomin^dd,^d))/dxb^d)+ixomin^d;
 
          ! do interpolation inside the grid
          {ixcmin^d=max(ixbmin^d,ixomin^d)\}
          {ixcmax^d=min(ixbmax^d,ixomax^d)\}
          xbmin^nd=rnode(rpxmin^nd_,igrid)
          xbmax^nd=rnode(rpxmax^nd_,igrid)
          ixmax^nd=floor((phys_trac_mask-xbmin^nd)/dxb^nd)+ixomin^nd
          if (xbmax^nd>phys_trac_mask) ixcmax^nd=min(ixmax^nd,ixcmax^nd)
          {do ixc^d=ixcmin^d,ixcmax^d\}
            ds=0.d0
            {ds=ds+(xfnow(^d)-ps(igrid)%x(ixc^dd,^d))**2\}
            ds=sqrt(ds)
            if(ds<1.0d-2*dxb1) then
              weight=(1/(1.0d-2*dxb1))**weightindex
            else
              weight=(1/ds)**weightindex
            endif
            ps(igrid)%wextra(ixc^d,tweight_)=ps(igrid)%wextra(ixc^d,tweight_)+weight
            ps(igrid)%wextra(ixc^d,tcoff_)=ps(igrid)%wextra(ixc^d,tcoff_)+weight*tcnow
          {enddo\}
          
          ! do interpolation in neighbor grids that have the same level and pe
          {
            if (ixbmin^d<ixomin^d) then
              idn^dd=0;
              idn^d=-1
              if (neighbor(2,idn^dd,igrid)==mype .and. neighbor_type(idn^dd,igrid)==neighbor_sibling) then
                igrid_nb=neighbor(1,idn^dd,igrid)
                ixcmin^dd=max(ixbmin^dd,ixomin^dd);
                ixcmax^dd=min(ixbmax^dd,ixomax^dd);
                ixcmin^d=ixomax^d+(ixbmin^d-ixomin^d)
                ixcmax^d=ixomax^d
                xbmin^nd=rnode(rpxmin^nd_,igrid_nb)
                xbmax^nd=rnode(rpxmax^nd_,igrid_nb)
                ixmax^nd=floor((phys_trac_mask-xbmin^nd)/dxb^nd)+ixomin^nd
                if (xbmax^nd>phys_trac_mask) ixcmax^nd=min(ixmax^nd,ixcmax^nd)
 
                {do ixc^dd=ixcmin^dd,ixcmax^dd;}
                  ds=0.d0
                  {ds=ds+(xfnow(^dd)-ps(igrid_nb)%x({ixc^dd},^dd))**2;}
                  ds=sqrt(ds)
                  if(ds<1.0d-2*dxb1) then
                    weight=(1/(1.0d-2*dxb1))**weightindex
                  else
                    weight=(1/ds)**weightindex
                  endif
                  ps(igrid_nb)%wextra(ixc^dd,tweight_)=ps(igrid_nb)%wextra(ixc^dd,tweight_)+weight
                  ps(igrid_nb)%wextra(ixc^dd,tcoff_)=ps(igrid_nb)%wextra(ixc^dd,tcoff_)+weight*tcnow
                {enddo;}
              endif
            endif
 
            if (ixbmax^d>ixomin^d) then
              idn^dd=0;
              idn^d=1
              if (neighbor(2,idn^dd,igrid)==mype .and. neighbor_type(idn^dd,igrid)==neighbor_sibling) then
                igrid_nb=neighbor(1,idn^dd,igrid)
                xbmin^nd=rnode(rpxmin^nd_,igrid_nb)
                if (xbmin^nd<phys_trac_mask) then
                  ixcmin^dd=max(ixbmin^dd,ixomin^dd);
                  ixcmax^dd=min(ixbmax^dd,ixomax^dd);
                  ixcmin^d=ixomin^d
                  ixcmax^d=ixomin^d+(ixbmax^d-ixomax^d)
                  xbmax^nd=rnode(rpxmax^nd_,igrid_nb)
                  ixmax^nd=floor((phys_trac_mask-xbmin^nd)/dxb^nd)+ixomin^nd
                  if (xbmax^nd>phys_trac_mask) ixcmax^nd=min(ixmax^nd,ixcmax^nd)
 
                  {do ixc^dd=ixcmin^dd,ixcmax^dd;}
                    ds=0.d0
                    {ds=ds+(xfnow(^dd)-ps(igrid_nb)%x({ixc^dd},^dd))**2;}
                    ds=sqrt(ds)
                    if(ds<1.0d-2*dxb1) then
                      weight=(1/(1.0d-2*dxb1))**weightindex
                    else
                      weight=(1/ds)**weightindex
                    endif
                    ps(igrid_nb)%wextra(ixc^dd,tweight_)=ps(igrid_nb)%wextra(ixc^dd,tweight_)+weight
                    ps(igrid_nb)%wextra(ixc^dd,tcoff_)=ps(igrid_nb)%wextra(ixc^dd,tcoff_)+weight*tcnow
                  {enddo;}
                endif
              endif
            endif
          \}
        enddo
 
        ilp=ipmax+1
      enddo
    enddo
 
 
    do j=1,ngrid_trac
      igrid=trac_grid(j)
      where(ps(igrid)%wextra(ixo^s,tweight_)>0.d0)
        ps(igrid)%wextra(ixo^s,tcoff_)=ps(igrid)%wextra(ixo^s,tcoff_)/ps(igrid)%wextra(ixo^s,tweight_)
      elsewhere
        ps(igrid)%wextra(ixo^s,tcoff_)=t_bott
      endwhere
    enddo
 
  end subroutine interp_tcoff
 
 
  subroutine trac_after_setdt(tco,trac_alfa,T_peak, T_bott_in)
       double precision, intent(in)    :: trac_alfa,tco,T_peak, T_bott_in
      ! default lower limit of cutoff temperature
      t_bott = t_bott_in
      select case(phys_trac_type)
      case(0)
        !> Test case, fixed cutoff temperature
        !> do nothing here
      case(1)
        !> 1D TRAC method
        call trac_simple(tco,trac_alfa,t_peak)
      case(2)
        !> LTRAC method, by iijima et al. 2021
        !> do nothing here 
        call ltrac(t_peak)
      case(3)
        !> 2D or 3D TRACL(ine) method
        call tracl(.false.,t_peak)
      case(4)
        !> 2D or 3D TRACB(lock) method
        call tracb(.false.,t_peak)
      case(5)
        !> 2D or 3D TRACL(ine) method with mask
        call tracl(.true.,t_peak)
      case(6)
        !> 2D or 3D TRACB(lock) method with mask
        call tracb(.true.,t_peak)
      case default
        call mpistop("undefined TRAC method type")
      end select
 
  end subroutine trac_after_setdt
 
  subroutine initialize_trac_after_settree
    use mod_global_parameters
    use mod_physics, only: phys_trac_after_setdt
  
    if(phys_trac) then
      phys_trac_after_setdt => trac_after_setdt
      if(phys_trac_type .eq. 3) then
        if(mype .eq. 0) write(*,*) 'Using TRACL(ine) global method'
        if(mype .eq. 0) write(*,*) 'By default, magnetic field lines are traced every 4 grid cells'
        call init_trac_line(.false.)
      end if
      if(phys_trac_type .eq. 4) then
        if(mype .eq. 0) write(*,*) 'Using TRACB(lock) global method'
        if(mype .eq. 0) write(*,*) 'Currently, only valid in Cartesian uniform settings'
        if(mype .eq. 0) write(*,*) 'By default, magnetic field lines are traced every 4 grid cells'
        call init_trac_block(.false.)
      end if
      if(phys_trac_type .eq. 5) then
        if(mype .eq. 0) write(*,*) 'Using TRACL(ine) method with a mask'
        if(mype .eq. 0) write(*,*) 'By default, magnetic field lines are traced every 4 grid cells'
        call init_trac_line(.true.)
      end if
      if(phys_trac_type .eq. 6) then
        if(mype .eq. 0) write(*,*) 'Using TRACB(lock) method with a mask'
        if(mype .eq. 0) write(*,*) 'Currently, only valid in Cartesian uniform settings'
        if(mype .eq. 0) write(*,*) 'By default, magnetic field lines are traced every 4 grid cells'
        call init_trac_block(.true.)
      end if
    end if
  
  end subroutine initialize_trac_after_settree
 
 
 
end module mod_trac