mod_fld.t

Go to the documentation of this file.

!> Nicolas Moens
!> Module for including flux limited diffusion (FLD)-approximation in Radiation-hydrodynamics simulations using mod_rhd
!> Based on Turner and stone 2001. See
!> [1]Moens, N., Sundqvist, J. O., El Mellah, I., Poniatowski, L., Teunissen, J., and Keppens, R.,
!> “Radiation-hydrodynamics with MPI-AMRVAC . Flux-limited diffusion”,
!> <i>Astronomy and Astrophysics</i>, vol. 657, 2022. doi:10.1051/0004-6361/202141023.
!> For more information.
 
module mod_fld
    implicit none
 
    !> source split for energy interact and radforce:
    logical :: fld_eint_split = .false.
    logical :: fld_radforce_split = .false.
 
    !> Opacity per unit of unit_density
    double precision, public :: fld_kappa0 = 0.34d0
 
    !> mean particle mass
    double precision, public :: fld_mu = 0.6d0
 
    !> Tolerance for bisection method for Energy sourceterms
    !> This is a percentage of the minimum of gas- and radiation energy
    double precision, public :: fld_bisect_tol = 1.d-4
 
    !> Tolerance for adi method for radiative Energy diffusion
    double precision, public :: fld_diff_tol = 1.d-4
 
    !> Number for splitting the diffusion module
    double precision, public :: diff_crit
 
    !> Use constant Opacity?
    character(len=8) :: fld_opacity_law = 'const'
    character(len=6) :: fld_opal_table = 'Y09800' !>'xxxxxx'
 
    !> Diffusion limit lambda = 0.33
    character(len=16) :: fld_fluxlimiter = 'Pomraning'
 
    !> diffusion coefficient for multigrid method
    integer :: i_diff_mg
 
    !> Which method to solve diffusion part
    character(len=8) :: fld_diff_scheme = 'mg'
 
    !> Which method to find the root for the energy interaction polynomial
    character(len=8) :: fld_interaction_method = 'Halley'
 
    !> Take running average for Diffusion coefficient
    logical :: diff_coef_filter = .false.
    integer :: size_d_filter = 1
 
    !> Take a running average over the fluxlimiter
    logical :: flux_lim_filter = .false.
    integer :: size_l_filter = 1
 
    !> Use or dont use lineforce opacities
    logical :: lineforce_opacities = .false.
 
    !> Resume run when multigrid returns error
    logical :: diffcrash_resume = .true.
 
    !> Index for Flux weighted opacities
    integer, allocatable, public :: i_opf(:)
 
    !> A copy of rhd_Gamma
    double precision, private, protected :: fld_gamma
 
    !> running timestep for diffusion solver, initialised as zero
    double precision :: dt_diff = 0.d0
 
    !> public methods
    !> these are called in mod_rhd_phys
    public :: get_fld_rad_force
    public :: get_fld_energy_interact
    public :: fld_radforce_get_dt
    public :: fld_init
    public :: fld_get_radflux
    public :: fld_get_radpress
    public :: fld_get_fluxlimiter
    public :: fld_get_opacity
 
  contains
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!! GENERAL
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  !> Reading in fld-list parameters from .par file
  subroutine fld_params_read(files)
    use mod_global_parameters, only: unitpar
    use mod_constants
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /fld_list/ fld_kappa0, fld_eint_split, fld_radforce_split, &
    fld_bisect_tol, fld_diff_tol,&
    fld_opacity_law, fld_fluxlimiter, fld_diff_scheme, fld_interaction_method, &
    diff_coef_filter, size_d_filter, flux_lim_filter, size_l_filter, &
    lineforce_opacities, diffcrash_resume, fld_opal_table
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status="old")
       read(unitpar, fld_list, end=111)
       111    close(unitpar)
    end do
 
  end subroutine fld_params_read
 
  !> Initialising FLD-module:
  !> Read opacities
  !> Initialise Multigrid
  !> adimensionalise kappa
  !> Add extra variables to w-array, flux, kappa, eddington Tensor
  !> Lambda and R
  !> ...
  subroutine fld_init(He_abundance, rhd_radiation_diffusion, rhd_gamma)
    use mod_global_parameters
    use mod_variables
    use mod_physics
    use mod_opal_opacity, only: init_opal
    use mod_multigrid_coupling
 
    double precision, intent(in) :: he_abundance, rhd_gamma
    logical, intent(in) :: rhd_radiation_diffusion
    double precision :: sigma_thomson
    integer :: idir,jdir
 
    character(len=1) :: ind_1
    character(len=1) :: ind_2
    character(len=2) :: cmp_f
    character(len=5) :: cmp_e
 
    !> read par files
    call fld_params_read(par_files)
 
    !> Set lineforce opacities as variable
    if (lineforce_opacities) then
      allocate(i_opf(ndim))
      do idir = 1,ndim
        write(ind_1,'(I1)') idir
        cmp_f = 'k' // ind_1
        i_opf(idir) = var_set_extravar(cmp_f,cmp_f)
      enddo
    endif
 
    if (rhd_radiation_diffusion) then
      if (fld_diff_scheme .eq. 'mg') then
 
        use_multigrid = .true.
 
        phys_implicit_update => diffuse_e_rad_mg
        phys_evaluate_implicit => evaluate_e_rad_mg
 
        mg%n_extra_vars = 1
        mg%operator_type = mg_vhelmholtz
 
      endif
    endif
 
    i_diff_mg = var_set_extravar("D", "D")
 
    !> Need mean molecular weight
    fld_mu = (1.+4*he_abundance)/(2.+3.*he_abundance)
 
    !> set rhd_gamma
    fld_gamma = rhd_gamma
 
    !> Read in opacity table if necesary
    if (fld_opacity_law .eq. 'opal') call init_opal(he_abundance,fld_opal_table)
    if ((fld_opacity_law .eq. 'thomson') .or. (fld_opacity_law .eq. 'fastwind'))  then
      sigma_thomson = 6.6524585d-25
      fld_kappa0 = sigma_thomson/const_mp * (1.+2.*he_abundance)/(1.+4.*he_abundance)
    endif
  end subroutine fld_init
 
  !> w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
  !> This subroutine handles the radiation force
  subroutine get_fld_rad_force(qdt,ixI^L,ixO^L,wCT,w,x,&
       energy,qsourcesplit,active)
    use mod_constants
    use mod_global_parameters
    use mod_usr_methods
    use mod_geometry
 
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: qdt, x(ixi^s,1:ndim)
    double precision, intent(in)    :: wct(ixi^s,1:nw)
    double precision, intent(inout) :: w(ixi^s,1:nw)
    logical, intent(in) :: energy,qsourcesplit
    logical, intent(inout) :: active
 
    double precision :: radiation_forcect(ixo^s,1:ndim)
    double precision :: kappact(ixo^s)
    double precision :: rad_fluxct(ixo^s,1:ndim)
 
    double precision :: div_v(ixi^s,1:ndim,1:ndim)
    double precision :: edd(ixo^s,1:ndim,1:ndim)
    double precision :: nabla_vp(ixo^s)
    double precision :: vel(ixi^s), grad_v(ixi^s), grad0_v(ixo^s)
    double precision :: grad_e(ixo^s)
 
    integer :: idir, jdir
 
    double precision :: fld_r(ixo^s), lambda(ixo^s)
 
    !> Calculate and add sourceterms
    if(qsourcesplit .eqv. fld_radforce_split) then
      active = .true.
 
      call fld_get_opacity(wct, x, ixi^l, ixo^l, kappact)
      call fld_get_radflux(wct, x, ixi^l, ixo^l, rad_fluxct)
      call fld_get_fluxlimiter(wct, x, ixi^l, ixo^l, lambda, fld_r)
 
      do idir = 1,ndim
        !> Radiation force = kappa*rho/c *Flux = lambda gradE
        radiation_forcect(ixo^s,idir) = kappact(ixo^s)*rad_fluxct(ixo^s,idir)/(const_c/unit_velocity)
 
        ! call gradientO(wCT(ixI^S,iw_r_e),x,ixI^L,ixO^L,idir,grad_E,nghostcells)
        ! radiation_forceCT(ixO^S,idir) = lambda(ixO^S)*grad_E(ixO^S)
 
        !> Momentum equation source term
        w(ixo^s,iw_mom(idir)) = w(ixo^s,iw_mom(idir)) &
            + qdt * wct(ixo^s,iw_rho)*radiation_forcect(ixo^s,idir)
        ! w(ixO^S,iw_mom(idir)) = w(ixO^S,iw_mom(idir)) &
        !     + qdt *radiation_forceCT(ixO^S,idir)
 
        ! if (energy .and. .not. block%e_is_internal) then
          !> Energy equation source term (kinetic energy)
          w(ixo^s,iw_e) = w(ixo^s,iw_e) &
              + qdt * wct(ixo^s,iw_mom(idir))*radiation_forcect(ixo^s,idir)
          ! w(ixO^S,iw_e) = w(ixO^S,iw_e) &
          !     + qdt * wCT(ixO^S,iw_mom(idir))/wCT(ixO^S,iw_rho)*radiation_forceCT(ixO^S,idir)
        ! endif
      enddo
 
      !> Photon tiring
      !> calculate tensor div_v
      !> !$OMP PARALLEL DO
      do idir = 1,ndim
        do jdir = 1,ndim
          vel(ixi^s) = wct(ixi^s,iw_mom(jdir))/wct(ixi^s,iw_rho)
 
          call gradient(vel,ixi^l,ixo^l,idir,grad_v)
          div_v(ixo^s,idir,jdir) = grad_v(ixo^s)
 
          ! call gradientO(vel,x,ixI^L,ixO^L,idir,grad0_v,nghostcells)
          ! div_v(ixO^S,idir,jdir) = grad0_v(ixO^S)
        enddo
      enddo
      !> !$OMP END PARALLEL DO
 
      call fld_get_eddington(wct, x, ixi^l, ixo^l, edd)
 
      !> VARIABLE NAMES DIV ARE ACTUALLY GRADIENTS
      {^ifoned
      nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1)  &
      }
 
      {^iftwod
      !>eq 34 Turner and stone (Only 2D)
      nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1)  &
                   + div_v(ixo^s,1,2)*edd(ixo^s,1,2)  &
                   + div_v(ixo^s,2,1)*edd(ixo^s,2,1)  &
                   + div_v(ixo^s,2,2)*edd(ixo^s,2,2)
      }
 
      {^ifthreed
      nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1)  &
                   + div_v(ixo^s,1,2)*edd(ixo^s,1,2)  &
                   + div_v(ixo^s,1,3)*edd(ixo^s,1,3)  &
                   + div_v(ixo^s,2,1)*edd(ixo^s,2,1)  &
                   + div_v(ixo^s,2,2)*edd(ixo^s,2,2)  &
                   + div_v(ixo^s,2,3)*edd(ixo^s,2,3)  &
                   + div_v(ixo^s,3,1)*edd(ixo^s,3,1)  &
                   + div_v(ixo^s,3,2)*edd(ixo^s,3,2)  &
                   + div_v(ixo^s,3,3)*edd(ixo^s,3,3)
      }
 
      w(ixo^s,iw_r_e) = w(ixo^s,iw_r_e) &
          - qdt * nabla_vp(ixo^s)*wct(ixo^s,iw_r_e)
 
    end if
 
  end subroutine get_fld_rad_force
 
 
  subroutine fld_radforce_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
    use mod_global_parameters
    use mod_usr_methods
 
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: dx^d, x(ixi^s,1:ndim), w(ixi^s,1:nw)
    double precision, intent(inout) :: dtnew
 
    double precision :: radiation_force(ixo^s,1:ndim)
    double precision :: rad_flux(ixo^s,1:ndim)
    double precision :: kappa(ixo^s)
    double precision                :: dxinv(1:ndim), max_grad
    integer                         :: idim
 
    ^d&dxinv(^d)=one/dx^d;
 
    call fld_get_opacity(w, x, ixi^l, ixo^l, kappa)
    call fld_get_radflux(w, x, ixi^l, ixo^l, rad_flux)
 
    do idim = 1, ndim
      radiation_force(ixo^s,idim) = w(ixo^s,iw_rho)*kappa(ixo^s)*rad_flux(ixo^s,idim)/(const_c/unit_velocity)
      max_grad = maxval(abs(radiation_force(ixo^s,idim)))
      max_grad = max(max_grad, epsilon(1.0d0))
      dtnew = min(dtnew, courantpar / sqrt(max_grad * dxinv(idim)))
    end do
 
  end subroutine fld_radforce_get_dt
 
  !> w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
  !> This subroutine handles the energy exchange between gas and radiation
  subroutine get_fld_energy_interact(qdt,ixI^L,ixO^L,wCT,w,x,&
       energy,qsourcesplit,active)
    use mod_constants
    use mod_global_parameters
    use mod_usr_methods
 
    use mod_physics, only: phys_get_pthermal  !needed to get temp
 
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: qdt, x(ixi^s,1:ndim)
    double precision, intent(in)    :: wct(ixi^s,1:nw)
    double precision, intent(inout) :: w(ixi^s,1:nw)
    logical, intent(in) :: energy,qsourcesplit
    logical, intent(inout) :: active
 
    !> Calculate and add sourceterms
    if(qsourcesplit .eqv. fld_eint_split) then
      active = .true.
      !> Add energy sourceterms
      call energy_interaction(w, w, x, ixi^l, ixo^l)
    end if
  end subroutine get_fld_energy_interact
 
  !> Sets the opacity in the w-array
  !> by calling mod_opal_opacity
  subroutine fld_get_opacity(w, x, ixI^L, ixO^L, fld_kappa)
    use mod_global_parameters
    use mod_physics, only: phys_get_pthermal
    use mod_physics, only: phys_get_tgas
    use mod_usr_methods
    use mod_opal_opacity
 
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out) :: fld_kappa(ixo^s)
    double precision :: temp(ixi^s), pth(ixi^s), a2(ixo^s)
    double precision :: rho0,temp0,n,sigma_b
    double precision :: akram, bkram
    double precision :: vth(ixo^s), gradv(ixi^s), eta(ixo^s), t(ixo^s)
 
    integer :: i,j,ix^d, idir
 
    select case (fld_opacity_law)
      case('const')
        fld_kappa(ixo^s) = fld_kappa0/unit_opacity
      case('thomson')
        fld_kappa(ixo^s) = fld_kappa0/unit_opacity
      case('kramers')
        rho0 = half !> Take lower value of rho in domain
        fld_kappa(ixo^s) = fld_kappa0/unit_opacity*((w(ixo^s,iw_rho)/rho0))
      case('bump')
        !> Opacity bump
        rho0 = 0.2d0 !0.5d-1
        n = 7.d0
        sigma_b = 2.d-2
        !fld_kappa(ixO^S) = fld_kappa0/unit_opacity*(one + n*dexp(-((rho0  - w(ixO^S,iw_rho))**two)/rho0))
        fld_kappa(ixo^s) = fld_kappa0/unit_opacity*(one + n*dexp(-one/sigma_b*(dlog(w(ixo^s,iw_rho)/rho0))**two))
      case('non_iso')
        call phys_get_pthermal(w,x,ixi^l,ixo^l,temp)
        temp(ixo^s)=temp(ixo^s)/w(ixo^s,iw_rho)
 
        rho0 = 0.5d0 !> Take lower value of rho in domain
        temp0 = one
        n = -7.d0/two
        fld_kappa(ixo^s) = fld_kappa0/unit_opacity*(w(ixo^s,iw_rho)/rho0)*(temp(ixo^s)/temp0)**n
      case('fastwind')
        call phys_get_pthermal(w,x,ixi^l,ixo^l,pth)
        a2(ixo^s) = pth(ixo^s)/w(ixo^s,iw_rho)*unit_velocity**2.d0
 
        akram = 13.1351597305
        bkram = -4.5182188206
 
        fld_kappa(ixo^s) = fld_kappa0/unit_opacity &
        * (1.d0+10.d0**akram*w(ixo^s,iw_rho)*unit_density*(a2(ixo^s)/1.d12)**bkram)
 
        {do ix^d=ixomin^d,ixomax^d\ }
          !> Hard limit on kappa
          fld_kappa(ix^d) = min(fld_kappa(ix^d),2.3d0*fld_kappa0/unit_opacity)
 
          !> Limit kappa through T
          ! fld_kappa(ix^D) = fld_kappa0/unit_opacity &
          ! * (1.d0+10.d0**akram*w(ix^D,iw_rho)*unit_density &
          ! * (max(a2(ix^D),const_kB*5.9d4/(fld_mu*const_mp))/1.d12)**bkram)
        {enddo\ }
 
      case('opal')
        call phys_get_tgas(w,x,ixi^l,ixo^l,temp)
        {do ix^d=ixomin^d,ixomax^d\ }
            rho0 = w(ix^d,iw_rho)*unit_density
            temp0 = temp(ix^d)*unit_temperature
            call set_opal_opacity(rho0,temp0,n)
            fld_kappa(ix^d) = n/unit_opacity
        {enddo\ }
 
      case('special')
        if (.not. associated(usr_special_opacity)) then
          call mpistop("special opacity not defined")
        endif
        call usr_special_opacity(ixi^l, ixo^l, w, x, fld_kappa)
 
      case default
        call mpistop("Doesn't know opacity law")
      end select
  end subroutine fld_get_opacity
 
  !> Calculate fld flux limiter
  !> This subroutine calculates flux limiter lambda using the prescription
  !> stored in fld_fluxlimiter.
  !> It also calculates the ratio of radiation scaleheight and mean free path
  subroutine fld_get_fluxlimiter(w, x, ixI^L, ixO^L, fld_lambda, fld_R)
    use mod_global_parameters
    use mod_geometry
    use mod_usr_methods
 
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out) :: fld_r(ixo^s), fld_lambda(ixo^s)
    double precision :: kappa(ixo^s)
    double precision ::  normgrad2(ixo^s)
    double precision :: grad_r_e(ixi^s), grad_r_eo(ixo^s)
    double precision :: rad_e(ixi^s)
    integer :: idir, ix^d
 
    double precision :: tmp_l(ixi^s), filtered_l(ixi^s)
    integer :: filter, idim
 
    select case (fld_fluxlimiter)
    case('Diffusion')
      fld_lambda(ixo^s) = 1.d0/3.d0
      fld_r(ixo^s) = zero
 
    case('FreeStream')
      !> Calculate R everywhere
      !> |grad E|/(rho kappa E)
      normgrad2(ixo^s) = 0.d0 !smalldouble
 
      rad_e(ixi^s) = w(ixi^s, iw_r_e)
      do idir = 1,ndim
        call gradiento(rad_e,x,ixi^l,ixo^l,idir,grad_r_eo,nghostcells)
        normgrad2(ixo^s) = normgrad2(ixo^s) + grad_r_eo(ixo^s)**2
 
        ! call gradient(rad_e,ixI^L,ixO^L,idir,grad_r_e)
        ! normgrad2(ixO^S) = normgrad2(ixO^S) + grad_r_e(ixO^S)**2
      end do
 
      call fld_get_opacity(w, x, ixi^l, ixo^l, kappa)
 
      fld_r(ixo^s) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
 
      !> Calculate the flux limiter, lambda
      fld_lambda(ixo^s) = one/fld_r(ixo^s)
 
    case('Pomraning')
      !> Calculate R everywhere
      !> |grad E|/(rho kappa E)
      normgrad2(ixo^s) = 0.d0 !smalldouble !*w(ixO^S,iw_r_e)
 
      rad_e(ixi^s) = w(ixi^s, iw_r_e)
      do idir = 1,ndim
        call gradiento(rad_e,x,ixi^l,ixo^l,idir,grad_r_eo,nghostcells)
        normgrad2(ixo^s) = normgrad2(ixo^s) + grad_r_eo(ixo^s)**2
 
        ! call gradient(rad_e,ixI^L,ixO^L,idir,grad_r_e)
        ! normgrad2(ixO^S) = normgrad2(ixO^S) + grad_r_e(ixO^S)**2
      end do
 
      call fld_get_opacity(w, x, ixi^l, ixo^l, kappa)
 
      fld_r(ixo^s) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
 
      !> Calculate the flux limiter, lambda
      !> Levermore and Pomraning: lambda = (2 + R)/(6 + 3R + R^2)
      fld_lambda(ixo^s) = (2.d0+fld_r(ixo^s))/(6.d0+3*fld_r(ixo^s)+fld_r(ixo^s)**2.d0)
 
    case('Pomraning2')
      call mpistop("Pomraning2 is not quite working, use Pomraning or Minerbo")
      !> Calculate R everywhere
      !> |grad E|/(rho kappa E)
      normgrad2(ixo^s) = 0.d0 !smalldouble
 
      rad_e(ixi^s) = w(ixi^s, iw_r_e)
      do idir = 1,ndim
        call gradiento(rad_e,x,ixi^l,ixo^l,idir,grad_r_eo,nghostcells)
        normgrad2(ixo^s) = normgrad2(ixo^s) + grad_r_eo(ixo^s)**2
 
        ! call gradient(rad_e,ixI^L,ixO^L,idir,grad_r_e)
        ! normgrad2(ixO^S) = normgrad2(ixO^S) + grad_r_e(ixO^S)**2
      end do
 
      call fld_get_opacity(w, x, ixi^l, ixo^l, kappa)
 
      fld_r(ixo^s) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
 
      !> Calculate the flux limiter, lambda
      !> Levermore and Pomraning: lambda = 1/R(coth(R)-1/R)
      fld_lambda(ixo^s) = one/fld_r(ixo^s)*(one/dtanh(fld_r(ixo^s)) - one/fld_r(ixo^s))
      ! fld_lambda(ixO^S) = one/fld_R(ixO^S)*(dcosh(fld_R(ixO^S))/dsinh(fld_R(ixO^S)) - one/fld_R(ixO^S))
 
      !>WHAT HAPPENS WHEN R=0 (full diffusion) => 1/R = NAN => dtanh(1/R) =????
      where(dtanh(fld_r(ixo^s)) .ne. dtanh(fld_r(ixo^s)))
        fld_lambda(ixo^s) = 1.d0/3.d0
      endwhere
 
    case('Minerbo')
      !> Calculate R everywhere
      !> |grad E|/(rho kappa E)
      normgrad2(ixo^s) = 0.d0 !smalldouble
 
      rad_e(ixi^s) = w(ixi^s, iw_r_e)
      do idir = 1,ndim
        call gradiento(rad_e,x,ixi^l,ixo^l,idir,grad_r_eo,nghostcells)
        normgrad2(ixo^s) = normgrad2(ixo^s) + grad_r_eo(ixo^s)**2
 
        ! call gradient(rad_e,ixI^L,ixO^L,idir,grad_r_e)
        ! normgrad2(ixO^S) = normgrad2(ixO^S) + grad_r_e(ixO^S)**2
      end do
 
      call fld_get_opacity(w, x, ixi^l, ixo^l, kappa)
 
      fld_r(ixo^s) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
 
      !> Calculate the flux limiter, lambda
      !> Minerbo:
      {do ix^d=ixomin^d,ixomax^d\ }
          if (fld_r(ix^d) .lt. 3.d0/2.d0) then
            fld_lambda(ix^d) = 2.d0/(3.d0 + dsqrt(9.d0 + 12.d0*fld_r(ix^d)**2.d0))
          else
            fld_lambda(ix^d) = 1.d0/(1.d0 + fld_r(ix^d) + dsqrt(1.d0 + 2.d0*fld_r(ix^d)))
          endif
      {enddo\}
 
    case('special')
      if (.not. associated(usr_special_fluxlimiter)) then
        call mpistop("special fluxlimiter not defined")
      endif
      call usr_special_fluxlimiter(ixi^l, ixo^l, w, x, fld_lambda, fld_r)
    case default
      call mpistop('Fluxlimiter unknown')
    end select
 
 
    if (flux_lim_filter) then
      if (size_l_filter .lt. 1) call mpistop("D filter of size < 1 makes no sense")
      if (size_l_filter .gt. nghostcells) call mpistop("D filter of size > nghostcells makes no sense")
 
      tmp_l(ixo^s) = fld_lambda(ixo^s)
      filtered_l(ixo^s) = zero
 
      do filter = 1,size_l_filter
        {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_l_filter\}
          do idim = 1,ndim
            filtered_l(ix^d) = filtered_l(ix^d) &
                             + tmp_l(ix^d+filter*kr(idim,^d)) &
                             + tmp_l(ix^d-filter*kr(idim,^d))
          enddo
        {enddo\}
      enddo
 
      {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_d_filter\}
        tmp_l(ix^d) = (tmp_l(ix^d)+filtered_l(ix^d))/(1+2*size_l_filter*ndim)
      {enddo\}
 
      fld_lambda(ixo^s) = tmp_l(ixo^s)
    endif
 
  end subroutine fld_get_fluxlimiter
 
  !> Calculate Radiation Flux
  !> stores radiation flux in w-array
  subroutine fld_get_radflux(w, x, ixI^L, ixO^L, rad_flux)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out) :: rad_flux(ixo^s, 1:ndim)
 
    double precision :: grad_r_e(ixi^s), grad_r_eo(ixo^s)
    double precision :: rad_e(ixi^s)
    double precision :: kappa(ixo^s), lambda(ixo^s), fld_r(ixo^s)
    integer :: ix^d, idir
 
    rad_e(ixi^s) = w(ixi^s, iw_r_e)
 
    call fld_get_opacity(w, x, ixi^l, ixo^l, kappa)
    call fld_get_fluxlimiter(w, x, ixi^l, ixo^l, lambda, fld_r)
 
    !> Calculate the Flux using the fld closure relation
    !> F = -c*lambda/(kappa*rho) *grad E
    do idir = 1,ndim
      call gradiento(rad_e,x,ixi^l,ixo^l,idir,grad_r_eo,nghostcells)
      rad_flux(ixo^s, idir) = -(const_c/unit_velocity)*lambda(ixo^s)/(kappa(ixo^s)*w(ixo^s,iw_rho))*grad_r_eo(ixo^s)
 
      ! call gradient(rad_e,ixI^L,ixO^L,idir,grad_r_e)
      ! rad_flux(ixO^S, idir) = -(const_c/unit_velocity)*lambda(ixO^S)/(kappa(ixO^S)*w(ixO^S,iw_rho))*grad_r_e(ixO^S)
    end do
 
  end subroutine fld_get_radflux
 
  !> Calculate Eddington-tensor
  !> Stores Eddington-tensor in w-array
  subroutine fld_get_eddington(w, x, ixI^L, ixO^L, eddington_tensor)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out) :: eddington_tensor(ixO^S,1:ndim,1:ndim)
    double precision :: tnsr2(ixO^S,1:ndim,1:ndim)
    double precision :: normgrad2(ixO^S), f(ixO^S)
    double precision :: grad_r_e(ixI^S, 1:ndim), rad_e(ixI^S)
    double precision :: grad_r_eO(ixO^S, 1:ndim)
    double precision :: lambda(ixO^S), fld_R(ixO^S)
    integer :: i,j, idir,jdir
 
    !> Calculate R everywhere
    !> |grad E|/(rho kappa E)
    normgrad2(ixo^s) = zero
 
    rad_e(ixi^s) = w(ixi^s, iw_r_e)
    do idir = 1,ndim
      call gradiento(rad_e,x,ixi^l,ixo^l,idir,grad_r_eo(ixo^s, idir),nghostcells)
      normgrad2(ixo^s) = normgrad2(ixo^s) + grad_r_eo(ixo^s,idir)**2
 
      ! call gradient(rad_e,ixI^L,ixO^L,idir,grad_r_e(ixI^S,idir))
      ! normgrad2(ixO^S) = normgrad2(ixO^S) + grad_r_e(ixO^S,idir)**two
    end do
 
    call fld_get_fluxlimiter(w, x, ixi^l, ixo^l, lambda, fld_r)
 
    !> Calculate radiation pressure
    !> P = (lambda + lambda^2 R^2)*E
    f(ixo^s) = lambda(ixo^s) + lambda(ixo^s)**two * fld_r(ixo^s)**two
    f(ixo^s) = one/two*(one-f(ixo^s)) + one/two*(3.d0*f(ixo^s) - one)
 
    {^ifoned
    eddington_tensor(ixo^s,1,1) = f(ixo^s)
    }
 
    {^nooned
    do idir = 1,ndim
      eddington_tensor(ixo^s,idir,idir) = half*(one-f(ixo^s))
    enddo
 
    do idir = 1,ndim
      do jdir = 1,ndim
        if (idir .ne. jdir) eddington_tensor(ixo^s,idir,jdir) = zero
        tnsr2(ixo^s,idir,jdir) =  half*(3.d0*f(ixo^s) - 1)&
        *grad_r_eo(ixo^s,idir)*grad_r_eo(ixo^s,jdir)/normgrad2(ixo^s)
      enddo
    enddo
 
    ! do idir = 1,ndim
    !   do jdir = 1,ndim
    !     if (idir .ne. jdir) eddington_tensor(ixO^S,idir,jdir) = zero
    !     tnsr2(ixO^S,idir,jdir) =  half*(3.d0*f(ixO^S) - 1)&
    !     *grad_r_e(ixO^S,idir)*grad_r_e(ixO^S,jdir)/normgrad2(ixO^S)
    !   enddo
    ! enddo
 
    do idir = 1,ndim
      do jdir = 1,ndim
        where ((tnsr2(ixo^s,idir,jdir) .eq. tnsr2(ixo^s,idir,jdir)) &
          .and. (normgrad2(ixo^s) .gt. smalldouble))
          eddington_tensor(ixo^s,idir,jdir) = eddington_tensor(ixo^s,idir,jdir) + tnsr2(ixo^s,idir,jdir)
        endwhere
      enddo
    enddo
    }
 
  end subroutine fld_get_eddington
 
  !> Calculate Radiation Pressure
  !> Returns Radiation Pressure as tensor
  subroutine fld_get_radpress(w, x, ixI^L, ixO^L, rad_pressure)
    use mod_global_parameters
 
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision             :: eddington_tensor(ixo^s,1:ndim,1:ndim)
    double precision, intent(out):: rad_pressure(ixo^s,1:ndim,1:ndim)
 
    integer i,j
 
    call fld_get_eddington(w, x, ixi^l, ixo^l, eddington_tensor)
 
    do i=1,ndim
      do j=1,ndim
        rad_pressure(ixo^s,i,j) = eddington_tensor(ixo^s,i,j)* w(ixo^s,iw_r_e)
      enddo
    enddo
  end subroutine fld_get_radpress
 
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!!!!!!!!!!!!!! Multigrid diffusion
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  !> Calling all subroutines to perform the multigrid method
  !> Communicates rad_e and diff_coeff to multigrid library
  subroutine diffuse_e_rad_mg(dtfactor,qdt,qtC,psa,psb)
    use mod_global_parameters
    use mod_forest
    ! use mod_ghostcells_update
    use mod_multigrid_coupling
    ! use mod_physics, only: phys_set_mg_bounds, phys_req_diagonal
 
    type(state), target :: psa(max_blocks) !< Advance psa=psb+dtfactor*qdt*F_im(psa)
    type(state), target :: psb(max_blocks)
    double precision, intent(in) :: qdt
    double precision, intent(in) :: qtC
    double precision, intent(in) :: dtfactor
 
    integer                      :: n
    double precision             :: res, max_residual, lambda
    integer, parameter           :: max_its = 100
 
    integer                        :: iw_to,iw_from
    integer                        :: iigrid, igrid, id
    integer                        :: nc, lvl, i
    type(tree_node), pointer       :: pnode
    real(dp)                       :: fac, facD
 
    ! print*, '1'
 
    !> Set diffusion timestep, add previous timestep if mg did not converge:
    if (it == 0) dt_diff = 0
    dt_diff = dt_diff + qdt
 
    ! Avoid setting a very restrictive limit to the residual when the time step
    ! is small (as the operator is ~ 1/(D * qdt))
    if (qdt < dtmin) then
        if(mype==0)then
            print *,'skipping implicit solve: dt too small!'
            print *,'Currently at time=',global_time,' time step=',qdt,' dtmin=',dtmin
        endif
        return
    endif
    ! max_residual = 1d-7/qdt
    max_residual = fld_diff_tol !1d-7/qdt
 
    mg%operator_type = mg_vhelmholtz
    mg%smoother_type = mg_smoother_gs
    call mg_set_methods(mg)
 
    if (.not. mg%is_allocated) call mpistop("multigrid tree not allocated yet")
 
