analytical_profiles_p Subroutine

  subroutine analytical_profiles_p(time,params,Y_R,Y_Z,P,F,ne,Te,Zeff,PSIp)
    !! @note Subroutine that calculates the analytical plasma profiles at
    !! the particles' position. @endnote
    TYPE(KORC_PARAMS), INTENT(IN)                           :: params
    REAL(rp), DIMENSION(8), INTENT(IN)  :: Y_R,Y_Z,PSIp
    REAL(rp), INTENT(IN)  :: time
    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]].
    TYPE(FIELDS), INTENT(IN)      :: F
    REAL(rp), DIMENSION(8),INTENT(OUT) :: ne
    !! Background electron density seen by simulated particles.
    REAL(rp), DIMENSION(8),INTENT(OUT) :: Te
    !! Backgroun temperature density seen by simulated particles.
    REAL(rp), DIMENSION(8),INTENT(OUT) :: Zeff
    !! Effective atomic charge seen by simulated particles.
    INTEGER(ip)                                        :: cc
    !! Particle iterator.
    REAL(rp) :: R0,Z0,a,ne0,n_ne,Te0,n_Te,Zeff0
    REAL(rp) :: R0_RE,Z0_RE,sigmaR_RE,sigmaZ_RE,psimax_RE
    REAL(rp) :: n_REr0,n_tauion,n_lamfront,n_lamback,n_lamshelf
    REAL(rp) :: n_psifront,n_psiback,n_psishelf
    REAL(rp) :: n_tauin,n_tauout,n_shelfdelay,n_shelf
    REAL(rp) :: n0t,n_taut
    REAL(rp) :: PSIp0,PSIp_lim,psiN_0
    REAL(rp), DIMENSION(8) :: r_a,rm,rm_RE,PSIpN,PSIp_temp

    R0=P%R0
    Z0=P%Z0
    a=P%a
    
    ne0=P%neo
    n_ne=P%n_ne

    Te0=P%Teo
    n_Te=P%n_Te

    Zeff0=P%Zeffo

    R0_RE=P%R0_RE
    Z0_RE=P%Z0_RE
    n_REr0=P%n_REr0
    n_tauion=P%n_tauion
    n_tauin=P%n_tauin
    n_tauout=P%n_tauout
    n_shelfdelay=P%n_shelfdelay
    n_lamfront=P%n_lamfront
    n_lamback=P%n_lamback
    n_lamshelf=P%n_lamshelf
    n_psifront=P%n_lamfront*params%cpp%length
    n_psiback=P%n_lamback*params%cpp%length
    n_psishelf=P%n_lamshelf*params%cpp%length
    n_shelf=P%n_shelf
    
    PSIp_lim=F%PSIp_lim
    PSIp0=F%PSIP_min
    psiN_0=P%psiN_0
    

!    write(6,*) 'PSIp',PSIp(1)*(params%cpp%Bo*params%cpp%length**2)
!    write(6,*) 'PSIp_lim',PSIp_lim*(params%cpp%Bo*params%cpp%length**2)   
!    write(6,*) 'PSIp0',PSIp0*(params%cpp%Bo*params%cpp%length**2)
    
!    write(6,'("R0_RE: "E17.10)') R0_RE
!    write(6,'("Z0_RE: "E17.10)') Z0_RE
!    write(6,'("n_REr0: "E17.10)') n_REr0

    
    SELECT CASE (TRIM(P%ne_profile))
    CASE('FLAT')
       !$OMP SIMD
       do cc=1_idef,8_idef
          ne(cc) = ne0
       end do
       !$OMP END SIMD

    CASE('SPONG')
       !$OMP SIMD
       do cc=1_idef,8_idef
          rm(cc)=sqrt((Y_R(cc)-R0)**2+(Y_Z(cc)-Z0)**2)
          r_a(cc)=rm(cc)/a
          ne(cc) = ne0*(1._rp-0.2*r_a(cc)**8)+n_ne
       end do
       !$OMP END SIMD
    CASE('RE-EVO')
       !$OMP SIMD
       do cc=1_idef,8_idef
          rm_RE(cc)=sqrt((Y_R(cc)-R0_RE)**2+(Y_Z(cc)-Z0_RE)**2)
          ne(cc) = (ne0-n_ne)/4._rp*(1+tanh((rm_RE(cc)+ &
               n_REr0*(time/n_tauion-1))/n_lamfront))* &
               (1+tanh(-(rm_RE(cc)-n_REr0)/n_lamback))+n_ne
       end do
       !$OMP END SIMD
    CASE('RE-EVO1')
       !$OMP SIMD
       do cc=1_idef,8_idef
          rm_RE(cc)=sqrt((Y_R(cc)-R0_RE)**2+(Y_Z(cc)-Z0_RE)**2)
          ne(cc) = (ne0-n_ne)/8._rp*(1+tanh((rm_RE(cc)+ &
               n_REr0*(time/n_tauion-1))/n_lamfront))* &
               (1+tanh(-(rm_RE(cc)-n_REr0)/n_lamback))* &
               (2*(n_shelf-n_ne)/(ne0-n_ne)+(ne0-n_shelf)/(ne0-n_ne)* &
               (1-tanh((rm_RE(cc)+n_REr0*((time-n_shelfdelay)/n_tauin-1))/ &
               n_lamshelf)))+n_ne
       end do
       !$OMP END SIMD

       
    CASE('RE-EVO-PSI')
       !$OMP SIMD
       do cc=1_idef,8_idef
          PSIpN(cc)=(PSIp(cc)-PSIp0)/(PSIp_lim-PSIp0)
          ne(cc) = (ne0-n_ne)/8._rp*(1+tanh((sqrt(abs(PSIpN(cc)))+ &
               sqrt(abs(psiN_0))*(time/n_tauion-1))/n_psifront))* &
               (1+tanh(-(sqrt(abs(PSIpN(cc)))-sqrt(abs(psiN_0)))/n_psiback))* &
               (2*(n_shelf-n_ne)/(ne0-n_ne)+(ne0-n_shelf)/(ne0-n_ne)* &
               (1-tanh((sqrt(abs(PSIpN(cc)))+ sqrt(abs(psiN_0))* &
               ((time-n_shelfdelay)/n_tauin-1))/n_psishelf)))+n_ne
       end do
       !$OMP END SIMD

       !       write(6,*) 'at time ',time*params%cpp%time, &
       !' ne: ',ne(1)/params%cpp%length**3
      
