Subroutine that calculates the analytical plasma profiles at the particles' position.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | pchunk | |||
| real(kind=rp), | intent(in) | :: | time | |||
| type(KORC_PARAMS), | intent(in) | :: | params | |||
| real(kind=rp), | intent(in), | DIMENSION(params%pchunk) | :: | Y_R | ||
| real(kind=rp), | intent(in), | DIMENSION(params%pchunk) | :: | Y_Z | ||
| 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(kind=rp), | intent(out), | DIMENSION(params%pchunk) | :: | ne | Background electron density seen by simulated particles. |
|
| real(kind=rp), | intent(out), | DIMENSION(params%pchunk) | :: | Te | Backgroun temperature density seen by simulated particles. |
|
| real(kind=rp), | intent(out), | DIMENSION(params%pchunk) | :: | Zeff | Effective atomic charge seen by simulated particles. |
|
| real(kind=rp), | intent(in), | DIMENSION(params%pchunk) | :: | PSIp |
subroutine analytical_profiles_p(pchunk,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
INTEGER, INTENT(IN) :: pchunk
REAL(rp), DIMENSION(params%pchunk), 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(params%pchunk),INTENT(OUT) :: ne
!! Background electron density seen by simulated particles.
REAL(rp), DIMENSION(params%pchunk),INTENT(OUT) :: Te
!! Backgroun temperature density seen by simulated particles.
REAL(rp), DIMENSION(params%pchunk),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,R0a
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(params%pchunk) :: r_a,rm,rm_RE,PSIpN,PSIp_temp
R0=P%R0
Z0=P%Z0
a=P%a
R0a=F%AB%Ro
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(output_unit_write,*) 'PSIp',PSIp(1)*(params%cpp%Bo*params%cpp%length**2)
! write(output_unit_write,*) 'PSIp_lim',PSIp_lim*(params%cpp%Bo*params%cpp%length**2)
! write(output_unit_write,*) 'PSIp0',PSIp0*(params%cpp%Bo*params%cpp%length**2)
! write(output_unit_write,'("R0_RE: "E17.10)') R0_RE
! write(output_unit_write,'("Z0_RE: "E17.10)') Z0_RE
! write(output_unit_write,'("n_REr0: "E17.10)') n_REr0
SELECT CASE (TRIM(P%ne_profile))
CASE('FLAT')
!$OMP SIMD
do cc=1_idef,pchunk
ne(cc) = ne0
end do
!$OMP END SIMD
CASE('FLAT-RAMP')
!$OMP SIMD
do cc=1_idef,pchunk
ne(cc) = n_ne+(ne0-n_ne)*time/n_tauion
end do
!$OMP END SIMD
CASE('TANH-RAMP')
!$OMP SIMD
do cc=1_idef,pchunk
ne(cc) = n_ne+(ne0-n_ne)/2*(tanh((time-n_shelfdelay)/n_tauion)+1._rp)
end do
!$OMP END SIMD
CASE('SINE')
!$OMP SIMD
do cc=1_idef,pchunk
ne(cc) = n_ne+(ne0-n_ne)*sin(time/n_tauion)
end do
!$OMP END SIMD
CASE('SPONG')
!$OMP SIMD
do cc=1_idef,pchunk
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('MST_FSA')
!$OMP SIMD
do cc=1_idef,pchunk
rm(cc)=sqrt((Y_R(cc)-R0a)**2+(Y_Z(cc)-Z0)**2)
r_a(cc)=rm(cc)/a
ne(cc) = (ne0-n_ne)*(1._rp-r_a(cc)**4._rp)**4._rp+n_ne
!write(6,*) 'R',Y_R*params%cpp%length,'R0',R0*params%cpp%length, &
! 'Z',Y_Z*params%cpp%length,'Z0',Z0*params%cpp%length, &
! 'a',a*params%cpp%length
!write(6,*) 'r_a',r_a,'ne',ne(cc)*params%cpp%density
end do
!$OMP END SIMD
CASE('RE-EVO')
!$OMP SIMD
do cc=1_idef,pchunk
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,pchunk
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,pchunk
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(output_unit_write,*) '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(output_unit_write,*) 'PSIp: ',PSIp(cc)
! write(output_unit_write,*) 'PSIp0: ',PSIp0
! write(output_unit_write,*) 'PSIp_lim: ',PSIp_lim
! write(output_unit_write,*) '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,pchunk
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(output_unit_write,*) '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,pchunk
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('RE-EVO-PSIP-G1')
! write(output_unit_write,*) 'time: ',time*params%cpp%time
n0t=(ne0-n_ne)/2._rp*(tanh((time-1.75*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,pchunk
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,pchunk
ne(cc) = ne0
end do
!$OMP END SIMD
END SELECT
SELECT CASE (TRIM(P%Te_profile))
CASE('FLAT')
!$OMP SIMD
do cc=1_idef,pchunk
Te(cc) = Te0
end do
!$OMP END SIMD
CASE('SPONG')
!$OMP SIMD
do cc=1_idef,pchunk
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('MST_FSA')
!$OMP SIMD
do cc=1_idef,pchunk
rm(cc)=sqrt((Y_R(cc)-R0a)**2+(Y_Z(cc)-Z0)**2)
r_a(cc)=rm(cc)/a
Te(cc) = (Te0-n_Te)*(1._rp-r_a(cc)**8._rp)**4._rp+n_Te
!write(6,*) 'T_e',Te(cc)*params%cpp%temperature/C_E
end do
!$OMP END SIMD
CASE DEFAULT
!$OMP SIMD
do cc=1_idef,pchunk
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,pchunk
Zeff(cc) = P%Zeffo
end do
!$OMP END SIMD
CASE('SPONG')
!$OMP SIMD
do cc=1_idef,pchunk
Zeff(cc) = P%Zeffo
end do
!$OMP END SIMD
CASE DEFAULT
!$OMP SIMD
do cc=1_idef,pchunk
Zeff(cc) = P%Zeffo
end do
!$OMP END SIMD
END SELECT
! write(output_unit_write,*) PSIpN(1)
! write(output_unit_write,'("ne: "E17.10)') ne(1)/params%cpp%length**3
! write(output_unit_write,'("rm_RE: "E17.10)') rm_RE(1)
end subroutine analytical_profiles_p