!   lambda = 1.d0/(dtfactor * qdt)
    lambda = 1.d0/(dtfactor * dt_diff)
    call vhelmholtz_set_lambda(lambda)
 
    call update_diffcoeff(psb)
 
    fac = 1.d0
    facd = 1.d0
 
    ! print*, '2'
 
    !This is mg_copy_to_tree from psb state
    call mg_copy_to_tree(i_diff_mg, mg_iveps, factor=facd, state_from=psb)
    !This is mg_copy_to_tree from psb state
    !!!  replaces::  call mg_copy_to_tree(su_, mg_irhs, factor=-lambda)
    ! iw_from=i_diff_mg
    ! iw_to=mg_iveps
    ! fac=facD
    ! do iigrid=1,igridstail; igrid=igrids(iigrid);
    !    pnode => igrid_to_node(igrid, mype)%node
    !    id    =  pnode%id
    !    lvl   =  mg%boxes(id)%lvl
    !    nc    =  mg%box_size_lvl(lvl)
    !    ! Include one layer of ghost cells on grid leaves
    !    {^IFONED
    !    mg%boxes(id)%cc(0:nc+1, iw_to) = fac * &
    !         psb(igrid)%w(ixMlo1-1:ixMhi1+1, iw_from)
    !    }
    !    {^IFTWOD
    !    mg%boxes(id)%cc(0:nc+1, 0:nc+1, iw_to) = fac * &
    !         psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, iw_from)
    !    }
    !    {^IFTHREED
    !    mg%boxes(id)%cc(0:nc+1, 0:nc+1, 0:nc+1, iw_to) = fac * &
    !         psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, &
    !         ixMlo3-1:ixMhi3+1, iw_from)
    !    }
    ! end do
 
    ! print*, '3'
 
    !This is mg_copy_to_tree from psb state
    call mg_copy_to_tree(iw_r_e, mg_iphi, factor=fac, state_from=psb)
    !This is mg_copy_to_tree from psb state
    !!!  replaces::  call mg_copy_to_tree(su_, mg_irhs, factor=-lambda)
    ! iw_from=iw_r_e
    ! iw_to=mg_iphi
    ! fac=fac
    ! do iigrid=1,igridstail; igrid=igrids(iigrid);
    !    pnode => igrid_to_node(igrid, mype)%node
    !    id    =  pnode%id
    !    lvl   =  mg%boxes(id)%lvl
    !    nc    =  mg%box_size_lvl(lvl)
    !    ! Include one layer of ghost cells on grid leaves
    !    {^IFONED
    !    mg%boxes(id)%cc(0:nc+1, iw_to) = fac * &
    !         psb(igrid)%w(ixMlo1-1:ixMhi1+1, iw_from)
    !    }
    !    {^IFTWOD
    !    mg%boxes(id)%cc(0:nc+1, 0:nc+1, iw_to) = fac * &
    !         psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, iw_from)
    !    }
    !    {^IFTHREED
    !    mg%boxes(id)%cc(0:nc+1, 0:nc+1, 0:nc+1, iw_to) = fac * &
    !         psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, &
    !         ixMlo3-1:ixMhi3+1, iw_from)
    !    }
    ! end do
 
    ! print*, '4'
 
 
    !>replace call set_rhs(mg, -1/dt, 0.0_dp)
    call mg_copy_to_tree(iw_r_e, mg_irhs, factor=-1/(dtfactor*dt_diff), state_from=psb)
    !This is mg_copy_to_tree from psb state
    !!!  replaces::  call mg_copy_to_tree(su_, mg_irhs, factor=-lambda)
    ! iw_from=iw_r_e
    ! iw_to=mg_irhs
    ! fac=-1/(dtfactor*dt_diff)
    ! do iigrid=1,igridstail; igrid=igrids(iigrid);
    !    pnode => igrid_to_node(igrid, mype)%node
    !    id    =  pnode%id
    !    lvl   =  mg%boxes(id)%lvl
    !    nc    =  mg%box_size_lvl(lvl)
    !    ! Include one layer of ghost cells on grid leaves
    !    {^IFONED
    !    mg%boxes(id)%cc(0:nc+1, iw_to) = fac * &
    !         psb(igrid)%w(ixMlo1-1:ixMhi1+1, iw_from)
    !    }
    !    {^IFTWOD
    !    mg%boxes(id)%cc(0:nc+1, 0:nc+1, iw_to) = fac * &
    !         psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, iw_from)
    !    }
    !    {^IFTHREED
    !    mg%boxes(id)%cc(0:nc+1, 0:nc+1, 0:nc+1, iw_to) = fac * &
    !         psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, &
    !         ixMlo3-1:ixMhi3+1, iw_from)
    !    }
    ! end do
 
    ! call phys_set_mg_bounds()
 
    ! print*, '5'
 
 
    if (time_advance)then
      call mg_restrict(mg, mg_iveps)
      call mg_fill_ghost_cells(mg, mg_iveps)
    endif
 
    ! print*, '6'
 
 
    call mg_fas_fmg(mg, .true., max_res=res)
    do n = 1, max_its
      !print*, n, res
      if (res < max_residual) exit
      call mg_fas_vcycle(mg, max_res=res)
    end do
 
    ! print*, '7'
 
 
    if (res .le. 0.d0) then
      if (diffcrash_resume) then
        if (mg%my_rank == 0) &
        write(*,*) it, ' resiudal zero ', res
        return
      endif
      if (mg%my_rank == 0) then
        print*, res
        error stop "Diffusion residual to zero"
      endif
    endif
 
    ! print*, '8'
 
 
    if (n == max_its + 1) then
      if (diffcrash_resume) then
        if (mg%my_rank == 0) &
        write(*,*) it, ' resiudal high ', res
        return
      endif
       if (mg%my_rank == 0) then
          print *, "Did you specify boundary conditions correctly?"
          print *, "Or is the variation in diffusion too large?"
          print*, n, res
          print *, mg%bc(1, mg_iphi)%bc_value, mg%bc(2, mg_iphi)%bc_value
       end if
       error stop "diffusion_solve: no convergence"
    end if
 
    ! print*, '9'
 
 
    !> Reset dt_diff when diffusion worked out