!       !$OMP SIMD
!       do cc=1_idef,8
!          if(isnan(ne(cc))) then
!             write(6,*) 'PSIp: ',PSIp(cc)
!             write(6,*) 'PSIp0: ',PSIp0
!             write(6,*) 'PSIp_lim: ',PSIp_lim
!             write(6,*) 'PSIpN: ',PSIpN(cc)

!             stop 'ne_eval is a NaN'
!          end if
!       end do
       !       !$OMP END SIMD

    CASE('RE-EVO-PSIN-SG')

       n0t=(ne0-n_ne)/2._rp*(tanh(time/n_tauin)- &
            tanh((time-n_shelfdelay)/n_tauin))
       n_taut=n_psishelf*erf((time+params%dt/100._rp)/n_tauion)
       
       !$OMP SIMD
       do cc=1_idef,8_idef
          PSIpN(cc)=(PSIp(cc)-PSIp0)/(PSIp_lim-PSIp0)
          ne(cc) = n0t*exp(-(sqrt(abs(PSIpN(cc)))-sqrt(abs(psiN_0)))**2._rp/ &
               (2._rp*n_taut**2._rp))*(1._rp+erf(-10._rp* &
               (sqrt(abs(PSIpN(cc)))-sqrt(abs(psiN_0)))/ &
               (sqrt(2._rp)*n_taut)))/2._rp+n_ne
       end do
       !$OMP END SIMD
       
    CASE('RE-EVO-PSIP-G')

!       write(6,*) 'time: ',time*params%cpp%time
       
       n0t=(ne0-n_ne)/2._rp*(tanh((time-n_tauin)/n_tauin)- &
            tanh((time-n_shelfdelay)/n_tauout))
       n_taut=n_psishelf*erf((time+params%dt/100._rp)/n_tauion)
       
       !$OMP SIMD
       do cc=1_idef,8_idef
          PSIp_temp(cc)=PSIp(cc)*(params%cpp%Bo*params%cpp%length**2)
          ne(cc) = n0t*exp(-(sqrt(abs(PSIp_temp(cc)))-sqrt(abs(psiN_0)))**2._rp/ &
               (2._rp*n_taut**2._rp))+n_ne
       end do
       !$OMP END SIMD
       
    CASE DEFAULT
       !$OMP SIMD
       do cc=1_idef,8_idef
          ne(cc) = ne0
       end do
       !$OMP END SIMD
    END SELECT

    SELECT CASE (TRIM(P%Te_profile))
    CASE('FLAT')
       !$OMP SIMD
       do cc=1_idef,8_idef
          Te(cc) = Te0
       end do
       !$OMP END SIMD
    CASE('SPONG')
       !$OMP SIMD
       do cc=1_idef,8_idef
          rm(cc)=sqrt((Y_R(cc)-R0)**2+(Y_Z(cc)-Z0)**2)
          r_a(cc)=rm(cc)/a
          Te(cc) = Te0*(1._rp-0.6*r_a(cc)**2)**2+Te0*n_Te
       end do
       !$OMP END SIMD
    CASE DEFAULT
       !$OMP SIMD
       do cc=1_idef,8_idef
          Te(cc) = P%Teo
       end do
       !$OMP END SIMD
    END SELECT

    SELECT CASE (TRIM(P%Zeff_profile))
    CASE('FLAT')
       !$OMP SIMD
       do cc=1_idef,8_idef
          Zeff(cc) = P%Zeffo
       end do
       !$OMP END SIMD
    CASE('SPONG')
       !$OMP SIMD
       do cc=1_idef,8_idef
          Zeff(cc) = P%Zeffo
       end do
       !$OMP END SIMD
    CASE DEFAULT
       !$OMP SIMD
       do cc=1_idef,8_idef
          Zeff(cc) = P%Zeffo
       end do
       !$OMP END SIMD
    END SELECT
          

!    write(6,*) PSIpN(1)
    
!    write(6,'("ne: "E17.10)') ne(1)/params%cpp%length**3
!    write(6,'("rm_RE: "E17.10)') rm_RE(1)
    
  end subroutine analytical_profiles_p