mod_fld.t

Go to the documentation of this file.

!> Module for flux limited diffusion (FLD)-approximation in Radiation-(Magneto)hydrodynamics simulations
!>
!> Full description of RHD-FLD in
!> Moens N., Sundqvist J.O., El Mellah I., Poniatowski L., Teunissen J. & Keppens R. 2022, A&A 657, A81
!> Radiation-hydrodynamics with MPI-AMRVAC . Flux-limited diffusion
!> doi:10.1051/0004-6361/202141023
!>
!> Full description for RMHD-FLD in
!> N. Narechania, R. Keppens, A. ud-Doula, N. Moens & J. Sundqvist 2025, A&A 696, A131
!> doi:10.1051/0004-6361/202452208
!> Radiation-magnetohydrodynamics with MPI-AMRVAC using flux-limited diffusion
 
module mod_fld
    use mod_comm_lib, only: mpistop
    use mod_geometry
    implicit none
    !> source split for energy interact and radforce:
    logical :: fld_radforce_split = .false.
    !> switch to handle photon tiring explicit or implicit
    logical :: fld_tiring_explicit = .false.
    !> Opacity value when using constant opacity
    double precision, public :: fld_kappa0 = 0.0d0
    !> 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 radiative Energy diffusion
    double precision, public :: fld_diff_tol = 1.d-4
    !> handling ramp-up phase with explicit diffusion dt limit, slowly boosted
    logical :: fld_slowsteps = .false.
    double precision, public :: fld_boost_dt = 0.0d0
    !> switches for using changed cmax-cmin bounds
    logical :: fld_bound_diff = .true.
    !> switches for local debug purposes
    logical :: fld_force_mg_converged = .true.
    logical :: fld_debug,fld_no_mg
    double precision, public :: fld_cnorm = 0.0d0
    !> switches for opacity
    character(len=40)  :: fld_opacity_law = 'const'
    character(len=40) :: fld_opal_table = 'Y09800' 
    !> flux limiter choice
    character(len=40) :: fld_fluxlimiter = 'Pomraning'
    !> diffusion coefficient for multigrid method
    integer :: i_diff_mg
    !> diffusion coefficient stencil control
    integer :: nth_for_diff_mg
    !> Which method to find the root for the energy interaction polynomial
    character(len=40) :: fld_interaction_method = 'Halley'
    !> Abstract interface for the gas-EoS getters the radiation fluid needs
    !> (same shape as thermal_conduction's get_var_subr).
    abstract interface
      subroutine fld_get_var(w,x,ixI^L,ixO^L,res)
        use mod_global_parameters
        integer, intent(in)          :: ixI^L, ixO^L
        double precision, intent(in) :: w(ixI^S,nw)
        double precision, intent(in) :: x(ixI^S,1:ndim)
        double precision, intent(out):: res(ixI^S)
      end subroutine fld_get_var
    end interface
 
    !> Radiation fluid object: gas-EoS callbacks the FLD module needs, wired by
    !> the physics module at link time (mirrors tc_fluid / rc_fluid). The
    !> instance lives in mod_(m)hd_phys and is threaded in here as `fl`.
    type fld_fluid
      !> adiabatic index (constant for FI; a snapshot of eos%gamma)
      double precision :: gamma
      !> gas temperature from primitive pressure, T = p/(R*rho)
      procedure(fld_get_var), pointer, nopass :: get_tgas => null()
      !> R factor (mean-molecular-weight gas constant) for the emission term
      procedure(fld_get_var), pointer, nopass :: get_rfactor => null()
    end type fld_fluid
 
    !> public methods
    !> these are called in mod_hd_phys or mod_mhd_phys
    public :: fld_fluid
    public :: fld_init
    public :: fld_get_radpress
    public :: fld_get_diffcoef_central
    public :: add_fld_rad_force
    public :: fld_radforce_get_dt
    public :: fld_get_local_invtauc
    !> wired to physics-module wrappers (which inject fld_fl) in mod_(m)hd_phys
    public :: fld_implicit_update
    public :: fld_evaluate_implicit
    public :: fld_set_mg_bounds
    ! these are made public for mod_usr purposes and diagnostics
    public :: fld_get_radflux
    public :: fld_get_fluxlimiter
    public :: fld_get_fluxlimiter_prim
    public :: fld_get_opacity_prim
  contains
 
  !> 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_radforce_split, &
    fld_bisect_tol, fld_diff_tol, fld_opacity_law, fld_fluxlimiter, &
    fld_interaction_method, fld_opal_table, nth_for_diff_mg, &
    fld_debug,fld_no_mg,fld_force_mg_converged, fld_slowsteps,fld_boost_dt, &
    fld_tiring_explicit,fld_cnorm,fld_bound_diff
 
    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 and adimensionalise kappa
  subroutine fld_init(gas_eos_type, gas_ion_energy, gas_pi_table)
    use mod_global_parameters
    use mod_variables
    use mod_physics
    use mod_opal_opacity, only: init_opal_table
    use mod_multigrid_coupling
    character(len=*), intent(in) :: gas_eos_type, gas_pi_table
    logical, intent(in) :: gas_ion_energy
 
    nth_for_diff_mg=1
    fld_debug=.false.
    fld_no_mg=.false.
    ! initialize constant opacity with free electron Thomson scattering value
    fld_kappa0=const_kappae
    call fld_params_read(par_files)
    ! Reject unsupported gas closures before loading opacity tables or using
    ! EOS callbacks. These are implementation limits, not physical exclusions.
    if(trim(gas_eos_type)=='LTE') call mpistop("FLD is not yet available with eos_type='LTE': "//&
         'the radiation exchange takes the gas temperature as (gamma-1) e/(rho R)')
    if(gas_ion_energy.and.trim(gas_eos_type)=='PI') then
      call mpistop('FLD currently does not support PI with ionE=T: '//&
           'the radiation exchange solver assumes ideal-gas internal energy')
    end if
    if(trim(fld_opacity_law)=='opal') then
      if(trim(gas_eos_type)=='PI'.and.trim(gas_pi_table)=='prominence') then
        call mpistop('FLD with PI prominence and OPAL opacity is not implemented: '//&
             'the opacity temperature callback is incompatible with primitive states')
      end if
    end if
    ! sanity checks on input
    if(fld_kappa0<smalldouble)then
       if(mype==0) print *,'fld_kappa0=',fld_kappa0
       call mpistop("please set the constant opacity to a reasonable value")
    endif
    if(fld_bisect_tol<smalldouble)then
       if(mype==0) print *,'fld_bisect_tol=',fld_bisect_tol
       call mpistop("convergence tolerance for root solver too strict")
    endif
    if(fld_tiring_explicit)then
       if(mype==0) print *,'Will do photon tiring explicit!!!'
    endif
    if(.not.fld_no_mg)then
       if(fld_diff_tol<smalldouble)then
         if(mype==0) print *,'fld_diff_tol=',fld_diff_tol
         call mpistop("convergence tolerance for MG solver too strict")
      endif
      select case(nth_for_diff_mg)
        case(1)
         ! no need for stencil extension
        case(2)
         ! need for stencil extension
         phys_wider_stencil=1 
       case default
         call mpistop("nth_for_diff_mg must be 1 or 2")
       end select
       ! phys_implicit_update / phys_evaluate_implicit are wired in the physics
       ! module (mod_(m)hd_phys) to wrappers that inject the fld_fl object;
       ! set_mg_bounds needs no fluid object so it is bound here directly.
       phys_set_mg_bounds     => fld_set_mg_bounds
       ! store the diffusion coefficient as extra variable (needed in mg vhelmholtz)
       i_diff_mg = var_set_extravar("D", "D")
       use_multigrid = .true.
       ! use multigrid to solve a helmholtz equation with variable coefficient
       ! this is stored also as extra variable in the mg solver as mg_iveps
       mg%n_extra_vars = 1
       mg%operator_type = mg_vhelmholtz
       ! choice of smoother: can be mg_smoother_gs or gsrb (latter recommended)
       mg%smoother_type = mg_smoother_gsrb
    endif
    !> Read in opacity table if necesary
    if(trim(fld_opacity_law) .eq. 'opal') then
      if(si_unit)call mpistop("adjust opal module with SI-cgs conversions for SI - or use cgs!")
      call init_opal_table(fld_opal_table)
    endif
  end subroutine fld_init
 
  !> Set the boundaries for the diffusion of E
  subroutine fld_set_mg_bounds
    use mod_global_parameters
    use mod_multigrid_coupling
    use mod_usr_methods
 
    integer :: ib
 
    ! Set boundary conditions for the multigrid solver
    do ib = 1, 2*ndim
       select case (typeboundary(iw_r_e, ib))
       case (bc_symm)
          ! d/dx u = 0
          mg%bc(ib, mg_iphi)%bc_type = mg_bc_neumann
          mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
       case (bc_asymm)
          ! u = 0
          mg%bc(ib, mg_iphi)%bc_type = mg_bc_dirichlet
          mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
       case (bc_cont)
          ! d/dx u = 0
          mg%bc(ib, mg_iphi)%bc_type = mg_bc_neumann
          mg%bc(ib, mg_iphi)%bc_value = 0.0_dp
       case (bc_periodic)
          ! Nothing to do here, this is picked up through periodB variable
       case (bc_special)
          if(mype==0)then
             print *,'Special boundary for Erad needs specific user-set MG BC treatment'
             print *,'                 and this could be through usr_special_mg_bc call'
          endif
          if (associated(usr_special_mg_bc)) then 
             call usr_special_mg_bc(ib)
          endif
       case default
          call mpistop("divE_multigrid warning: unknown b.c. ")
       end select
       ! Neumann on diffusion coefficient is needed on all lower grid levels
       ! d/dx u = 0
       mg%bc(ib, mg_iveps)%bc_type = mg_bc_neumann
       mg%bc(ib, mg_iveps)%bc_value = 0.0_dp
    end do
 
  end subroutine fld_set_mg_bounds
 
 
  !> 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 and its work added explicitly
  !> and the energy interaction term combined with photon tiring using an implicit update
  subroutine add_fld_rad_force(qdt,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active,fl)
    use mod_constants
    use mod_global_parameters
    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),wctprim(ixi^s,1:nw)
    double precision, intent(inout) :: w(ixi^s,1:nw)
    logical, intent(in) :: qsourcesplit
    logical, intent(inout) :: active
    type(fld_fluid), intent(in)     :: fl
 
    integer :: idir,jdir,nth_for_fld,ix^d
    double precision, dimension(ixI^S) :: a1,a2,a3,c0,c1,kappa
    double precision, dimension(ixI^S) :: e_gas,e_rad,tmp
    double precision, dimension(ixI^S,1:ndim,1:ndim) :: div_v,edd
 
    !> Calculate and add sourceterms
    if(qsourcesplit .eqv. fld_radforce_split) then
      active = .true.
      nth_for_fld=2
      ! store here lambda in a1 and fld_R in a2
      call fld_get_eddington(wctprim,x,ixi^l,ixo^l,edd,a1,a2,nth_for_fld,fl)
 
      !> Photon tiring : calculate tensor grad v (named div_v here)
      ! NOTE: This is ok for uniform Cartesian only!!!!!
      ! TODO: introduce gradient of vector in geometry module and call that one
      do idir = 1,ndim
        do jdir = 1,ndim
          call gradient(wctprim(ixi^s,iw_mom(jdir)),ixi^l,ixo^l,idir,tmp)
          div_v(ixo^s,idir,jdir) = tmp(ixo^s)
        enddo
      enddo
      ! perform contraction fe : grad(v) with fe eddington tensor
      {^ifoned
      a3(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1)
      }
      {^iftwod
      a3(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
      a3(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)
      }
 
      do idir = 1,ndim
        call gradient(wctprim(ixi^s,iw_r_e),ixi^l,ixo^l,idir,tmp,nth_for_fld)
        ! Radiation force = kappa*rho/c *Flux = lambda gradE
        ! recycle grad E to store -lambda (grad E)_i
        tmp(ixo^s) = -a1(ixo^s)*tmp(ixo^s)
        !> Momentum equation source term
        w(ixo^s,iw_mom(idir)) = w(ixo^s,iw_mom(idir))+ qdt*tmp(ixo^s)
        !> Energy equation source term 
        w(ixo^s,iw_e) = w(ixo^s,iw_e) + qdt*wctprim(ixo^s,iw_mom(idir))*tmp(ixo^s)
      enddo
      !> photon tiring when handled explicitly
      if(fld_tiring_explicit)then
          w(ixo^s,iw_r_e) = w(ixo^s,iw_r_e) - qdt*wctprim(ixo^s,iw_r_e)*a3(ixo^s)
          a3(ixo^s)=0.0d0
      endif
 
      call get_and_check_egas_erad_from_conserved(w,ixi^l,ixo^l,e_gas,e_rad)
      ! BEGIN radiative exchange part with or without photon tiring included
         call fld_get_opacity_prim(wctprim,x,ixi^l,ixo^l,kappa,fl)
         !> Coefficients for the polynomial in Moens et al. 2022, eq 37. but with photon tiring (a3)
         !> FI emission term: a1 carries arad*T^4 with T=(gamma-1)*e_gas/(rho*R) folded
         !> into the e_gas^4 root-find below. gamma and R come from the EoS via fl.
         call fl%get_Rfactor(wct,x,ixi^l,ixo^l,tmp)
         a1(ixo^s) = qdt*kappa(ixo^s)*c_norm*arad_norm*(fl%gamma-one)**4/(wct(ixo^s,iw_rho)**3*tmp(ixo^s)**4)
         a2(ixo^s) = c_norm*kappa(ixo^s)*wct(ixo^s,iw_rho)*qdt
         a3(ixo^s) = a3(ixo^s)*qdt
 
         c0(ixo^s) = ((one+a2(ixo^s)+a3(ixo^s))*e_gas(ixo^s)+a2(ixo^s)*e_rad(ixo^s))/a1(ixo^s)/(one+a3(ixo^s))
         c1(ixo^s) = (one+a2(ixo^s)+a3(ixo^s))/a1(ixo^s)/(one+a3(ixo^s))
 
         !> 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),c0(ix^d),c1(ix^d))
           case('Newton')
            call newton_method(e_gas(ix^d),c0(ix^d),c1(ix^d))
           case('Halley')
            call halley_method(e_gas(ix^d),c0(ix^d),c1(ix^d))
           case default
            call mpistop('root-method not known')
           end select 
         {enddo\}
 
         e_rad(ixo^s) = (a1(ixo^s)*e_gas(ixo^s)**4.d0+e_rad(ixo^s))/(one+a2(ixo^s)+a3(ixo^s))
 
         if(check_small_values.and..not.fix_small_values)then
          {do ix^db= ixomin^db,ixomax^db\}
           if(e_gas(ix^d)<small_e.or.e_rad(ix^d)<small_r_e) then
              write(*,*) "Error in FLD add_fld_rad_force: small value"
              write(*,*) "of internal or radiation energy density after exchange"
              write(*,*) "Iteration: ", it, " Time: ", global_time
              write(*,*) "Location: ", x(ix^d,:)
              write(*,*) "Cell number: ", ix^d
              write(*,*) "internal energy density  is=",e_gas(ix^d)," versus   small_e=",small_e
              write(*,*) "radiation energy density is=",e_rad(ix^d)," versus small_r_e=",small_r_e
              call mpistop("FLD error:May need to turn on fixes")
           end if
          {end do\}
         endif
 
        if(fix_small_values)then
          {do ix^d = ixomin^d,ixomax^d\ }
            e_gas(ix^d) = max(e_gas(ix^d),small_e)
            e_rad(ix^d) = max(e_rad(ix^d),small_r_e)
          {enddo\}
        endif
 
      ! END radiative exchange part with or without photon tiring included
 
      !> Update gas-energy in w, internal + kinetic
      w(ixo^s,iw_e) = e_gas(ixo^s)
      w(ixo^s,iw_e) = w(ixo^s,iw_e)+half*sum(w(ixo^s,iw_mom(:))**2,dim=ndim+1)/w(ixo^s,iw_rho)
      if(allocated(iw_mag)) then
        w(ixo^s,iw_e) = w(ixo^s,iw_e)+half*sum(w(ixo^s,iw_mag(:))**2,dim=ndim+1)
      endif
      !> Update rad-energy in w
      w(ixo^s,iw_r_e) = e_rad(ixo^s)
    end if
  end subroutine add_fld_rad_force
 
  !> maybe needed for restricting the cmax/cmin in diffusion limit
  !> CURRENTLY UNUSED
  !> NOTE: w is primitive on entry
  subroutine fld_get_local_invtauc(w,ixI^L,ixO^L,dx^D,x,invtauc,fl)
    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(out)   :: invtauc(ixi^s)
    type(fld_fluid), intent(in)     :: fl
 
    double precision :: dxinv(1:ndim)
    double precision :: kappa(ixi^s)
 
    call fld_get_opacity_prim(w, x, ixi^l, ixo^l, kappa, fl)
    ^d&dxinv(^d)=one/dx^d;
    invtauc(ixo^s)=4.0d0*minval(dxinv(1:ndim))/(3.0d0*kappa(ixo^s)*w(ixo^s,iw_rho))
  end subroutine fld_get_local_invtauc
 
  !> get dt limit for radiation force and FLD explicit source additions
  !> NOTE: w is primitive on entry
  subroutine fld_radforce_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x,fl)
    use mod_global_parameters
    use mod_usr_methods
    use mod_geometry
    use mod_physics, only: phys_get_csrad2
 
    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
    type(fld_fluid), intent(in)     :: fl
 
    integer          :: idim,idims,nth_for_fld
    double precision :: dxinv(1:ndim), max_grav
    double precision :: lambda(ixi^s),fld_r(ixi^s)
    double precision :: tmp(ixi^s),boostfactor
    double precision :: max_diff,max_diff_coef,max_diff_dim,dtdifflimit
    double precision :: cmax(ixi^s),cmaxtot(ixi^s),courantmaxtots
 
    if(fld_debug)print *,'DT limit on entry to radforce_get_dt=',dtnew
    nth_for_fld=2
    call fld_get_fluxlimiter_prim(w,x,ixi^l,ixo^l,lambda,fld_r,nth_for_fld,fl)
    if(slab_uniform) then
      ^d&dxinv(^d)=one/dx^d;
      do idim = 1, ndim
        call gradient(w(ixi^s,iw_r_e),ixi^l,ixo^l,idim,tmp,nth_for_fld)
        max_grav = maxval(dabs(-lambda(ixo^s)*tmp(ixo^s)/w(ixo^s,iw_rho)))
        max_grav = max(max_grav, epsilon(1.0d0))
        dtnew = min(dtnew, 1.0d0 / dsqrt(max_grav * dxinv(idim)))
      end do
    else
      do idim = 1, ndim
        call gradient(w(ixi^s,iw_r_e),ixi^l,ixo^l,idim,tmp,nth_for_fld)
        max_grav = maxval(dabs(-lambda(ixo^s)*tmp(ixo^s)/w(ixo^s,iw_rho))/block%ds(ixo^s,idim))
        max_grav = max(max_grav, epsilon(1.0d0))
        dtnew = min(dtnew, 1.0d0 / dsqrt(max_grav))
      end do
    endif
    if(fld_debug)print *,'DT limit after RADFORCE eff grav=',dtnew
 
    if(fld_slowsteps) then
      boostfactor=1.0d0
      call fld_get_opacity_prim(w, x, ixi^l, ixo^l, tmp, fl)
      max_diff=0.0d0
      if(slab_uniform) then
         max_diff_coef = maxval(c_norm*lambda(ixo^s)/(w(ixo^s,iw_rho)*tmp(ixo^s)))
         ^d&dxinv(^d)=one/dx^d;
         do idim = 1, ndim
            max_diff_dim = max(2.0d0*ndim*max_diff_coef*dxinv(idim)**2, epsilon(1.0d0))
            max_diff = max(max_diff,max_diff_dim)
         end do
      else
         do idim = 1, ndim
            max_diff_coef = maxval(c_norm*lambda(ixo^s)/(w(ixo^s,iw_rho)*tmp(ixo^s))/block%ds(ixo^s,idim)**2)
            max_diff_dim = max(2.0d0*ndim*max_diff_coef, epsilon(1.0d0))
            max_diff = max(max_diff,max_diff_dim)
         end do
      endif
      dtdifflimit=1.0d0/max_diff
      if(fld_debug) print *,'DT limit from diffusion is actually=',dtdifflimit
      if(slowsteps>it-it_init+1) then
         boostfactor=1.0d0+fld_boost_dt*(dble(it-it_init+1)/dble(slowsteps))
      else
         boostfactor=1.0d0+fld_boost_dt
      endif
      dtdifflimit=dtdifflimit*boostfactor
      if(fld_debug) print *,'DT limit from diffusion is boosted to=',dtdifflimit
      dtnew=min(dtnew,dtdifflimit)
      if(fld_debug) print *,'DT limit after diffusion is =',dtnew
    endif
 
    ! here we interface back to fld_get_radpress
    call phys_get_csrad2(w,x,ixi^l,ixo^l,tmp)
    if(slab_uniform) then
       ^d&dxinv(^d)=one/dx^d;
       do idims=1,ndim
          cmax(ixo^s)=dabs(w(ixo^s,iw_mom(idims)))+dsqrt(tmp(ixo^s))
          if(idims==1) then
            cmaxtot(ixo^s)=cmax(ixo^s)*dxinv(idims)
          else
            cmaxtot(ixo^s)=cmaxtot(ixo^s)+cmax(ixo^s)*dxinv(idims)
          end if
       end do
    else
       do idims=1,ndim
          cmax(ixo^s)=dabs(w(ixo^s,iw_mom(idims)))+dsqrt(tmp(ixo^s))
          if(idims==1) then
            cmaxtot(ixo^s)=cmax(ixo^s)/block%ds(ixo^s,idims)
          else
            cmaxtot(ixo^s)=cmaxtot(ixo^s)+cmax(ixo^s)/block%ds(ixo^s,idims)
          end if
       end do
    end if
    ! courantmaxtots='max(summed c/dx)'
    courantmaxtots=maxval(cmaxtot(ixo^s))
    if(fld_debug)print *,'DT limit RADFORCE CSRAD=',courantpar/courantmaxtots
    if(courantmaxtots>smalldouble) dtnew=min(dtnew,courantpar/courantmaxtots)
    if(fld_debug)print *,'DT limit FINALLY ENFORCED IS NOW=',dtnew
 
  end subroutine fld_radforce_get_dt
 
  !> Sets the opacity in the w-array
  !> by calling mod_opal_opacity
  !> NOTE: assumes primitives in w
  !> NOTE: assuming opacity is local, not ok with cak line force
  subroutine fld_get_opacity_prim(w, x, ixI^L, ixO^L, fld_kappa, fl)
    use mod_global_parameters
    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(ixi^s)
    type(fld_fluid), intent(in)  :: fl
 
    integer :: ix^d
    double precision :: rho0,temp0,kapp0
    double precision :: temp(ixi^s)
 
    select case (trim(fld_opacity_law))
      case('const_norm')
        fld_kappa(ixo^s) = fld_kappa0
      case('const')
        fld_kappa(ixo^s) = fld_kappa0/unit_opacity
      case('opal')
        call fl%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,kapp0)
          fld_kappa(ix^d) = kapp0/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_prim
 
  !> Returns Radiation Pressure as tensor
  !> NOTE: w is primitive on entry
  subroutine fld_get_radpress(w, x, ixI^L, ixO^L, rad_pressure, fl)
    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, intent(out):: rad_pressure(ixi^s,1:ndim,1:ndim)
    type(fld_fluid), intent(in)  :: fl
    integer :: i,j,nth
    double precision             :: eddington_tensor(ixi^s,1:ndim,1:ndim)
    double precision             :: lambda(ixi^s),fld_r(ixi^s)
 
    ! always use 4th order CD here
    nth=2
    call fld_get_eddington(w, x, ixi^l, ixo^l, eddington_tensor, lambda, fld_r, nth, fl)
    if(fld_debug.and..false.)then
       print *,'In get_radPress with nth=',nth,' on ixO=',ixo^l
       print *,'Max and Min value of fe'
       print *,maxval(eddington_tensor(ixo^s,1:ndim,1:ndim))
       print *,minval(eddington_tensor(ixo^s,1:ndim,1:ndim))
       print *,'Max and Min value of Erad'
       print *,maxval(w(ixo^s,iw_r_e))
       print *,minval(w(ixo^s,iw_r_e))
       print *,'End get_radPress'
    endif
    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
 
  !> This subroutine calculates flux limiter lambda according to fld_fluxlimiter
  !> It also calculates fld_R which is ratio of radiation scaleheight and mean free path
  !> NOTE: nth and ixI and ixO not free to choose here: TODO
  subroutine fld_get_fluxlimiter(w,x,ixI^L,ixO^L,fld_lambda,fld_R,nth,fl)
    use mod_global_parameters
    use mod_geometry
    use mod_usr_methods
    use mod_physics, only: phys_to_primitive
    integer, intent(in)           :: ixi^l,ixo^l,nth
    double precision, intent(in)  :: w(ixi^s,1:nw)
    double precision, intent(in)  :: x(ixi^s,1:ndim)
    double precision, intent(out) :: fld_r(ixi^s),fld_lambda(ixi^s)
    type(fld_fluid), intent(in)   :: fl
 
    double precision :: wprim(ixi^s,1:nw)
 
    wprim(ixi^s,1:nw)=w(ixi^s,1:nw)
    call phys_to_primitive(ixi^l,ixi^l,wprim,x)
    call fld_get_fluxlimiter_prim(wprim,x,ixi^l,ixo^l,fld_lambda,fld_r,nth,fl)
 
  end subroutine fld_get_fluxlimiter
 
  !> This subroutine calculates flux limiter lambda according to fld_fluxlimiter
  !> It also calculates fld_R which is ratio of radiation scaleheight and mean free path
  !> NOTE: this one operates on primitives
  !> NOTE: nth and ixI and ixO not free to choose 
  subroutine fld_get_fluxlimiter_prim(w,x,ixI^L,ixO^L,fld_lambda,fld_R,nth,fl)
    use mod_global_parameters
    use mod_geometry
    use mod_usr_methods
    integer, intent(in)           :: ixi^l,ixo^l,nth
    double precision, intent(in)  :: w(ixi^s,1:nw)
    double precision, intent(in)  :: x(ixi^s,1:ndim)
    double precision, intent(out) :: fld_r(ixi^s),fld_lambda(ixi^s)
    type(fld_fluid), intent(in)   :: fl
 
    integer :: idir, ix^d
    double precision :: kappa(ixi^s),normgrad2(ixi^s)
    double precision :: grad_r_e(ixi^s)
 
    select case(fld_fluxlimiter)
    case('Diffusion')
      ! optically thick limit
      fld_lambda(ixo^s) = 1.d0/3.d0
      fld_r(ixo^s) = zero
    case('FreeStream')
      ! optically thin limit
      normgrad2(ixo^s) = zero
      do idir=1,ndim
        call gradient(w(ixi^s,iw_r_e),ixi^l,ixo^l,idir,grad_r_e,nth)
        normgrad2(ixo^s) = normgrad2(ixo^s)+grad_r_e(ixo^s)**2
      end do
      call fld_get_opacity_prim(w,x,ixi^l,ixo^l,kappa,fl)
      ! Calculate R everywhere
      ! |grad E|/(rho kappa E)
      fld_r(ixo^s) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
      where(normgrad2(ixo^s)<smalldouble**2)
        ! treat uniform case as diffusion limit
        fld_r(ixo^s)=zero
        fld_lambda(ixo^s) = 1.0d0/3.0d0
      elsewhere
        fld_lambda(ixo^s) = one/fld_r(ixo^s)
      endwhere
    case('Pomraning')
      ! Calculate R everywhere
      ! |grad E|/(rho kappa E)
      normgrad2(ixo^s) = zero
      do idir = 1,ndim
        call gradient(w(ixi^s,iw_r_e),ixi^l,ixo^l,idir,grad_r_e,nth)
        normgrad2(ixo^s) = normgrad2(ixo^s) + grad_r_e(ixo^s)**2
      end do
      call fld_get_opacity_prim(w,x,ixi^l,ixo^l,kappa,fl)
      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)
    case('Minerbo')
      ! Calculate R everywhere
      ! |grad E|/(rho kappa E)
      normgrad2(ixo^s) = zero
      do idir = 1,ndim
        call gradient(w(ixi^s,iw_r_e),ixi^l,ixo^l,idir,grad_r_e,nth)
        normgrad2(ixo^s) = normgrad2(ixo^s) + grad_r_e(ixo^s)**2
      end do
      call fld_get_opacity_prim(w, x, ixi^l, ixo^l, kappa, fl)
      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))
        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
 
  end subroutine fld_get_fluxlimiter_prim
 
  !> Calculate Radiation Flux 
  !> NOTE: only for diagnostics purposes (w conservative on entry)
  !> This returns cell centered values for radiation flux 
  subroutine fld_get_radflux(w, x, ixI^L, ixO^L, rad_flux, fl)
    use mod_global_parameters
    use mod_geometry
    use mod_physics, only: phys_to_primitive
    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(ixi^s, 1:ndim)
    type(fld_fluid), intent(in)  :: fl
 
    integer :: idir,nth_for_fld
    double precision :: wprim(ixi^s,1:nw)
    double precision :: grad_r_e(ixi^s)
    double precision :: kappa(ixi^s), lambda(ixi^s), fld_r(ixi^s)
 
    wprim(ixi^s,1:nw)=w(ixi^s,1:nw)
    call phys_to_primitive(ixi^l,ixi^l,wprim,x)
 
    call fld_get_opacity_prim(wprim, x, ixi^l, ixo^l, kappa, fl)
    ! always use 4th order CD here
    nth_for_fld=2
    call fld_get_fluxlimiter_prim(wprim, x, ixi^l, ixo^l, lambda, fld_r, nth_for_fld, fl)
    !> Calculate the Flux using the fld closure relation
    !> F = -c*lambda/(kappa*rho) *grad E
    do idir = 1,ndim
      call gradient(wprim(ixi^s,iw_r_e),ixi^l,ixo^l,idir,grad_r_e,nth_for_fld)
      rad_flux(ixo^s,idir)=-(c_norm*lambda(ixo^s)/(kappa(ixo^s)*wprim(ixo^s,iw_rho)))*grad_r_e(ixo^s)
    end do
  end subroutine fld_get_radflux
 
  !> Calculate Eddington-tensor (where w is primitive)
  !> also feeds back the flux limiter lambda and R
  subroutine fld_get_eddington(w, x, ixI^L, ixO^L, eddington_tensor, lambda, fld_R, nth, fl)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in)          :: ixI^L, ixO^L, nth
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out) :: eddington_tensor(ixI^S,1:ndim,1:ndim)
    double precision, intent(out) :: lambda(ixI^S),fld_R(ixI^S)
    type(fld_fluid), intent(in)  :: fl
 
    integer :: idir,jdir
    double precision :: normgrad2(ixI^S)
    double precision :: tmp(ixI^S),grad_r_e(ixI^S,1:ndim)
    double precision :: nn_regularized(ixI^S,1:ndim,1:ndim)
 
    normgrad2(ixo^s) = zero
    do idir = 1,ndim
      call gradient(w(ixi^s, iw_r_e),ixi^l,ixo^l,idir,tmp,nth)
      grad_r_e(ixo^s,idir)=tmp(ixo^s)
      normgrad2(ixo^s)=normgrad2(ixo^s)+tmp(ixo^s)**2
    end do
    do idir = 1,ndim
      do jdir = 1,ndim
       if(idir==jdir)then
          nn_regularized(ixo^s,idir,jdir)=(grad_r_e(ixo^s,idir)*grad_r_e(ixo^s,jdir)+smalldouble**2)/(normgrad2(ixo^s)+smalldouble**2)
       else
          nn_regularized(ixo^s,idir,jdir)=(grad_r_e(ixo^s,idir)*grad_r_e(ixo^s,jdir))/(normgrad2(ixo^s)+smalldouble**2)
       endif
      enddo
    enddo
    ! get lambda and R
    call fld_get_fluxlimiter_prim(w,x,ixi^l,ixo^l,lambda,fld_r,nth,fl)
    ! store f_e= lambda + lambda^2 R^2
    tmp(ixo^s) = lambda(ixo^s)+(lambda(ixo^s)*fld_r(ixo^s))**2
    do idir = 1,ndim
      ! first compute the isotropic (diagonal) part
      eddington_tensor(ixo^s,idir,idir) = half*(one-tmp(ixo^s))
    enddo
    do idir = 1,ndim
      do jdir = 1,ndim
        ! initialize off-diagonal part here
        if(idir .ne. jdir) eddington_tensor(ixo^s,idir,jdir) = zero
        ! add part depending on unit vectors along gradient E
        eddington_tensor(ixo^s,idir,jdir) = eddington_tensor(ixo^s,idir,jdir)+&
            half*(3.d0*tmp(ixo^s)-one)*nn_regularized(ixo^s,idir,jdir)
      enddo
    enddo
  end subroutine fld_get_eddington
 
  subroutine get_and_check_egas_erad_from_conserved(w,ixI^L,ixO^L,e_gas,E_rad)
    use mod_global_parameters
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(out) :: e_gas(ixI^S),E_rad(ixI^S)
    integer :: ix^D
 
      !> e_gas is the INTERNAL ENERGY without KINETIC ENERGY
      e_gas(ixo^s) = w(ixo^s,iw_e)-half*sum(w(ixo^s,iw_mom(:))**2,dim=ndim+1)/w(ixo^s,iw_rho)
      if(allocated(iw_mag)) then
        e_gas(ixo^s) = e_gas(ixo^s)-half*sum(w(ixo^s,iw_mag(:))**2,dim=ndim+1)
      endif
      e_rad(ixo^s) = w(ixo^s,iw_r_e)
 
      if(check_small_values.and..not.fix_small_values)then
         {do ix^db= ixomin^db,ixomax^db\}
           if(e_gas(ix^d)<small_e.or.e_rad(ix^d)<small_r_e) then
              write(*,*) "Error in FLD get_egas_Erad: small value"
              write(*,*) "Iteration: ", it, " Time: ", global_time
              write(*,*) "Cell number: ", ix^d
              write(*,*) "internal energy density  is=",e_gas(ix^d)," versus   small_e=",small_e
              write(*,*) "radiation energy density is=",e_rad(ix^d)," versus small_r_e=",small_r_e
              call mpistop("FLD error:May need to turn on fixes")
           end if
         {end do\}
      endif
 
      if(fix_small_values)then
        {do ix^d = ixomin^d,ixomax^d\ }
          e_gas(ix^d) = max(e_gas(ix^d),small_e)
          e_rad(ix^d) = max(e_rad(ix^d),small_r_e)
        {enddo\}
      endif
 
  end subroutine get_and_check_egas_erad_from_conserved
 
  !> Calling all subroutines to perform the multigrid method
  !> Communicates rad_e and diff_coeff to multigrid library
  !> Advance psa=psb+dtfactor*qdt*F_im(psa)
  subroutine fld_implicit_update(dtfactor,qdt,qtC,psa,psb,fl)
    use mod_global_parameters
    use mod_forest
    use mod_multigrid_coupling
    use mod_input_output, only: get_global_maxima,get_global_minima
 
    type(state), target :: psa(max_blocks)
    type(state), target :: psb(max_blocks)
    double precision, intent(in) :: qdt
    double precision, intent(in) :: qtc
    double precision, intent(in) :: dtfactor
    type(fld_fluid), intent(in)  :: fl
 
    integer, parameter           :: max_its = 100
    integer                      :: n,ixo^l,ix^d
    double precision             :: res, max_residual, mg_lambda, fac
    double precision :: wmax(nw),wmin(nw)
    double precision :: wmaxb(nw),wminb(nw)
    integer :: iigrid, igrid
 
    if(fld_no_mg)return
    
    ! we need first to compute the (variable) diffusion coefficient on entire grid
    !   this must be done in mesh+1 ghostcell layer
    call update_diffcoeff(psa,fl)
 
    ! now we multiply the diffusion coefficient with dtfactor*dt on entire mesh+1 domain
    ixo^l=ixm^ll^ladd1;
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       {do ix^d = ixomin^d,ixomax^d\ }
          psa(igrid)%w(ix^d,i_diff_mg)= psa(igrid)%w(ix^d,i_diff_mg)*dtfactor*qdt
       {enddo\}
    end do
 
    fac = 1.0d0
    max_residual = fld_diff_tol
 
    if(fld_debug)then
       call get_global_maxima(wmax,psa)
       call get_global_minima(wmin,psa)
       call get_global_maxima(wmaxb,psb)
       call get_global_minima(wminb,psb)
       if(mype==0)then
          ! the MG needs to be scaled such that everything is order unity
          print *,'Currently at time=',global_time,' time step=',qdt,' dtfactor=',dtfactor
          print *,'at start of MG solver, we have fld_diff_tol       =',fld_diff_tol
          print *,'at start of MG solver, we have LHS E_rad range as :',wmax(iw_r_e),wmin(iw_r_e)
          print *,'at start of MG solver, we have Diff coeff range as:',wmax(i_diff_mg),wmin(i_diff_mg)
          print *,'at start of MG solver, we have density range as   :',wmax(iw_rho),wmin(iw_rho)
          print *,'at start of MG solver, we have RHS E_rad range as :',wmaxb(iw_r_e),wminb(iw_r_e)
          print *,'at start of MG solver, we have qdt as',qdt,' and max_residual=',max_residual
          print *,'at start of MG solver, ratio coeffs on level 1  =',wmax(i_diff_mg)/(dx(1,1)**2)
          print *,'at start of MG solver, ratio coeffs on level max=',wmax(i_diff_mg)/(dx(1,refine_max_level)**2)
       endif
    endif
 
    call mg_set_methods(mg)
    if(.not. mg%is_allocated) call mpistop("multigrid tree not allocated yet")
 
    ! Here we handle the global helmholtz problem with variable coefficient
    ! The equation we solve is div([D]^n nabla Erad^(n+1)) -(1/dt)Erad^(n+1)=-(1/dt)Erad^n
    ! we reformulate to div([D x dt]^n nabla Erad^(n+1)) -Erad^(n+1)=-Erad^n
    ! Helmholtz equation is div(eps nabla phi) -mg_lambda phi = f
    !    hence eps is our variable coefficient dt x D=fld_lambda*c/(kappa*rho)
    !    hence phi is Erad and mg_lambda is unity
    call vhelmholtz_set_lambda(fac)
    ! copy in the (variable) diffusion coefficient in mg_iveps
    ! NOTE: this copies also the ghostcell values for this coefficient to mg
    call mg_copy_to_tree(i_diff_mg, mg_iveps, factor=fac, state_from=psa)
    ! copy in the Erad variable in mg_iphi (the one we solve for)
    call mg_copy_to_tree(iw_r_e, mg_iphi, factor=fac, state_from=psa)
    ! copy in RHS f factor as Erad with factor -1
    call mg_copy_to_tree(iw_r_e, mg_irhs, factor=-fac, state_from=psb)
    ! becuase the variable coefficient is needed on lower grid levels in mg
    ! we need to restrict and adopt BCs for this variable: Neumann is set 
    call mg_restrict(mg, mg_iveps)
    call mg_fill_ghost_cells(mg, mg_iveps)
    ! Now try solving with MG
    call mg_fas_fmg(mg, .true., max_res=res)
    if(fld_debug.and.mype==0)print *,'MG residual obtained with FMG is =',res
    do n = 1, max_its
      if(res < max_residual) exit
      call mg_fas_vcycle(mg, max_res=res)
      if(fld_debug.and.mype==0)print *,'MG residual obtained with Vcycle =',res,' at step =',n
    end do
    if(fld_debug.and.mype==0)print *,'FINAL MG residual obtained is =',res
    if(res .ge. max_residual) then
      if (mype == 0) then 
        write(*,*) it, ' residual from MG ', res
        write(*,*) it, ' max_residual in MG ', max_residual
        write(*,*) it, ' dtfactor*qdt in MG ', qdt*dtfactor
        print *,'Currently at time=',global_time,' time step=',qdt,' dtfactor=',dtfactor
      endif
      if(fld_force_mg_converged) then
          call mpistop("no convergence in MG")
      else
          if(mype==0) write(*,*) 'WARNING for it=',it,' NO CONVERGENCE IN MG but still carry on'
      endif
    end if
    ! copy back the Erad variable in iw_r_e
    call mg_copy_from_tree_gc(mg_iphi, iw_r_e, state_to=psa)
 
   if(check_small_values.and..not.fix_small_values)then
       ixo^l=ixm^ll;
       do iigrid=1,igridstail; igrid=igrids(iigrid);
         {do ix^db= ixomin^db,ixomax^db\}
           if(psa(igrid)%w(ix^d,iw_r_e)<small_r_e) then
              write(*,*) "Error in FLD fld_implicit_update: small value"
              write(*,*) "of radiation energy density after MG"
              write(*,*) "Iteration: ", it, " Time: ", global_time
              write(*,*) "Location: ", psa(igrid)%x(ix^d,:)," on grid",igrid
              write(*,*) "Cell number: ", ix^d
              write(*,*) "radiation energy density is=",psa(igrid)%w(ix^d,iw_r_e)," versus small_r_e=",small_r_e
              call mpistop("FLD error:May need to turn on fixes")
           end if
         {end do\}
       end do
    endif
 
    if(fix_small_values)then
       ixo^l=ixm^ll;
       do iigrid=1,igridstail; igrid=igrids(iigrid);
         {do ix^d = ixomin^d,ixomax^d\ }
           psa(igrid)%w(ix^d,iw_r_e)= max(psa(igrid)%w(ix^d,iw_r_e),small_r_e)
         {enddo\}
       end do
    endif
 
  end subroutine fld_implicit_update
 
  subroutine update_diffcoeff(psa,fl)
    use mod_global_parameters
    type(state), target :: psa(max_blocks)
    type(fld_fluid), intent(in) :: fl
 
    integer :: iigrid, igrid, ixO^L
 
    ! we will need diffusion coefficients in i-1 i i+1
    ixo^l=ixm^ll^ladd1;
    !$OMP PARALLEL DO PRIVATE(igrid)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
        call fld_get_diffcoef_central(psa(igrid)%w, psa(igrid)%x, ixg^ll, ixo^l, fl)
    end do
    !$OMP END PARALLEL DO
 
  end subroutine update_diffcoeff
 
  !> inplace update of psa==>F_im(psa)
  subroutine fld_evaluate_implicit(qtC,psa,fl)
    use mod_global_parameters
    type(state), target :: psa(max_blocks)
    double precision, intent(in) :: qtc
    type(fld_fluid), intent(in)  :: fl
    integer :: iigrid, igrid
    integer :: ixo^l
 
    if(fld_no_mg)return
 
    call update_diffcoeff(psa,fl)
 
    ixo^l=ixm^ll;
    !$OMP PARALLEL DO PRIVATE(igrid)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
       call evaluate_diffterm_onegrid(ixg^ll,ixo^l,psa(igrid)%w,psa(igrid)%x)
    end do
    !$OMP END PARALLEL DO
 
  end subroutine fld_evaluate_implicit
 
  !> inplace update of psa==>F_im(psa)
  subroutine evaluate_diffterm_onegrid(ixI^L,ixO^L,w,x)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(inout) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
 
    double precision :: divF(ixI^S)
    integer :: idir, jxO^L, hxO^L
 
    if(.not.slab)call mpistop("laplacian coded up for uniform cartesian grid")
 
    ! Here we use diffusion coefficients in positions i-1 i i+1 and exploit harmonic means i.e. 2ab/(a+b)
    ! since this is how the multigrid library handles the div(eps nabla phi) term in Helmholtz equation
    !   div(eps nabla phi) - mg_lambda phi = f
    divf(ixo^s) = 0.d0
    do idir = 1,ndim
       hxo^l=ixo^l-kr(idir,^d);
       jxo^l=ixo^l+kr(idir,^d);
       divf(ixo^s) = divf(ixo^s) + &
           (w(jxo^s,iw_r_e)*two*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(ixo^s,iw_r_e)*(two*w(ixo^s,i_diff_mg)*w(jxo^s,i_diff_mg)/(w(ixo^s,i_diff_mg) + w(jxo^s,i_diff_mg)) &
                           +two*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)*two*w(ixo^s,i_diff_mg)*w(hxo^s,i_diff_mg)/(w(ixo^s,i_diff_mg) + w(hxo^s,i_diff_mg)))/dxlevel(idir)**2
       ! below uses artihmetic mean, different from mg method
       !divF(ixO^S) = divF(ixO^S) + &
       !    (w(jxO^S,iw_r_e)*half*(w(ixO^S,i_diff_mg) + w(jxO^S,i_diff_mg)) &
       !   -w(ixO^S,iw_r_e)*(half*(w(ixO^S,i_diff_mg) + w(jxO^S,i_diff_mg)) &
       !                    +half*(w(ixO^S,i_diff_mg) + w(hxO^S,i_diff_mg))) &
       !    +w(hxO^S,iw_r_e)*half*(w(ixO^S,i_diff_mg) + w(hxO^S,i_diff_mg)))/dxlevel(idir)**2
    enddo
    ! only the E variable is handled implicitly, all else must be zero here
    w(ixo^s,1:nw)=zero
    w(ixo^s,iw_r_e) = divf(ixo^s)
 
  end subroutine evaluate_diffterm_onegrid
 
  !> Calculates cell-centered diffusion coefficient to be used in multigrid
  subroutine fld_get_diffcoef_central(w, x, ixI^L, ixO^L, fl)
    use mod_global_parameters
    use mod_geometry
    use mod_usr_methods
    use mod_physics, only: phys_to_primitive
    integer, intent(in)          :: ixi^l, ixo^l
    double precision, intent(in) :: x(ixi^s,1:ndim)
    double precision, intent(inout) :: w(ixi^s,1:nw)
    type(fld_fluid), intent(in)  :: fl
 
    double precision :: wprim(ixi^s,1:nw)
    double precision :: kappa(ixi^s),lambda(ixi^s),fld_r(ixi^s)
 
    wprim(ixi^s,1:nw)=w(ixi^s,1:nw)
    ! ensure entries in entire ixI range
    call phys_to_primitive(ixi^l,ixi^l,wprim,x)
    call fld_get_opacity_prim(wprim, x, ixi^l, ixo^l, kappa, fl)
    ! note that we use central difference here (last argument is 1 or 2)
    call fld_get_fluxlimiter_prim(wprim, x, ixi^l, ixo^l, lambda, fld_r, nth_for_diff_mg, fl)
    w(ixo^s,i_diff_mg) = c_norm*lambda(ixo^s)/(kappa(ixo^s)*wprim(ixo^s,iw_rho))
    if(associated(usr_special_diffcoef)) call usr_special_diffcoef(w, wprim, x, ixi^l, ixo^l)
    if(minval(w(ixo^s,i_diff_mg))<smalldouble) then
        print *,'min diffcoef=',minval(w(ixo^s,i_diff_mg))
        call mpistop("too small diffusion coefficient")
    endif
    if(maxval(w(ixo^s,i_diff_mg))>bigdouble)  call mpistop("too large diffusion coefficient")
 
    if(fld_debug.and..false.)then
      print *,'setting diffcoefs with data on',ixi^l
      print *,'min diffcoef=',minval(w(ixo^s,i_diff_mg))
      print *,'min lambda kappa rho fld_R'
      print *,minval(lambda(ixo^s))
      print *,minval(kappa(ixo^s))
      print *,minval(wprim(ixo^s,iw_rho))
      print *,minval(fld_r(ixo^s))
      print *,'max diffcoef=',maxval(w(ixo^s,i_diff_mg))
      print *,'max lambda kappa rho fld_R'
      print *,maxval(lambda(ixo^s))
      print *,maxval(kappa(ixo^s))
      print *,maxval(wprim(ixo^s,iw_rho))
      print *,maxval(fld_r(ixo^s))
      print *,'done setting diffcoefs in slot',i_diff_mg,' on range',ixo^l
    endif
 
  end subroutine fld_get_diffcoef_central
 
 
  !> Find the root of the 4th degree polynomial using the bisection method
  subroutine bisection_method(e_gas, c0, c1)
    use mod_global_parameters
    double precision, intent(in)    :: c0, c1
    double precision, intent(inout) :: e_gas
    double precision :: bisect_a, bisect_b, bisect_c
    integer :: n, max_its
 
    n = 0
    max_its = 100
    bisect_a = zero
    bisect_b = min(dabs(c0/c1),dabs(c0)**(1.d0/4.d0))+smalldouble
    do while(dabs(bisect_b-bisect_a) .ge. fld_bisect_tol)
      bisect_c = (bisect_a + bisect_b)/two
      n = n +1
      if(n .gt. max_its) then
        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
        if(fld_debug)print*, "IGNORING GAS-RAD ENERGY EXCHANGE ", c0, c1
        if(fld_debug)print*, polynomial_bisection(bisect_a, c0, c1), polynomial_bisection(bisect_b, c0, c1)
        call mpistop('issues in bisection scheme')
        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 method
  subroutine newton_method(e_gas, c0, c1)
    use mod_global_parameters
    double precision, intent(in)    :: c0, c1
    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(dabs(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
        if(fld_debug)print*, 'skip to bisection algorithm'
        call bisection_method(e_gas, 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, c0, c1)
    use mod_global_parameters
    double precision, intent(in)    :: c0, c1
    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 (dabs(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
        if(fld_debug)print*, 'skip to Newton algorithm'
        call newton_method(e_gas, 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
end module mod_fld