Zav_of_te Function

  function Zav_of_te(imp_element_index,te)

    !! Electron temperature dependent average atomic number
    !! author: H Lux, CCFE, Culham Science Centre
    !! imp_element : input imp_dat : impurity element
    !! te : input real : electron temperature (keV)
    !! This routine returns the interpolated average atomic
    !! charge for a given electron temperature.
    !! <P>The Zav versus temperature data is interpolated from
    !! lookup tables from the ADAS data base provided by Martin O'Mullane.
    !! None
    !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       implicit none

    real(dp) :: Zav_of_te

    !  Arguments

    integer,  intent(in) :: imp_element_index
    real(dp), intent(in) :: te

    !  Local variables

    integer :: i
    real(dp) :: xi, yi, c

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !  Interpolate tabulated data

    !  Temperatures lower than the minimum, or higher than the maximum,
    !  are dealt with by taking the Zav value at the nearest tabulated
    !  temperature point

    if (te <= impurity_arr_Temp_keV(imp_element_index, 1)) then
       ! This should not be too unreasonable.
       Zav_of_te = impurity_arr_Zav(imp_element_index,1)

    else if (te >= impurity_arr_Temp_keV(imp_element_index, impurity_arr_len_tab(imp_element_index))) then
       !  This should be okay, as most elements are fully ionised by now.
       Zav_of_te = impurity_arr_Zav(imp_element_index, impurity_arr_len_tab(imp_element_index))

    else

       do i = 1, impurity_arr_len_tab(imp_element_index)-1

          !  Linear interpolation in log-lin space

          if ( (te > impurity_arr_Temp_keV(imp_element_index, i)) .and. &
               (te <= impurity_arr_Temp_keV(imp_element_index, i+1)) ) then

             yi = impurity_arr_Zav(imp_element_index, i)
             xi = log(impurity_arr_Temp_keV(imp_element_index, i))
             c  = (impurity_arr_Zav(imp_element_index, i+1) - yi) / &
                  (log(impurity_arr_Temp_keV(imp_element_index,i+1)) - xi)
             Zav_of_te = yi + c * (log(te) - xi)
             exit
          end if

       end do

    end if

  end function Zav_of_te