SUBROUTINE integrated_spectral_density(params,spp)
IMPLICIT NONE
TYPE(KORC_PARAMS), INTENT(IN) :: params
TYPE(SPECIES), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: spp
CHARACTER(MAX_STRING_LENGTH) :: var_name
REAL(rp), DIMENSION(3) :: binorm
REAL(rp), DIMENSION(3) :: X, V, B, E
REAL(rp), DIMENSION(:), ALLOCATABLE :: P
REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: P_lambda_p
REAL(rp), DIMENSION(:,:), ALLOCATABLE :: np_lambda_p
REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: P_lambda_t
REAL(rp), DIMENSION(:,:), ALLOCATABLE :: np_lambda_t
REAL(rp), DIMENSION(:), ALLOCATABLE :: zeta
REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: np_p
REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: np_t
REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: Psyn_lambda_p
REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: PTot_p
REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: Psyn_lambda_t
REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: PTot_t
REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: array3D
REAL(rp), DIMENSION(:,:), ALLOCATABLE :: array2D
REAL(rp), DIMENSION(:), ALLOCATABLE :: array1D
REAL(rp) :: Dtor
REAL(rp) :: phi
REAL(rp) :: Rpol, Zpol
REAL(rp) :: q, m, k, u, g, lc
REAL(rp), DIMENSION(:), ALLOCATABLE :: photon_energy
INTEGER :: ii,jj,kk,ll,ss,pp
REAL(rp), DIMENSION(:), ALLOCATABLE :: send_buffer, receive_buffer
INTEGER :: numel, mpierr
REAL(rp) :: units
ALLOCATE(Psyn_lambda_p(pplane%grid_dims(1),pplane%grid_dims(2),cam%ntor_sections,params%num_species))
ALLOCATE(PTot_p(pplane%grid_dims(1),pplane%grid_dims(2),cam%ntor_sections,params%num_species))
ALLOCATE(Psyn_lambda_t(pplane%grid_dims(1),pplane%grid_dims(2),cam%ntor_sections,params%num_species))
ALLOCATE(PTot_t(pplane%grid_dims(1),pplane%grid_dims(2),cam%ntor_sections,params%num_species))
ALLOCATE(np_p(pplane%grid_dims(1),pplane%grid_dims(2),cam%ntor_sections,params%num_species))
ALLOCATE(np_t(pplane%grid_dims(1),pplane%grid_dims(2),cam%ntor_sections,params%num_species))
ALLOCATE(P(cam%Nlambda))
ALLOCATE(P_lambda_p(cam%Nlambda,cam%ntor_sections,params%num_species))
ALLOCATE(np_lambda_p(cam%ntor_sections,params%num_species))
ALLOCATE(P_lambda_t(cam%Nlambda,cam%ntor_sections,params%num_species))
ALLOCATE(np_lambda_t(cam%ntor_sections,params%num_species))
ALLOCATE(zeta(cam%Nlambda))
ALLOCATE(photon_energy(cam%Nlambda))
np_p = 0.0_rp
np_t = 0.0_rp
Psyn_lambda_p = 0.0_rp
PTot_p = 0.0_rp
Psyn_lambda_t = 0.0_rp
PTot_t = 0.0_rp
P = 0.0_rp
P_lambda_p = 0.0_rp
P_lambda_t = 0.0_rp
np_lambda_p = 0.0_rp
np_lambda_t = 0.0_rp
photon_energy = C_h*C_C/cam%lambda
Dtor = 2.0_rp*C_PI/REAL(cam%ntor_sections,rp)
do ss=1_idef,params%num_species
q = ABS(spp(ss)%q)*params%cpp%charge
m = spp(ss)%m*params%cpp%mass
!$OMP PARALLEL DO FIRSTPRIVATE(q,m,photon_energy,Dtor) &
!$OMP& PRIVATE(binorm,X,V,B,E,k,u,g,lc,ii,jj,kk,ll,pp,zeta,P, &
!$OMP& Rpol,Zpol,phi) &
!$OMP& SHARED(params,spp,ss,Psyn_lambda_p,PTot_p,Psyn_lambda_t, &
!$OMP& PTot_t,np_p,np_t,P_lambda_p,np_lambda_p,P_lambda_t,np_lambda_t)
do pp=1_idef,spp(ss)%ppp
if ( spp(ss)%vars%flag(pp) .EQ. 1_idef ) then
V = spp(ss)%vars%V(pp,:)*params%cpp%velocity
X = spp(ss)%vars%X(pp,:)*params%cpp%length
g = spp(ss)%vars%g(pp)
B = spp(ss)%vars%B(pp,:)*params%cpp%Bo
E = spp(ss)%vars%E(pp,:)*params%cpp%Eo
binorm = cross(V,E) + cross(V,cross(V,B))
u = SQRT(DOT_PRODUCT(V,V))
k = q*SQRT(DOT_PRODUCT(binorm,binorm))/(spp(ss)%vars%g(pp)*m*u**3)
lc = (4.0_rp*C_PI/3.0_rp)/(k*g**3)
zeta = lc/cam%lambda
Rpol = SQRT(SUM(X(1:2)**2))
Zpol = X(3)
ii = FLOOR((Rpol - pplane%Rmin)/pplane%DR) + 1_idef
jj = FLOOR((Zpol + ABS(pplane%Zmin))/pplane%DZ) + 1_idef
phi = ATAN2(X(2),X(1))
if (phi.LT.0.0_rp) phi = 2.0_rp*C_PI + phi
kk = floor(phi/Dtor) + 1_idef
call P_integral(zeta,P)
P = (C_C*C_E**2)*P/(SQRT(3.0_rp)*C_E0*g**2*cam%lambda**3)
if (spp(ss)%vars%eta(pp) .LT. 90.0_rp) then
P_lambda_p(:,kk,ss) = P_lambda_p(:,kk,ss) + P
np_lambda_p(kk,ss) = np_lambda_p(kk,ss) + 1.0_rp
Psyn_lambda_p(ii,jj,kk,ss) = Psyn_lambda_p(ii,jj,kk,ss) + trapz(cam%lambda,P)
PTot_p(ii,jj,kk,ss) = PTot_p(ii,jj,kk,ss) + spp(ss)%vars%Prad(pp);
np_p(ii,jj,kk,ss) = np_p(ii,jj,kk,ss) + 1_idef
else
P_lambda_t(:,kk,ss) = P_lambda_t(:,kk,ss) + P
np_lambda_t(kk,ss) = np_lambda_t(kk,ss) + 1.0_rp
Psyn_lambda_t(ii,jj,kk,ss) = Psyn_lambda_t(ii,jj,kk,ss) + trapz(cam%lambda,P)
PTot_t(ii,jj,kk,ss) = PTot_t(ii,jj,kk,ss) + spp(ss)%vars%Prad(pp);
np_t(ii,jj,kk,ss) = np_t(ii,jj,kk,ss) + 1_idef
end if
end if ! if confined
end do ! particles
!$OMP END PARALLEL DO
end do ! species
! * * * * * * * IMPORTANT * * * * * * *
! * * * * * * * * * * * * * * * * * * *
! Here Psyn_lambda has units of Watts/m
! or (photons/s)(m^-1). See above.
! * * * * * * * * * * * * * * * * * * *
! * * * * * * * IMPORTANT * * * * * * *
if (params%mpi_params%rank.EQ.0_idef) then
if (.NOT.cam%toroidal_sections) then
ALLOCATE(array3D(pplane%grid_dims(1),pplane%grid_dims(2),params%num_species))
ALLOCATE(array2D(cam%Nlambda,params%num_species))
ALLOCATE(array1D(params%num_species))
end if
end if
if (params%mpi_params%nmpi.GT.1_idef) then
numel = pplane%grid_dims(1)*pplane%grid_dims(2)*cam%ntor_sections*params%num_species
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(Psyn_lambda_p,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
Psyn_lambda_p = &
RESHAPE(receive_buffer,(/cam%num_pixels(1),cam%num_pixels(2),cam%ntor_sections,params%num_species/))
var_name = 'Psyn_p_pplane'
if (cam%toroidal_sections) then
call save_snapshot_var(params,Psyn_lambda_p,var_name)
else
array3D = SUM(Psyn_lambda_p,3)
call save_snapshot_var(params,array3D,var_name)
end if
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(Psyn_lambda_t,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
Psyn_lambda_t = &
RESHAPE(receive_buffer,(/cam%num_pixels(1),cam%num_pixels(2),cam%ntor_sections,params%num_species/))
var_name = 'Psyn_t_pplane'
if (cam%toroidal_sections) then
call save_snapshot_var(params,Psyn_lambda_t,var_name)
else
array3D = SUM(Psyn_lambda_t,3)
call save_snapshot_var(params,array3D,var_name)
end if
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(np_p,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
np_p = RESHAPE(receive_buffer,(/cam%num_pixels(1),cam%num_pixels(2),cam%ntor_sections,params%num_species/))
var_name = 'np_p_pplane'
if (cam%toroidal_sections) then
call save_snapshot_var(params,np_p,var_name)
else
array3D = SUM(np_p,3)
call save_snapshot_var(params,array3D,var_name)
end if
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(np_t,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
np_t = RESHAPE(receive_buffer,(/cam%num_pixels(1),cam%num_pixels(2),cam%ntor_sections,params%num_species/))
var_name = 'np_t_pplane'
if (cam%toroidal_sections) then
call save_snapshot_var(params,np_t,var_name)
else
array3D = SUM(np_t,3)
call save_snapshot_var(params,array3D,var_name)
end if
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(PTot_p,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
PTot_p = RESHAPE(receive_buffer,(/cam%num_pixels(1),cam%num_pixels(2),cam%ntor_sections,params%num_species/))
units = params%cpp%mass*(params%cpp%velocity**3)/params%cpp%length
PTot_p = units*PTot_p ! (Watts)
var_name = 'PTot_p_pplane'
if (cam%toroidal_sections) then
call save_snapshot_var(params,PTot_p,var_name)
else
array3D = SUM(PTot_p,3)
call save_snapshot_var(params,array3D,var_name)
end if
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(PTot_t,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
PTot_t = RESHAPE(receive_buffer,(/cam%num_pixels(1),cam%num_pixels(2),cam%ntor_sections,params%num_species/))
units = params%cpp%mass*(params%cpp%velocity**3)/params%cpp%length
PTot_t = units*PTot_t ! (Watts)
var_name = 'PTot_t_pplane'
if (cam%toroidal_sections) then
call save_snapshot_var(params,PTot_t,var_name)
else
array3D = SUM(PTot_t,3)
call save_snapshot_var(params,array3D,var_name)
end if
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
numel = cam%ntor_sections*params%num_species
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(np_lambda_p,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
np_lambda_p = RESHAPE(receive_buffer,(/cam%ntor_sections,params%num_species/))
var_name = 'np_p_lambda'
if (cam%toroidal_sections) then
call save_snapshot_var(params,np_lambda_p,var_name)
else
array1D = SUM(np_lambda_p,1)
call save_snapshot_var(params,array1D,var_name)
end if
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(np_lambda_t,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
np_lambda_t = RESHAPE(receive_buffer,(/cam%ntor_sections,params%num_species/))
var_name = 'np_t_lambda'
if (cam%toroidal_sections) then
call save_snapshot_var(params,np_lambda_t,var_name)
else
array1D = SUM(np_lambda_t,1)
call save_snapshot_var(params,array1D,var_name)
end if
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
numel = cam%Nlambda*cam%ntor_sections*params%num_species
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(P_lambda_p,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
P_lambda_p = RESHAPE(receive_buffer,(/cam%Nlambda,cam%ntor_sections,params%num_species/))
var_name = 'P_p_lambda'
if (cam%toroidal_sections) then
call save_snapshot_var(params,P_lambda_p,var_name)
else
array2D = SUM(P_lambda_p,2)
call save_snapshot_var(params,array2D,var_name)
end if
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(P_lambda_t,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
P_lambda_t = RESHAPE(receive_buffer,(/cam%Nlambda,cam%ntor_sections,params%num_species/))
var_name = 'P_t_lambda'
if (cam%toroidal_sections) then
call save_snapshot_var(params,P_lambda_t,var_name)
else
array2D = SUM(P_lambda_t,2)
call save_snapshot_var(params,array2D,var_name)
end if
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
CALL MPI_BARRIER(MPI_COMM_WORLD,mpierr)
else
var_name = 'np_p_pplane'
if (cam%toroidal_sections) then
call save_snapshot_var(params,np_p,var_name)
else
array3D = SUM(np_p,3)
call save_snapshot_var(params,array3D,var_name)
end if
var_name = 'np_t_pplane'
if (cam%toroidal_sections) then
call save_snapshot_var(params,np_t,var_name)
else
array3D = SUM(np_t,3)
call save_snapshot_var(params,array3D,var_name)
end if
var_name = 'Psyn_p_pplane'
if (cam%toroidal_sections) then
call save_snapshot_var(params,Psyn_lambda_p,var_name)
else
array3D = SUM(Psyn_lambda_p,3)
call save_snapshot_var(params,array3D,var_name)
end if
var_name = 'Psyn_t_pplane'
if (cam%toroidal_sections) then
call save_snapshot_var(params,Psyn_lambda_t,var_name)
else
array3D = SUM(Psyn_lambda_t,3)
call save_snapshot_var(params,array3D,var_name)
end if
var_name = 'P_p_lambda'
if (cam%toroidal_sections) then
call save_snapshot_var(params,P_lambda_p,var_name)
else
array2D = SUM(P_lambda_p,2)
call save_snapshot_var(params,array2D,var_name)
end if
var_name = 'P_t_lambda'
if (cam%toroidal_sections) then
call save_snapshot_var(params,P_lambda_t,var_name)
else
array2D = SUM(P_lambda_t,2)
call save_snapshot_var(params,array2D,var_name)
end if
var_name = 'np_p_lambda'
if (cam%toroidal_sections) then
call save_snapshot_var(params,np_lambda_p,var_name)
else
array1D = SUM(np_lambda_p,1)
call save_snapshot_var(params,array1D,var_name)
end if
var_name = 'np_t_lambda'
if (cam%toroidal_sections) then
call save_snapshot_var(params,np_lambda_t,var_name)
else
array1D = SUM(np_lambda_t,1)
call save_snapshot_var(params,array1D,var_name)
end if
units = params%cpp%mass*(params%cpp%velocity**3)/params%cpp%length
PTot_p = units*PTot_p ! (Watts)
var_name = 'PTot_p_pplane'
if (cam%toroidal_sections) then
call save_snapshot_var(params,PTot_p,var_name)
else
array3D = SUM(PTot_p,3)
call save_snapshot_var(params,array3D,var_name)
end if
units = params%cpp%mass*(params%cpp%velocity**3)/params%cpp%length
PTot_t = units*PTot_t ! (Watts)
var_name = 'PTot_t_pplane'
if (cam%toroidal_sections) then
call save_snapshot_var(params,PTot_t,var_name)
else
array3D = SUM(PTot_t,3)
call save_snapshot_var(params,array3D,var_name)
end if
end if
DEALLOCATE(np_p)
DEALLOCATE(np_t)
DEALLOCATE(Psyn_lambda_p)
DEALLOCATE(PTot_p)
DEALLOCATE(Psyn_lambda_t)
DEALLOCATE(PTot_t)
DEALLOCATE(P)
DEALLOCATE(P_lambda_p)
DEALLOCATE(np_lambda_p)
DEALLOCATE(P_lambda_t)
DEALLOCATE(np_lambda_t)
DEALLOCATE(zeta)
DEALLOCATE(photon_energy)
if (ALLOCATED(array3D)) DEALLOCATE(array3D)
if (ALLOCATED(array2D)) DEALLOCATE(array2D)
if (ALLOCATED(array1D)) DEALLOCATE(array1D)
END SUBROUTINE integrated_spectral_density