program minimal_example
  use, intrinsic :: IEEE_arithmetic

  integer, parameter :: pReal = selected_real_kind(15,307)
  real, dimension(120,5) :: state = 123.0_pReal
  integer, dimension(*), parameter :: mob = [1,2,3,4]
  integer, dimension(*), parameter :: dip = [9,10]

  type :: tParameters
    real(pReal) :: &
      rho_min, &
      rho_significant
    integer :: &
      sum_N_sl
  end type tParameters

  type(tParameters), dimension(1) :: param

  param(1)%sum_N_sl = 12
  param(1)%rho_min = 1.0_pReal
  param(1)%rho_significant = 1.0_pReal

  print*, getRho(1,2)
  
  contains
  function getRho(ph,en) result(getRho_out)
  
    
    integer, intent(in) :: ph, en
    real(pReal), dimension(param(ph)%sum_N_sl,10) :: getRho_out
  
  
    associate(prm => param(ph))
  
      getRho_out = reshape(state(:,en),[prm%sum_N_sl,10])
  
      ! ensure positive densities (not for imm, they have a sign)
      getRho_out(:,mob) = max(getRho_out(:,mob),0.0_pReal)
      getRho_out(:,dip) = max(getRho_out(:,dip),0.0_pReal)
  
      where(abs(getRho_out) < max(prm%rho_min/1e-6_pReal**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
        getRho_out = 0.0_pReal
  
    end associate
  
  end function getRho

end program minimal_example
