get_analytical_profiles Subroutine

  subroutine get_analytical_profiles(P,Y,ne,Te,Zeff,flag)
    !! @note Subroutine that calculates the analytical plasma profiles at
    !! the particles' position. @endnote
    TYPE(PROFILES), INTENT(IN)                         :: P
    !! An instance of KORC's derived type PROFILES containing all the
    !! information about the plasma profiles used in the simulation.
    !! See [[korc_types]] and [[korc_profiles]].
    REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN)  :: Y
    !! Particles' position in toroidal coordinates; Y(1,:) = \(r\), Y(2,:)
    !! = \(\theta\), Y(3,:) = \(\zeta\).
    REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: ne
    !! Background electron density seen by simulated particles.
    REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: Te
    !! Backgroun temperature density seen by simulated particles.
    REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: Zeff
    !! Effective atomic charge seen by simulated particles.
    INTEGER(is), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: flag
    !! Flag for each particle to decide whether it is being followed
    !! (flag=T) or not (flag=F).
    REAL(rp)                                           :: r_a
    !! Normalized toroidal radial position of simulated particles
    !! \(r/a\), where \(a\) is the plasma radius.
    REAL(rp)                                           :: fr
    !! Calculated radial profile.
    INTEGER(ip)                                        :: pp
    !! Particle iterator.
    INTEGER(ip)                                        :: ss
    !! Species iterator.

    if (Y(2,1).eq.0) then
       ss=1_idef
    else
       ss = size(Y,1)
    end if
   

    !$OMP PARALLEL DO FIRSTPRIVATE(ss) PRIVATE(pp,fr,r_a) &
    !$OMP& SHARED(P,Y,ne,Te,Zeff,flag)
    do pp=1_idef,ss
       if ( flag(pp) .EQ. 1_is ) then

          r_a = Y(pp,1)/P%a

!          write(output_unit_write,'("r: ",E17.10)') r_a
          
          SELECT CASE (TRIM(P%ne_profile))
          CASE('TANH')
             fr = 1_ip - TANH(2.0_rp*r_a)**P%n_ne
             ne(pp) = P%neo*fr
          CASE('FLAT')
             ne(pp) = P%neo
          CASE('POLYNOMIAL')
             fr = P%a_ne(1) + P%a_ne(2)*r_a + P%a_ne(3)*r_a**2 + &
                  P%a_ne(4)*r_a**3
             ne(pp) = P%neo*fr
          CASE('SPONG')
             fr = P%neo*(1._rp-0.2*r_a**8)+P%neo*P%n_ne
             ne(pp) = P%neo*fr
          CASE DEFAULT
             ne(pp) = P%neo
          END SELECT

          SELECT CASE (TRIM(P%Te_profile))
          CASE('TANH')
             fr = 1_ip - TANH(2.0_rp*r_a)**P%n_Te
             Te(pp) = P%Teo*fr
          CASE('FLAT')
             Te(pp) = P%Teo
          CASE('POLYNOMIAL')
             fr = P%a_Te(1) + P%a_Te(2)*r_a + P%a_Te(3)*r_a**2 + &
                  P%a_Te(4)*r_a**3
             Te(pp) = P%Teo*fr
          CASE('SPONG')
             fr = P%Teo*(1._rp-0.6*r_a**2)**2+P%Teo*P%n_Te
             ne(pp) = P%neo*fr
          CASE DEFAULT
             Te(pp) = P%Teo
          END SELECT

          SELECT CASE (TRIM(P%Zeff_profile))
          CASE('TANH')
             fr = 1_ip - TANH(2.0_rp*r_a)**P%n_Zeff
             Zeff(pp) = P%Zeffo*fr
          CASE('FLAT')
             Zeff(pp) = P%Zeffo
          CASE('POLYNOMIAL')
             fr = P%a_Zeff(1) + P%a_Zeff(2)*r_a + P%a_Zeff(3)*r_a**2 + &
                  P%a_Zeff(4)*r_a**3
             Zeff(pp) = P%Zeffo*fr
          CASE('SPONG')
             Zeff(pp) = P%Zeffo
          CASE DEFAULT
             Zeff(pp) = P%Zeffo
          END SELECT
       end if
    end do
    !$OMP END PARALLEL DO
  end subroutine get_analytical_profiles