!0887 dt_diff = 0.d0
    dt_diff = 0.d0
 
    ! !This is mg_copy_from_tree_gc for psa state
    call mg_copy_from_tree_gc(mg_iphi, iw_r_e, state_to=psa)
    !This is mg_copy_from_tree_gc for psa state
    !!! replaces:: call mg_copy_from_tree_gc(mg_iphi, su_)
    ! iw_from=mg_iphi
    ! iw_to=iw_r_e
    ! do iigrid=1,igridstail; igrid=igrids(iigrid);
    !    pnode => igrid_to_node(igrid, mype)%node
    !    id    =  pnode%id
    !    lvl   =  mg%boxes(id)%lvl
    !    nc    =  mg%box_size_lvl(lvl)
    !    ! Include one layer of ghost cells on grid leaves
    !    {^IFONED
    !    psa(igrid)%w(ixMlo1-1:ixMhi1+1, iw_to) = &
    !         mg%boxes(id)%cc(0:nc+1, iw_from)
    !    }
    !    {^IFTWOD
    !    psa(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, iw_to) = &
    !         mg%boxes(id)%cc(0:nc+1, 0:nc+1, iw_from)
    !    }
    !    {^IFTHREED
    !    psa(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, &
    !         ixMlo3-1:ixMhi3+1, iw_to) = &
    !         mg%boxes(id)%cc(0:nc+1, 0:nc+1, 0:nc+1, iw_from)
    !    }
    ! end do
 
    ! print*, '10'
 
 
  end subroutine diffuse_e_rad_mg
 
 
  !> inplace update of psa==>F_im(psa)
  subroutine evaluate_e_rad_mg(qtC,psa)
    use mod_global_parameters
    use mod_ghostcells_update
    use mod_physics, only: phys_req_diagonal
 
    type(state), target :: psa(max_blocks)
    double precision, intent(in) :: qtC
 
    integer :: iigrid, igrid, level
    integer :: ixO^L
 
    call update_diffcoeff(psa)
 
    ixo^l=ixg^ll^lsub1;
    !$OMP PARALLEL DO PRIVATE(igrid)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
       ! call fld_get_diffcoef_central(psa(igrid)%w, psa(igrid)%w, psa(igrid)%x, ixG^LL, ixO^L)
       call put_diffterm_onegrid(ixg^ll,ixo^l,psa(igrid)%w)
    end do
    !$OMP END PARALLEL DO
 
    ! enforce boundary conditions for psa
    call getbc(qtc,0.d0,psa,1,nwflux+nwaux,phys_req_diagonal)
 
  end subroutine evaluate_e_rad_mg
 
  !> inplace update of psa==>F_im(psa)
  subroutine put_diffterm_onegrid(ixI^L,ixO^L,w)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(inout) :: w(ixI^S, 1:nw)
 
    double precision :: gradE(ixO^S),divF(ixO^S)
    double precision :: divF_h(ixO^S),divF_j(ixO^S)
    double precision :: diff_term(ixO^S)
 
    integer                       :: idir, jxO^L, hxO^L
 
    ! call mpistop("phys_evaluate_implicit not implemented for FLD")
 
    divf(ixo^s) = 0.d0
    do idir = 1,ndim
      hxo^l=ixo^l-kr(idir,^d);
      jxo^l=ixo^l+kr(idir,^d);
 
      divf_h(ixo^s) = w(ixo^s,i_diff_mg)*w(hxo^s,i_diff_mg)/(w(ixo^s,i_diff_mg) + w(hxo^s,i_diff_mg))*(w(hxo^s,iw_r_e) - w(ixo^s,iw_r_e))
      divf_j(ixo^s) = w(ixo^s,i_diff_mg)*w(jxo^s,i_diff_mg)/(w(ixo^s,i_diff_mg) + w(jxo^s,i_diff_mg))*(w(jxo^s,iw_r_e) - w(ixo^s,iw_r_e))
      divf(ixo^s) = divf(ixo^s) + 2.d0*(divf_h(ixo^s) + divf_j(ixo^s))/dxlevel(idir)**2
    enddo
 
    w(ixo^s,iw_r_e) = divf(ixo^s)
 
  end subroutine put_diffterm_onegrid
 
 
  !> Calculates cell-centered diffusion coefficient to be used in multigrid
  subroutine fld_get_diffcoef_central(w, wCT, x, ixI^L, ixO^L)
    use mod_global_parameters
    use mod_geometry
    use mod_usr_methods
 
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(inout) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: wCT(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
 
    double precision :: kappa(ixO^S), lambda(ixO^S), fld_R(ixO^S)
 
    double precision :: max_D(ixI^S), grad_r_e(ixI^S), rad_e(ixI^S)
    integer :: idir,i,j, ix^D
 
 
    call fld_get_opacity(wct, x, ixi^l, ixo^l, kappa)
    call fld_get_fluxlimiter(wct, x, ixi^l, ixo^l, lambda, fld_r)
 
    !> calculate diffusion coefficient
    w(ixo^s,i_diff_mg) = (const_c/unit_velocity)*lambda(ixo^s)/(kappa(ixo^s)*wct(ixo^s,iw_rho))
 
    where (w(ixo^s,i_diff_mg) .lt. 0.d0) &
      w(ixo^s,i_diff_mg) = smalldouble
 
    if (diff_coef_filter) then
      !call mpistop('Hold your bloody horses, not implemented yet ')
      call fld_smooth_diffcoef(w, ixi^l, ixo^l)
    endif
 
    if (associated(usr_special_diffcoef)) &
      call usr_special_diffcoef(w, wct, x, ixi^l, ixo^l)
 
  end subroutine fld_get_diffcoef_central
 
  subroutine update_diffcoeff(psa)
    use mod_global_parameters
 
    type(state), target :: psa(max_blocks)
 
    ! double precision :: wCT(ixG^LL,1:nw)
    integer :: iigrid, igrid, level
    integer :: ixO^L
 
    ixo^l=ixg^ll^lsub1;
    !$OMP PARALLEL DO PRIVATE(igrid)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
 
       ! wCT = psa(igrid)%w
        call fld_get_diffcoef_central(psa(igrid)%w, psa(igrid)%w, psa(igrid)%x, ixg^ll, ixo^l)
    end do
    !$OMP END PARALLEL DO
 
  end subroutine update_diffcoeff
 
  !> Use running average on Diffusion coefficient
  subroutine fld_smooth_diffcoef(w, ixI^L, ixO^L)
    use mod_global_parameters
 
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(inout) :: w(ixI^S, 1:nw)
 
    double precision :: tmp_D(ixI^S), filtered_D(ixI^S)
    integer :: ix^D, filter, idim
 
    if (size_d_filter .lt. 1) call mpistop("D filter of size < 1 makes no sense")
    if (size_d_filter .gt. nghostcells) call mpistop("D filter of size > nghostcells makes no sense")
 
    tmp_d(ixo^s) = w(ixo^s,i_diff_mg)
    filtered_d(ixo^s) = zero
 
    do filter = 1,size_d_filter
      {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_d_filter\}
        do idim = 1,ndim
          filtered_d(ix^d) = filtered_d(ix^d) &
                           + tmp_d(ix^d+filter*kr(idim,^d)) &
                           + tmp_d(ix^d-filter*kr(idim,^d))
        enddo
      {enddo\}
    enddo
 
    {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_d_filter\}
      tmp_d(ix^d) = (tmp_d(ix^d)+filtered_d(ix^d))/(1+2*size_d_filter*ndim)
    {enddo\}
 
    w(ixo^s,i_diff_mg) = tmp_d(ixo^s)
  end subroutine fld_smooth_diffcoef
 
 
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!!!!!!!!!!!!!! Gas-Rad Energy interaction
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  !> This subroutine calculates the radiation heating, radiation cooling
  !> and photon tiring using an implicit scheme.
  !> These sourceterms are applied using the root-finding of a 4th order polynomial
  !> This routine loops over every cell in the domain
  !> and computes the coefficients of the polynomials in every cell
  subroutine energy_interaction(w, wCT, x, ixI^L, ixO^L)
    use mod_global_parameters
    use mod_geometry
    use mod_physics
    use mod_usr_methods
 
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    double precision, intent(in)    :: wCT(ixI^S,1:nw)
    double precision, intent(inout) :: w(ixI^S,1:nw)
 
    double precision :: a1(ixO^S), a2(ixO^S)
    double precision :: c0(ixO^S), c1(ixO^S)
    double precision :: e_gas(ixO^S), E_rad(ixO^S)
    double precision :: kappa(ixO^S)
    double precision :: sigma_b
 
    integer :: i,j,idir,ix^D
 
    !> e_gas is the INTERNAL ENERGY without KINETIC ENERGY
    ! if (.not. block%e_is_internal) then
      e_gas(ixo^s) = wct(ixo^s,iw_e) - half*sum(wct(ixo^s, iw_mom(:))**2, dim=ndim+1)/wct(ixo^s, iw_rho)
    ! else
    !   e_gas(ixO^S) = wCT(ixO^S,iw_e)
    ! endif
 
    {do ix^d=ixomin^d,ixomax^d\ }
      e_gas(ix^d) = max(e_gas(ix^d),small_pressure/(fld_gamma-1.d0))
    {enddo\}
 
    e_rad(ixo^s) = wct(ixo^s,iw_r_e)
 
    if (associated(usr_special_opacity_qdot)) then
      call usr_special_opacity_qdot(ixi^l,ixo^l,w,x,kappa)
    else
      call fld_get_opacity(wct, x, ixi^l, ixo^l, kappa)
    endif
 
    sigma_b = const_rad_a*const_c/4.d0*(unit_temperature**4.d0)/(unit_velocity*unit_pressure)
 
    if (fld_interaction_method .eq. 'Instant') then
      a1(ixo^s) = const_rad_a*(fld_mu*const_mp/const_kb*(fld_gamma-1))**4/(wct(ixo^s,iw_rho)*unit_density)**4 &
                  /unit_pressure**3
      a2(ixo^s) = e_gas(ixo^s) + e_rad(ixo^s)
 
      c0(ixo^s) = a2(ixo^s)/a1(ixo^s)
      c1(ixo^s) = 1.d0/a1(ixo^s)
    else
      !> Calculate coefficients for polynomial
      a1(ixo^s) = 4*kappa(ixo^s)*sigma_b*(fld_gamma-one)**4/wct(ixo^s,iw_rho)**3*dt
      a2(ixo^s) = (const_c/unit_velocity)*kappa(ixo^s)*wct(ixo^s,iw_rho)*dt
 
      c0(ixo^s) = ((one + a2(ixo^s))*e_gas(ixo^s) + a2(ixo^s)*e_rad(ixo^s))/a1(ixo^s)
      c1(ixo^s) = (one + a2(ixo^s))/a1(ixo^s)
    endif
 
    !> Loop over every cell for rootfinding method
    {do ix^d=ixomin^d,ixomax^d\ }
      select case(fld_interaction_method)
      case('Bisect')
        call bisection_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
      case('Newton')
        call newton_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
      case('Halley')
        call halley_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
      case('Instant')
        call halley_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
      case default
        call mpistop('root-method not known')
      end select
    {enddo\}
 
    if (fld_interaction_method .eq. 'Instant') then
      e_rad(ixo^s) = a2(ixo^s) - e_gas(ixo^s)
    else
      !> advance E_rad
      e_rad(ixo^s) = (a1(ixo^s)*e_gas(ixo^s)**4.d0 + e_rad(ixo^s))/(one + a2(ixo^s))
    endif
 
    !> new w = w + dt f(wCT)
    !> e_gas,E_rad = wCT + dt f(wCT)
    !> dt f(wCT) = e_gas,E_rad - wCT
    !> new w = w +  e_gas,E_rad - wCT
 
    !> WAIT A SECOND?! DOES THIS EVEN WORK WITH THESE KIND OF IMPLICIT METHODS?
    !> NOT QUITE SURE TO BE HONEST
    !> IS IT POSSIBLE TO SHUT DOWN SOURCE SPLITTING FOR JUST THIS TERM?
    !> FIX BY PERFORMING Energy_interaction on (w,w,...)
 
    ! call mpistop('This still has to be fixed somehow')
 
    !> Update gas-energy in w, internal + kinetic
    ! w(ixO^S,iw_e) = w(ixO^S,iw_e) + e_gas(ixO^S) - wCT(ixO^S,iw_e)
    ! {do ix^D=ixOmin^D,ixOmax^D\ }
    !   e_gas(ix^D) = max(e_gas(ix^D),small_pressure/(fld_gamma-1.d0))
    ! {enddo\}
 
    !{do ix^D=ixOmin^D,ixOmax^D\ }
    !  w(ix^D,iw_e) = max(e_gas(ix^D),small_pressure/(fld_gamma - 1))
    !{enddo\}
    w(ixo^s,iw_e) = e_gas(ixo^s)
 
    !> Beginning of module substracted wCT Ekin
    !> So now add wCT Ekin
    ! if (.not. block%e_is_internal) then
    w(ixo^s,iw_e) = w(ixo^s,iw_e) + half*sum(wct(ixo^s, iw_mom(:))**2, dim=ndim+1)/wct(ixo^s, iw_rho)
    ! else
    !   w(ixO^S,iw_e) = w(ixO^S,iw_e)
    ! endif
 
    !> Update rad-energy in w
    ! w(ixO^S,iw_r_e) = w(ixO^S,iw_r_e) + E_rad(ixO^S) - wCT(ixO^S,iw_r_e)
    w(ixo^s,iw_r_e) = e_rad(ixo^s)
 
  end subroutine energy_interaction
 
 
  !> Find the root of the 4th degree polynomial using the bisection method
  subroutine bisection_method(e_gas, E_rad, c0, c1)
    use mod_global_parameters
 
    double precision, intent(in)    :: c0, c1
    double precision, intent(in)    :: E_rad
    double precision, intent(inout) :: e_gas
 
    double precision :: bisect_a, bisect_b, bisect_c
    integer :: n, max_its
 
    n = 0
    max_its = 1d2
 
    bisect_a = zero
    bisect_b = min(abs(c0/c1),abs(c0)**(1.d0/4.d0))+smalldouble
 
    ! do while (abs(Polynomial_Bisection(bisect_b, c0, c1)-Polynomial_Bisection(bisect_a, c0, c1))&
    !    .ge. fld_bisect_tol*min(e_gas,E_rad))
    do while (abs(bisect_b-bisect_a) .ge. fld_bisect_tol*min(e_gas,e_rad))
      bisect_c = (bisect_a + bisect_b)/two
 
      n = n +1
      if (n .gt. max_its) then
        goto 2435
        call mpistop('No convergece in bisection scheme')
      endif
 
      if (polynomial_bisection(bisect_a, c0, c1)*&
      polynomial_bisection(bisect_b, c0, c1) .lt. zero) then
 
        if (polynomial_bisection(bisect_a, c0, c1)*&
        polynomial_bisection(bisect_c, c0, c1) .lt. zero) then
          bisect_b = bisect_c
        elseif (polynomial_bisection(bisect_b, c0, c1)*&
        polynomial_bisection(bisect_c, c0, c1) .lt. zero) then
          bisect_a = bisect_c
        elseif (polynomial_bisection(bisect_a, c0, c1) .eq. zero) then
          bisect_b = bisect_a
          bisect_c = bisect_a
          goto 2435
        elseif (polynomial_bisection(bisect_b, c0, c1) .eq. zero) then
          bisect_a = bisect_b
          bisect_c = bisect_b
          goto 2435
        elseif (polynomial_bisection(bisect_c, c0, c1) .eq. zero) then
          bisect_a = bisect_c
          bisect_b = bisect_c
          goto 2435
        else
          call mpistop("Problem with fld bisection method")
        endif
      elseif (polynomial_bisection(bisect_a, c0, c1) &
        - polynomial_bisection(bisect_b, c0, c1) .lt. fld_bisect_tol*polynomial_bisection(bisect_a, c0, c1)) then
        goto 2435
      else
        bisect_a = e_gas
        bisect_b = e_gas
        print*, "IGNORING GAS-RAD ENERGY EXCHANGE ", c0, c1
 
        print*, polynomial_bisection(bisect_a, c0, c1), polynomial_bisection(bisect_b, c0, c1)
 
        if (polynomial_bisection(bisect_a, c0, c1) .le. smalldouble) then
          bisect_b = bisect_a
        elseif (polynomial_bisection(bisect_a, c0, c1) .le. smalldouble) then
          bisect_a = bisect_b
        endif
 
        goto 2435
 
      endif
    enddo
 
      2435 e_gas = (bisect_a + bisect_b)/two
  end subroutine bisection_method
 
  !> Find the root of the 4th degree polynomial using the Newton-Ralphson method
  subroutine newton_method(e_gas, E_rad, c0, c1)
    use mod_global_parameters
 
    double precision, intent(in)    :: c0, c1
    double precision, intent(in)    :: E_rad
    double precision, intent(inout) :: e_gas
 
    double precision :: xval, yval, der, deltax
 
    integer :: ii
 
    yval = bigdouble
    xval = e_gas
    der = one
    deltax = one
 
    ii = 0
    !> Compare error with dx = dx/dy dy
    do while (abs(deltax) .gt. fld_bisect_tol)
      yval = polynomial_bisection(xval, c0, c1)
      der = dpolynomial_bisection(xval, c0, c1)
      deltax = -yval/der
      xval = xval + deltax
      ii = ii + 1
      if (ii .gt. 1d3) then
        print*, 'skip to bisection algorithm'
        call bisection_method(e_gas, e_rad, c0, c1)
        return
      endif
    enddo
 
    e_gas = xval
  end subroutine newton_method
 
  !> Find the root of the 4th degree polynomial using the Halley method
  subroutine halley_method(e_gas, E_rad, c0, c1)
    use mod_global_parameters
 
    double precision, intent(in)    :: c0, c1
    double precision, intent(in)    :: E_rad
    double precision, intent(inout) :: e_gas
 
    double precision :: xval, yval, der, dder, deltax
 
    integer :: ii
 
    yval = bigdouble
    xval = e_gas
    der = one
    dder = one
    deltax = one
 
    ii = 0
    !> Compare error with dx = dx/dy dy
    do while (abs(deltax) .gt. fld_bisect_tol)
      yval = polynomial_bisection(xval, c0, c1)
      der = dpolynomial_bisection(xval, c0, c1)
      dder = ddpolynomial_bisection(xval, c0, c1)
      deltax = -two*yval*der/(two*der**2 - yval*dder)
      xval = xval + deltax
      ii = ii + 1
      if (ii .gt. 1d3) then
        ! call mpistop('Halley did not convergggge')
        call newton_method(e_gas, e_rad, c0, c1)
        return
      endif
    enddo
 
    e_gas = xval
  end subroutine halley_method
 
  !> Evaluate polynomial at argument e_gas
  function polynomial_bisection(e_gas, c0, c1) result(val)
    use mod_global_parameters
 
    double precision, intent(in) :: e_gas
    double precision, intent(in) :: c0, c1
    double precision :: val
 
    val = e_gas**4.d0 + c1*e_gas - c0
  end function polynomial_bisection
 
  !> Evaluate first derivative of polynomial at argument e_gas
  function dpolynomial_bisection(e_gas, c0, c1) result(der)
    use mod_global_parameters
 
    double precision, intent(in) :: e_gas
    double precision, intent(in) :: c0, c1
    double precision :: der
 
    der = 4.d0*e_gas**3.d0 + c1
  end function dpolynomial_bisection
 
  !> Evaluate second derivative of polynomial at argument e_gas
  function ddpolynomial_bisection(e_gas, c0, c1) result(dder)
    use mod_global_parameters
 
    double precision, intent(in) :: e_gas
    double precision, intent(in) :: c0, c1
    double precision :: dder
 
    dder = 4.d0*3.d0*e_gas**2.d0
  end function ddpolynomial_bisection
 
  !> Calculate gradient of a scalar q within ixL in direction idir
  !> difference with gradient is gradq(ixO^S), NOT gradq(ixI^S)
  subroutine gradiento(q,x,ixI^L,ixO^L,idir,gradq,n)
    use mod_global_parameters
 
    integer, intent(in)             :: ixI^L, ixO^L, idir
    integer, intent(in)             :: n
 
    double precision, intent(in)    :: q(ixI^S), x(ixI^S,1:ndim)
    double precision, intent(out)   :: gradq(ixO^S)
    integer                         :: jxO^L, hxO^L
 
    ! hxO^L=ixO^L-n*kr(idir,^D);
    ! jxO^L=ixO^L+n*kr(idir,^D);
    !
    ! if (n .gt. nghostcells) call mpistop("gradientO stencil too wide")
    !
    ! gradq(ixO^S)=(q(jxO^S)-q(hxO^S))/(2*n*dxlevel(idir))
    ! gradq(ixO^S)=(q(jxO^S)-q(hxO^S))/(x(jxO^S,idir)-x(hxO^S,idir))
 
    !> Using higher order derivatives with wider stencil according to:
    !> https://en.wikipedia.org/wiki/Finite_difference_coefficient
 
     if (n .gt. nghostcells) then
       call mpistop("gradientO stencil too wide")
     elseif (n .eq. 1) then
       hxo^l=ixo^l-kr(idir,^d);
       jxo^l=ixo^l+kr(idir,^d);
       gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,idir)-x(hxo^s,idir))
     elseif (n .eq. 2) then
       gradq(ixo^s) = 0.d0
       !> coef 2/3
       hxo^l=ixo^l-kr(idir,^d);
       jxo^l=ixo^l+kr(idir,^d);
       gradq(ixo^s) = gradq(ixo^s) + 2.d0/3.d0*(q(jxo^s)-q(hxo^s))
       !> coef -1/12
       hxo^l=ixo^l-2*kr(idir,^d);
       jxo^l=ixo^l+2*kr(idir,^d);
       gradq(ixo^s) = gradq(ixo^s) - 1.d0/12.d0*(q(jxo^s)-q(hxo^s))
       !> divide by dx
       gradq(ixo^s) = gradq(ixo^s)/dxlevel(idir)
     elseif (n .eq. 3) then
       gradq(ixo^s) = 0.d0
       !> coef 3/4
       hxo^l=ixo^l-kr(idir,^d);
       jxo^l=ixo^l+kr(idir,^d);
       gradq(ixo^s) = gradq(ixo^s) + 3.d0/4.d0*(q(jxo^s)-q(hxo^s))
       !> coef -3/20
       hxo^l=ixo^l-2*kr(idir,^d);
       jxo^l=ixo^l+2*kr(idir,^d);
       gradq(ixo^s) = gradq(ixo^s) - 3.d0/20.d0*(q(jxo^s)-q(hxo^s))
       !> coef 1/60
       hxo^l=ixo^l-3*kr(idir,^d);
       jxo^l=ixo^l+3*kr(idir,^d);
       gradq(ixo^s) = gradq(ixo^s) + 1.d0/60.d0*(q(jxo^s)-q(hxo^s))
       !> divide by dx
       gradq(ixo^s) = gradq(ixo^s)/dxlevel(idir)
     else
       call mpistop("gradientO stencil unknown")
     endif
 
  end subroutine gradiento
 
end module mod_fld