mod_ne_te_profile.f90

Go to the documentation of this file.

module mod_ne_te_profile

  use itm_types
  use mod_dete
  use mod_helena_io, only : out_he

  implicit none

  contains

!---------------------------------------------------------------------------
  subroutine ne_te_profile(dt, psi, f, df, const, df_in)
!---------------------------------------------------------------------------
! calculates the values of T_e and n_i and the derivative in psi
!---------------------------------------------------------------------------

    implicit none

    type(ne_te), intent(in) :: dt
    real(r8), intent(in) :: psi
    real(r8), intent(out) :: f, df
    real(r8), intent(out), optional :: const
    real(r8), intent(in), optional :: df_in

    real(r8) :: ex, ext, denom, ztmp

    if (.not. present(df_in)) then

      ex = exp(psi / dt%delta)
      ext = exp(dt%psi_m / dt%delta)
      denom = 2.0_r8 * dt%delta * (ex + ext)
      ztmp = (dt%psi_m - psi) / (2.0_r8 * dt%delta)

      f = (0.5_r8 * (dt%h - dt%h_0) * (((1.0_r8 + dt%slope * ztmp) &
       * exp(ztmp) - exp(-ztmp)) / (exp(ztmp) + exp(-ztmp)) + 1.0_r8) &
       + dt%h_0) * dt%fac

      df = (dt%h - dt%h_0) * ext * (dt%delta * ((4.0_r8 &
       + dt%slope) * ex + dt%slope * ext) - dt%slope *ex &
       * (psi - dt%psi_m)) / (denom * denom)

      if (.not. present(const)) return

      if (dt%shape == 1) then
        const = f + dt%fac * df * (dt%alpha(1) + (1._r8 - dt%alpha(1)) &
         / (dt%alpha(2) + 1.0_r8)) * psi
      else
        const = f + dt%fac * df * ((2.0_r8 / 3.0_r8 * dt%alpha(1) &
         + 0.5_r8 * dt%alpha(2)) / (dt%alpha(1) + dt%alpha(2)) &
         + 0.5_r8 * dt%alpha(3))
      end if

    else

      if (dt%shape == 1) then
!       df = -dt%fac * df_in * (dt%alpha(1) + (1.0_r8 - dt%alpha(1)) &
!       * (psi &
        df = df_in * (dt%alpha(1) + (1.0_r8 - dt%alpha(1)) * (psi &
         / dt%psi_0) **dt%alpha(2))
        f = -dt%fac * df_in * (dt%alpha(1) * psi + (1.0_r8 - dt%alpha(1)) &
         / (dt%alpha(2) + 1.0_r8) / dt%psi_0**dt%alpha(2) &
         * psi**(dt%alpha(2) + 1.0_r8)) + const
      else
!       df = -dt%fac * df_in / dt%psi_0 * ((dt%alpha(1) * (psi &
!        / dt%psi_0)**0.5_r8 &
        df = df_in / dt%psi_0 * ((dt%alpha(1) * (psi / dt%psi_0)**0.5_r8 &
         + dt%alpha(2) * (psi / dt%psi_0)) / (dt%alpha(1) + dt%alpha(2)) & 
         + dt%alpha(3) * (1.0_r8 - psi / dt%psi_0))
        f = -dt%fac * df_in * ((2.0_r8 / 3.0_r8 * dt%alpha(1) * (psi &
         / dt%psi_0) **1.5_r8 + 0.5_r8 * dt%alpha(2) * (psi &
         / dt%psi_0)**2) / (dt%alpha(1) + dt%alpha(2)) + dt%alpha(3) * (psi &
         / dt%psi_0 - 0.5_r8 * (psi / dt%psi_0)**2)) + const
      end if

    end if
    return
  end subroutine ne_te_profile   

end module mod_ne_te_profile