diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 8beff240b7b7194f9d0efd8bcb8780c4715c1b2c..50d383b3ada49c5ca740ea87839482e16f7da3a4 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -210,34 +210,34 @@ end
%% Build grids
-Nkx = numel(kx);
-if Nkx > 1
- dkx = kx(2);
- Lx = 2*pi/dkx;
+Nky = numel(ky);
+if Nky > 1
+ dky = ky(2);
+ Ly = 2*pi/dky;
else
- dkx = 0;
- Lx = 0;
+ dky = 0;
+ Ly = 0;
end
-[~,ikx0] = min(abs(kx));
-Nx = 2*Nkx-1;
-x = linspace(-Lx/2,Lx/2,Nx+1); x = x(1:end-1);
+[~,iky0] = min(abs(ky));
+Ny = 2*Nky-1;
+y = linspace(-Ly/2,Ly/2,Ny+1); y = y(1:end-1);
-Nky = numel(ky);
-if Nky > 1
- dky = ky(2);
- Ly = 2*pi/dky;
+Nkx = numel(kx);
+if Nkx > 1
+ dkx = kx(2);
+ Lx = 2*pi/dkx;
else
- dky = 0;
- Ly = 0;
+ dkx = 0;
+ Lx = 0;
end
-[~,iky0] = min(abs(ky));
-Ny = Nky;
-y = linspace(-Ly/2,Ly/2,Ny+1); y = y(1:end-1);
+[~,ikx0] = min(abs(kx));
+Nx = Nkx;
+x = linspace(-Lx/2,Lx/2,Nx+1); x = x(1:end-1);
Nz = numel(z);
-[KY,KX] = meshgrid(ky,kx);
-KPERP2 = KY.^2+KX.^2;
+[KX,KY] = meshgrid(kx,ky);
+KPERP2 = KX.^2+KY.^2;
%% Add everything in output structure
% scaling
DATA.scale = (1/Nx/Ny)^2;
diff --git a/matlab/compute/compute_3D_zpinch_growth_rate.m b/matlab/compute/compute_3D_zpinch_growth_rate.m
index 1cfeed825d85e0dc2d7235bc3e2bc7b2422afce0..3a053fa6ae60ad57daaac0652adfc894d9844b4c 100644
--- a/matlab/compute/compute_3D_zpinch_growth_rate.m
+++ b/matlab/compute/compute_3D_zpinch_growth_rate.m
@@ -7,20 +7,20 @@ t = DATA.Ts3D;
nkx = DATA.Nkx; nky = DATA.Nky; nkz = DATA.Nz;
% Remove AA part
if DATA.Nx > 1
- ikxnz = abs(DATA.PHI(:,1,1,1)) > 0;
+ ikxnz = abs(DATA.PHI(1,:,1,1)) > 0;
else
- ikxnz = abs(DATA.PHI(:,1,1,1)) > -1;
+ ikxnz = abs(DATA.PHI(1,:,1,1)) > -1;
end
-ikynz = abs(DATA.PHI(1,:,1,1)) > 0;
+ikynz = abs(DATA.PHI(:,1,1,1)) > 0;
-phi = fft(DATA.PHI(ikxnz,ikynz,:,:),[],3);
+phi = fft(DATA.PHI(ikynz,ikxnz,:,:),[],3);
-omega = zeros(nkx,nky,nkz);
+omega = zeros(nky,nkx,nkz);
for iz = 1:nkz
for iy = 1:nky
for ix = 1:nkx
- omega(ix,iy,iz) = LinearFit_s(t(its:ite),squeeze(abs(phi(ix,iy,iz,its:ite))));
+ omega(iy,ix,iz) = LinearFit_s(t(its:ite),squeeze(abs(phi(iy,ix,iz,its:ite))));
end
end
end
@@ -34,14 +34,14 @@ poskz = kz>=0;
kxeq0 = kx==0;
kzeq0 = kz==0;
-omega = omega(poskx,posky,poskz);
+omega = omega(posky,poskx,poskz);
FIGURE.fig = figure;
nplots = OPTIONS.kxky + OPTIONS.kzkx + OPTIONS.kzky;
iplot = 1;
if OPTIONS.kxky
-[Y_XY,X_XY] = meshgrid(ky(posky),kx(poskx));
+[X_XY,Y_XY] = meshgrid(ky(posky),kx(poskx));
subplot(1,nplots,iplot)
if ~OPTIONS.keq0
toplot = squeeze(max(real(omega(:,:,:)),[],3));
diff --git a/matlab/compute/compute_fluxtube_growth_rate.m b/matlab/compute/compute_fluxtube_growth_rate.m
index b92ecbe9d8da2065017cfd5abbfc911ddbc08fce..fad8aecfae4bb2d75d83d0b704de79ed0c620f67 100644
--- a/matlab/compute/compute_fluxtube_growth_rate.m
+++ b/matlab/compute/compute_fluxtube_growth_rate.m
@@ -2,13 +2,13 @@ function [ linear_gr ] = compute_fluxtube_growth_rate(DATA, TRANGE, PLOT)
% Remove AA part
if DATA.Nx > 1
- ikxnz = abs(DATA.PHI(:,1,1,1)) > 0;
+ ikxnz = abs(DATA.PHI(1,:,1,1)) > 0;
else
- ikxnz = abs(DATA.PHI(:,1,1,1)) > -1;
+ ikxnz = abs(DATA.PHI(1,:,1,1)) > -1;
end
-ikynz = abs(DATA.PHI(1,:,1,1)) > 0;
+ikynz = abs(DATA.PHI(:,1,1,1)) > 0;
-phi = DATA.PHI(ikxnz,ikynz,:,:);
+phi = DATA.PHI(ikynz,ikxnz,:,:);
t = DATA.Ts3D;
[~,its] = min(abs(t-TRANGE(1)));
@@ -22,13 +22,13 @@ for it = its+1:ite
phi_n = phi(:,:,:,it);
phi_nm1 = phi(:,:,:,it-1);
dt = t(it)-t(it-1);
- ZS = sum(sum(phi_nm1,1),3);
+ ZS = sum(sum(phi_nm1,2),3);
wl = log(phi_n./phi_nm1)/dt;
- w_ky(:,is) = sum(sum(wl.*phi_nm1,1),3)./ZS;
+ w_ky(:,is) = squeeze(sum(sum(wl.*phi_nm1,2),3)./ZS);
for iky = 1:numel(w_ky(:,is))
- ce(iky,is) = abs(sum(sum(abs(w_ky(iky,is)-wl(:,iky,:)).^2.*phi_nm1(:,iky,:),1),3)./ZS(:,iky,:));
+ ce(iky,is) = abs(sum(sum(abs(w_ky(iky,is)-wl(iky,:,:)).^2.*phi_nm1(iky,:,:),2),3)./ZS(iky,:,:));
end
is = is + 1;
end
diff --git a/matlab/compute/ifourier_GENE.m b/matlab/compute/ifourier_GENE.m
index f5d98f75ede81706f721ddcb2204ee46c91c2dc9..3c3a454a2032a384d73508347e8301d7f618454e 100644
--- a/matlab/compute/ifourier_GENE.m
+++ b/matlab/compute/ifourier_GENE.m
@@ -5,43 +5,43 @@ function [ field_r ] = ifourier_GENE( field_c )
% comparison purpose.
%% Original
-% [nx,nky,nz]=size(field_c);
-% %extend to whole ky by imposing reality condition
-% ny=2*nky-1;
-%
-% if ny~=1
-% %note, we need one extra point which we set to zero for the ifft
-% spectrumKxKyZ=zeros(nx,ny,nz);
-% spectrumKxKyZ(:,1:nky,:)=field_c(:,:,:);
-% spectrumKxKyZ(1,(nky+1):(ny),:)=conj(field_c(1,nky:-1:2,:));
-% spectrumKxKyZ(2:nx,(nky+1):(ny),:)=conj(field_c(nx:-1:2,nky:-1:2,:));
-% else
-% %pad with zeros to interpolate on fine scale
-% ny=20;
-% spectrumKxKyZ=zeros(nx,ny,nz);
-% spectrumKxKyZ(:,2,:)=field_c(:,:,:);
-% end
-%
-% %inverse fft, symmetric as we are using real data
-% spectrumXKyZ=nx*ifft(spectrumKxKyZ,[],1);
-% field_r=ny*ifft(spectrumXKyZ,[],2,'symmetric');
-% clear spectrumKxKyZ
-
-%% Adapted for HeLaZ representation
-[nkx,ny,nz]=size(field_c);
+[nky,nx,nz]=size(field_c);
%extend to whole ky by imposing reality condition
-nx=2*nkx-1;
+ny=2*nky-1;
-%note, we need one extra point which we set to zero for the ifft
-spectrumKxKyZ=zeros(nx,ny,nz);
-spectrumKxKyZ(1:nkx,:,:)=field_c(:,:,:);
-spectrumKxKyZ((nkx+1):(nx),1,:)=conj(field_c(nkx:-1:2,1,:));
-spectrumKxKyZ((nkx+1):(nx),2:ny,:)=conj(field_c(nkx:-1:2,ny:-1:2,:));
+if ny~=1
+ %note, we need one extra point which we set to zero for the ifft
+ spectrumKxKyZ=zeros(ny,nx,nz);
+ spectrumKxKyZ(1:nky,:,:)=field_c(:,:,:);
+ spectrumKxKyZ((nky+1):(ny),1,:)=conj(field_c(nky:-1:2,1,:));
+ spectrumKxKyZ((nky+1):(ny),2:nx,:)=conj(field_c(nky:-1:2,nx:-1:2,:));
+else
+ %pad with zeros to interpolate on fine scale
+ ny=20;
+ spectrumKxKyZ=zeros(nx,ny,nz);
+ spectrumKxKyZ(:,2,:)=field_c(:,:,:);
+end
%inverse fft, symmetric as we are using real data
spectrumXKyZ=nx*ifft(spectrumKxKyZ,[],1);
field_r=ny*ifft(spectrumXKyZ,[],2,'symmetric');
clear spectrumKxKyZ
+
+%% Adapted for HeLaZ old representation
+% [nkx,ny,nz]=size(field_c);
+% %extend to whole ky by imposing reality condition
+% nx=2*nkx-1;
+%
+% %note, we need one extra point which we set to zero for the ifft
+% spectrumKxKyZ=zeros(nx,ny,nz);
+% spectrumKxKyZ(1:nkx,:,:)=field_c(:,:,:);
+% spectrumKxKyZ((nkx+1):(nx),1,:)=conj(field_c(nkx:-1:2,1,:));
+% spectrumKxKyZ((nkx+1):(nx),2:ny,:)=conj(field_c(nkx:-1:2,ny:-1:2,:));
+%
+% %inverse fft, symmetric as we are using real data
+% spectrumXKyZ=nx*ifft(spectrumKxKyZ,[],1);
+% field_r=ny*ifft(spectrumXKyZ,[],2,'symmetric');
+% clear spectrumKxKyZ
end
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index 5dc0088697a234b24d65ec6602cd5d99aeed74e2..d0c1a9ce80f8ee8571994e42775f32d3fb397191 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -4,7 +4,7 @@ NORMALIZED = OPTIONS.NORMALIZED;
Nma = OPTIONS.NMA; %Number moving average
t = OPTIONS.TIME;
iz = OPTIONS.iz;
-[~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(:,1,1,:))),2)));
+[~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(1,:,1,:))),2)));
FRAMES = zeros(size(OPTIONS.TIME));
@@ -19,18 +19,18 @@ MODES_SELECTOR = i; %(1:Zonal, 2: NZonal, 3: ky=kx)
if MODES_SELECTOR == 1
if NORMALIZED
- plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES)))./max(abs(squeeze(x(ik,1,iz,FRAMES)))),Nma);
+ plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES)))./max(abs(squeeze(x(1,ik,iz,FRAMES)))),Nma);
else
- plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES))),Nma);
+ plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES))),Nma);
end
kstr = 'k_x';
k = DATA.kx;
MODESTR = 'Zonal modes';
elseif MODES_SELECTOR == 2
if NORMALIZED
- plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES)))./max(abs(squeeze(x(1,ik,iz,FRAMES)))),Nma);
+ plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES)))./max(abs(squeeze(x(ik,1,iz,FRAMES)))),Nma);
else
- plt = @(x,ik) movmean(abs(squeeze(x(1,ik,iz,FRAMES))),Nma);
+ plt = @(x,ik) movmean(abs(squeeze(x(ik,1,iz,FRAMES))),Nma);
end
kstr = 'k_y';
k = DATA.ky;
diff --git a/matlab/load/load_3D_data.m b/matlab/load/load_3D_data.m
index 2921908adebde09d65ccbde6c0124a7e92ca416b..0a4fee17d6e512396b661d009476b1d94a54a42f 100644
--- a/matlab/load/load_3D_data.m
+++ b/matlab/load/load_3D_data.m
@@ -7,7 +7,7 @@ function [ data, time, dt ] = load_3D_data( filename, variablename )
dt = h5readatt(filename,'/data/input','dt');
cstart= h5readatt(filename,'/data/input','start_iframe3d');
- data = zeros(numel(kx),numel(ky),numel(z),numel(time));
+ data = zeros(numel(ky),numel(kx),numel(z),numel(time));
for it = 1:numel(time)
tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
diff --git a/matlab/load/load_5D_data.m b/matlab/load/load_5D_data.m
index 2de05533001857e8dd4d28d5b8709a66f715d701..bc96f113bba9decc602b02db6c46ed0cafbcac34 100644
--- a/matlab/load/load_5D_data.m
+++ b/matlab/load/load_5D_data.m
@@ -15,7 +15,7 @@ function [ data, time, dt ] = load_5D_data( filename, variablename )
dt = h5readatt(filename,'/data/input','dt');
cstart= h5readatt(filename,'/data/input','start_iframe5d');
- data = zeros(numel(p),numel(j),numel(kx),numel(ky),numel(z),numel(time));
+ data = zeros(numel(p),numel(j),numel(ky),numel(kx),numel(z),numel(time));
for it = 1:numel(time)
tmp = h5read(filename,['/data/var5d/', variablename,'/', num2str(cstart+it,'%06d')]);
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index be4680f652dbfb54368a6d1a0419ba53c7c0c120..1957f424a2f33896164166879060e516f13ef60c 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -10,13 +10,13 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, TAVG_0, TAVG_1,stf
Qx_infty_std = std (DATA.HFLUX_X(its0D:ite0D))*SCALE;
[~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(:,1,1,:))),2)));
Ns3D = numel(DATA.Ts3D);
- [KY, KX] = meshgrid(DATA.ky, DATA.kx);
+ [KX, KY] = meshgrid(DATA.kx, DATA.ky);
z = DATA.z;
%% computations
% Compute Gamma from ifft matlab
- Gx = zeros(DATA.Nx,DATA.Ny,numel(DATA.Ts3D));
+ Gx = zeros(DATA.Ny,DATA.Nx,numel(DATA.Ts3D));
for it = 1:numel(DATA.Ts3D)
for iz = 1:DATA.Nz
Gx(:,:,it) = Gx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
@@ -26,7 +26,7 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, TAVG_0, TAVG_1,stf
end
Gx_t_mtlb = squeeze(mean(mean(Gx,1),2));
% Compute Heat flux from ifft matlab
- Qx = zeros(DATA.Nx,DATA.Ny,numel(DATA.Ts3D));
+ Qx = zeros(DATA.Ny,DATA.Nx,numel(DATA.Ts3D));
for it = 1:numel(DATA.Ts3D)
for iz = 1:DATA.Nz
Qx(:,:,it) = Qx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
@@ -71,7 +71,7 @@ mvm = @(x) movmean(x,Nmvm);
% computation
Ns3D = numel(DATA.Ts3D);
[KY, KX] = meshgrid(DATA.ky, DATA.kx);
- plt = @(x) mean(x(:,:,1,:),2);
+ plt = @(x) mean(x(:,:,1,:),1);
kycut = max(DATA.ky);
kxcut = max(DATA.kx);
LP = (abs(KY)<kycut).*(abs(KX)<kxcut); %Low pass filter
@@ -87,7 +87,7 @@ mvm = @(x) movmean(x,Nmvm);
end
switch stfname
case 'phi'
- phi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
+ phi = zeros(DATA.Ny,DATA.Nx,1,Ns3D);
for it = 1:numel(DATA.Ts3D)
phi(:,:,1,it) = ifourier_GENE(compz(DATA.PHI(:,:,:,it)));
end
diff --git a/matlab/process_field.m b/matlab/process_field.m
index 5762fc84b51c6057ac03c95b353644ae0c11827a..5258ef7cdc949185f2d72f9d52af005053424220 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -17,28 +17,28 @@ LTXNAME = OPTIONS.NAME;
switch OPTIONS.PLAN
case 'xy'
XNAME = '$x$'; YNAME = '$y$';
- [Y,X] = meshgrid(DATA.y,DATA.x);
+ [X,Y] = meshgrid(DATA.x,DATA.y);
REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
case 'xz'
XNAME = '$x$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.x);
- REALP = 1; COMPDIM = 2; SCALE = 0;
+ REALP = 1; COMPDIM = 1; SCALE = 0;
case 'yz'
XNAME = '$y$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.y);
- REALP = 1; COMPDIM = 1; SCALE = 0;
+ REALP = 1; COMPDIM = 2; SCALE = 0;
case 'kxky'
XNAME = '$k_x$'; YNAME = '$k_y$';
- [Y,X] = meshgrid(DATA.ky,DATA.kx);
+ [X,Y] = meshgrid(DATA.kx,DATA.ky);
REALP = 0; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
case 'kxz'
XNAME = '$k_x$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.kx);
- REALP = 0; COMPDIM = 2; POLARPLOT = 0; SCALE = 0;
+ REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0;
case 'kyz'
XNAME = '$k_y$'; YNAME = '$z$';
[Y,X] = meshgrid(DATA.z,DATA.ky);
- REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0;
+ REALP = 0; COMPDIM = 2; POLARPLOT = 0; SCALE = 0;
case 'sx'
XNAME = '$v_\parallel$'; YNAME = '$\mu$';
[Y,X] = meshgrid(OPTIONS.XPERP,OPTIONS.SPAR);
@@ -110,7 +110,7 @@ end
switch REALP
case 1 % Real space plot
INTERP = OPTIONS.INTERP;
- process = @(x) real(fftshift(ifft2(x,Nx,Ny)));
+ process = @(x) real(fftshift(ifft2(x,Ny,Nx)));
shift_x = @(x) x;
shift_y = @(x) x;
case 0 % Frequencies plot
diff --git a/src/advance_field.F90 b/src/advance_field.F90
index a26bda44ebef6ee7a4ed842ba7646536e0f9fa3c..ae48815ed2887bc62d7ad082fe217dab25f008a4 100644
--- a/src/advance_field.F90
+++ b/src/advance_field.F90
@@ -32,7 +32,7 @@ CONTAINS
p_int = parray_e(ip)
DO ij=ijs_e,ije_e
IF((CLOS .NE. 1) .OR. (ip-1+2*(ij-1)+1 .LE. dmaxe))&
- CALL advance_field(moments_e(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,:), moments_rhs_e(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,:))
+ CALL advance_field(moments_e(ip,ij,ikys:ikye,ikxs:ikxe,izs:ize,:), moments_rhs_e(ip,ij,ikys:ikye,ikxs:ikxe,izs:ize,:))
ENDDO
ENDDO
ENDIF
@@ -41,7 +41,7 @@ CONTAINS
DO ij=ijs_i,ije_i
j_int = jarray_i(ij)
IF((CLOS .NE. 1) .OR. (ip-1+2*(ij-1)+1 .LE. dmaxi))&
- CALL advance_field(moments_i(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,:), moments_rhs_i(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,:))
+ CALL advance_field(moments_i(ip,ij,ikys:ikye,ikxs:ikxe,izs:ize,:), moments_rhs_i(ip,ij,ikys:ikye,ikxs:ikxe,izs:ize,:))
ENDDO
ENDDO
! Execution time end
@@ -60,21 +60,21 @@ CONTAINS
use initial_par, ONLY: ACT_ON_MODES
IMPLICIT NONE
- COMPLEX(dp), DIMENSION ( ikxs:ikxe, ikys:ikye, izs:ize, ntimelevel ) :: f
- COMPLEX(dp), DIMENSION ( ikxs:ikxe, ikys:ikye, izs:ize, ntimelevel ) :: f_rhs
+ COMPLEX(dp), DIMENSION ( ikys:ikye, ikxs:ikxe, izs:ize, ntimelevel ) :: f
+ COMPLEX(dp), DIMENSION ( ikys:ikye, ikxs:ikxe, izs:ize, ntimelevel ) :: f_rhs
INTEGER :: istage
SELECT CASE (updatetlevel)
CASE(1)
DO istage=1,ntimelevel
- f(ikxs:ikxe,ikys:ikye,izs:ize,1) = f(ikxs:ikxe,ikys:ikye,izs:ize,1) &
- + dt*b_E(istage)*f_rhs(ikxs:ikxe,ikys:ikye,izs:ize,istage)
+ f(ikys:ikye,ikxs:ikxe,izs:ize,1) = f(ikys:ikye,ikxs:ikxe,izs:ize,1) &
+ + dt*b_E(istage)*f_rhs(ikys:ikye,ikxs:ikxe,izs:ize,istage)
END DO
CASE DEFAULT
- f(ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel) = f(ikxs:ikxe,ikys:ikye,izs:ize,1);
+ f(ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel) = f(ikys:ikye,ikxs:ikxe,izs:ize,1);
DO istage=1,updatetlevel-1
- f(ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel) = f(ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel) + &
- dt*A_E(updatetlevel,istage)*f_rhs(ikxs:ikxe,ikys:ikye,izs:ize,istage)
+ f(ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel) = f(ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel) + &
+ dt*A_E(updatetlevel,istage)*f_rhs(ikys:ikye,ikxs:ikxe,izs:ize,istage)
END DO
END SELECT
END SUBROUTINE advance_field
diff --git a/src/advance_field_adapt.F90 b/src/advance_field_adapt.F90
deleted file mode 100644
index aed2a7e3d243058e6c05cb152617437cf7e59f80..0000000000000000000000000000000000000000
--- a/src/advance_field_adapt.F90
+++ /dev/null
@@ -1,80 +0,0 @@
-MODULE advance_field_routine
-
-USE prec_const
-implicit none
-
-CONTAINS
-
- SUBROUTINE advance_time_level
-
- USE basic
- USE time_integration
- use prec_const
- IMPLICIT NONE
-
- CALL set_updatetlevel(mod(updatetlevel,ntimelevel)+1)
-
- END SUBROUTINE advance_time_level
-
-
-
- SUBROUTINE advance_field( f, f_rhs )
-
- USE basic
- USE time_integration
- USE array
- USE grid
- use prec_const
- IMPLICIT NONE
-
- COMPLEX(dp), DIMENSION ( ikxs:ikxe, ikys:ikye, ntimelevel ) :: f
- COMPLEX(dp), DIMENSION ( ikxs:ikxe, ikys:ikye, ntimelevel ) :: f_rhs
- REAL(dp) :: error
- INTEGER :: istage
-
- SELECT CASE (updatetlevel)
-
- CASE(1)
- SELECT CASE (numerical_scheme)
-
- CASE ('DOPRI5_ADAPT')
- error = 0._dp
- DO iky=ikys,ikye
- DO ikx=ikxs,ikxe
- fs = f(ikx,iky,1)
- DO istage=1,ntimelevel
- f(ikx,iky,1) = f(ikx,iky,1) + dt*b_E(istage)*f_rhs(ikx,iky,istage)
- fs = fs + dt*b_Es(istage)*f_rhs(ikx,iky,istage)
- END DO
- IF ( ABS(f(ikx,iky,1) - fs) .GT. error ) THEN
- error = ABS(f(ikx,iky,1) - fs)
- ENDIF
- END DO
- END DO
- IF (error > TOL)
- CASE DEFAULT
- DO iky=ikys,ikye
- DO ikx=ikxs,ikxe
- fs = f(ikx,iky,1)
- DO istage=1,ntimelevel
- f(ikx,iky,1) = f(ikx,iky,1) + dt*b_E(istage)*f_rhs(ikx,iky,istage)
- END DO
- END DO
- END DO
- END SELECT
-
- CASE DEFAULT
- DO iky=ikys,ikye
- DO ikx=ikxs,ikxe
- f(ikx,iky,updatetlevel) = f(ikx,iky,1);
- DO istage=1,updatetlevel-1
- f(ikx,iky,updatetlevel) = f(ikx,iky,updatetlevel) + &
- dt*A_E(updatetlevel,istage)*f_rhs(ikx,iky,istage)
- END DO
- END DO
- END DO
- END SELECT
-
- END SUBROUTINE advance_field
-
-END MODULE advance_field_routine
diff --git a/src/auxval.F90 b/src/auxval.F90
index cebb78c0d262f8ab02dc75e508cb8a1125623691..7d6fd969feaa1e1c98888118650c5604f9bc583e 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -15,9 +15,11 @@ subroutine auxval
IF (my_id .EQ. 0) WRITE(*,*) '=== Set auxiliary values ==='
IF (LINEARITY .NE. 'linear') THEN
+ write(*,*) 'FFTW3 y-grid distribution'
CALL init_grid_distr_and_plans(Nx,Ny)
ELSE
CALL init_1Dgrid_distr
+ write(*,*) 'Manual y-grid distribution'
ENDIF
! Init the grids
CALL set_pgrid ! parallel kin (MPI distributed)
@@ -52,10 +54,10 @@ subroutine auxval
IF (my_id .EQ. 0) WRITE(*,*) ''
IF (my_id .EQ. 0) WRITE(*,*) '--------- Parallel environement ----------'
IF (my_id .EQ. 0) WRITE(*,'(A12,I3)') 'n_procs ', num_procs
- IF (my_id .EQ. 0) WRITE(*,'(A12,I3,A14,I3,A14,I3)') 'num_procs_p = ', num_procs_p, ', num_procs_kx = ', num_procs_kx, ', num_procs_z = ', num_procs_z
+ IF (my_id .EQ. 0) WRITE(*,'(A12,I3,A14,I3,A14,I3)') 'num_procs_p = ', num_procs_p, ', num_procs_ky = ', num_procs_ky, ', num_procs_z = ', num_procs_z
IF (my_id .EQ. 0) WRITE(*,*) ''
WRITE(*,'(A9,I3,A10,I3,A10,I3,A9,I3)')&
- 'my_id = ', my_id, ', rank_p = ', rank_p, ', rank_kx = ', rank_kx,', rank_z = ', rank_z
+ 'my_id = ', my_id, ', rank_p = ', rank_p, ', rank_ky = ', rank_ky,', rank_z = ', rank_z
WRITE(*,'(A22,I3,A11,I3)')&
' ips_e = ', ips_e, ', ipe_e = ', ipe_e
WRITE(*,'(A22,I3,A11,I3)')&
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index 8bbb5253e7f06d3bd3be668f9164f14e3ce3c8c3..ab7869f8b8e94743acf5b8d46ec0122f7c2acdca 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -14,8 +14,8 @@ MODULE basic
INTEGER :: comm0 ! Default communicator with a topology
INTEGER :: rank_0 ! Ranks in comm0
! Communicators for 1-dim cartesian subgrids of comm0
- INTEGER :: comm_p, comm_kx, comm_z
- INTEGER :: rank_p, rank_kx, rank_z! Ranks
+ INTEGER :: comm_p, comm_ky, comm_z
+ INTEGER :: rank_p, rank_ky, rank_z! Ranks
INTEGER :: commr_p0 ! Communicators along kx for only rank 0 on p
INTEGER :: jobnum = 0 ! Job number
@@ -28,7 +28,7 @@ MODULE basic
INTEGER :: my_id ! Rank in COMM_WORLD
INTEGER :: num_procs ! number of MPI processes
INTEGER :: num_procs_p ! Number of processes in p
- INTEGER :: num_procs_kx ! Number of processes in r
+ INTEGER :: num_procs_ky ! Number of processes in r
INTEGER :: num_procs_z ! Number of processes in z
INTEGER :: nbr_L, nbr_R ! Left and right neighbours (along p)
INTEGER :: nbr_T, nbr_B ! Top and bottom neighbours (along kx)
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index e0e66585b07938bdef6df1d0a31411fa31643249..b316e4425c0f496167abd87df7198268496f4cea 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -28,21 +28,21 @@ SUBROUTINE apply_closure_model
! all Napj s.t. p+2j <= 3
! -> (p,j) allowed are (0,0),(1,0),(0,1),(2,0),(1,1),(3,0)
! =>> Dmax is Pmax, condition is p+2j<=Pmax
+ DO iz = izs,ize
DO ikx = ikxs,ikxe
- DO iky = ikys,ikye
- DO iz = izs,ize
+ DO iky = ikys,ikye
IF(KIN_E) THEN
DO ip = ipgs_e,ipge_e
DO ij = ijgs_e,ijge_e
IF ( parray_e(ip)+2*jarray_e(ip) .GT. dmaxe) &
- moments_e(ip,ij,ikx,iky,iz,updatetlevel) = 0._dp
+ moments_e(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp
ENDDO
ENDDO
ENDIF
DO ip = ipgs_i,ipge_i
DO ij = ijgs_i,ijge_i
IF ( parray_i(ip)+2*jarray_i(ip) .GT. dmaxi) &
- moments_i(ip,ij,ikx,iky,iz,updatetlevel) = 0._dp
+ moments_i(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp
ENDDO
ENDDO
ENDDO
@@ -71,15 +71,15 @@ SUBROUTINE ghosts_upper_truncation
! applies only for the process that has largest j index
IF(ije_e .EQ. Jmaxe+1) THEN
DO ip = ipgs_e,ipge_e
- moments_e(ip,ije_e+1,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ moments_e(ip,ije_e+1,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
ENDDO
ENDIF
! applies only for the process that has largest p index
IF(ipe_e .EQ. Pmaxe+1) THEN
DO ij = ijgs_e,ijge_e
- moments_e(ipe_e+1,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ moments_e(ipe_e+1,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
IF(deltape .EQ. 1) THEN ! Must truncate the second stencil
- moments_e(ipe_e+2,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ moments_e(ipe_e+2,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
ENDIF
ENDDO
ENDIF
@@ -89,15 +89,15 @@ SUBROUTINE ghosts_upper_truncation
! applies only for the process that has largest j index
IF(ije_i .EQ. Jmaxi+1) THEN
DO ip = ipgs_i,ipge_i
- moments_i(ip,ije_i+1,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ moments_i(ip,ije_i+1,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
ENDDO
ENDIF
! applies only for the process that has largest p index
IF(ipe_i .EQ. Pmaxi+1) THEN
DO ij = ijgs_i,ijge_i
- moments_i(ipe_i+1,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ moments_i(ipe_i+1,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
IF(deltape .EQ. 1) THEN ! Must truncate the second stencil
- moments_i(ipe_i+2,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ moments_i(ipe_i+2,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
ENDIF
ENDDO
ENDIF
@@ -114,15 +114,15 @@ SUBROUTINE ghosts_lower_truncation
! applies only for the process that has lowest j index
IF(ijs_e .EQ. 1) THEN
DO ip = ipgs_e,ipge_e
- moments_e(ip,ijs_e-1,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ moments_e(ip,ijs_e-1,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
ENDDO
ENDIF
! applies only for the process that has lowest p index
IF(ips_e .EQ. 1) THEN
DO ij = ijgs_e,ijge_e
- moments_e(ips_e-1,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ moments_e(ips_e-1,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
IF(deltape .EQ. 1) THEN ! Must truncate the second stencil
- moments_e(ips_e-2,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ moments_e(ips_e-2,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
ENDIF
ENDDO
ENDIF
@@ -132,15 +132,15 @@ SUBROUTINE ghosts_lower_truncation
IF(ijs_i .EQ. 1) THEN
! applies only for the process that has lowest j index
DO ip = ipgs_i,ipge_i
- moments_i(ip,ijs_i-1,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ moments_i(ip,ijs_i-1,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
ENDDO
ENDIF
! applies only for the process that has lowest p index
IF(ips_i .EQ. 1) THEN
DO ij = ijgs_i,ijge_i
- moments_i(ips_i-1,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ moments_i(ips_i-1,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
IF(deltape .EQ. 1) THEN ! Must truncate the second stencil
- moments_i(ips_i-2,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ moments_i(ips_i-2,ij,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
ENDIF
ENDDO
ENDIF
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index 33954f78a6fecc715610def3db4fb0e243036092..ea43b1d3fc804b528b8e74fe5b54cd657ca7bfc3 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -117,15 +117,15 @@ CONTAINS
IF (KIN_E) THEN
DO ip = ips_e,ipe_e;DO ij = ijs_e,ije_e
DO ikx = ikxs, ikxe;DO iky = ikys, ikye; DO iz = izs,ize
- CALL LenardBernstein_e(ip,ij,ikx,iky,iz,TColl_)
- TColl_e(ip,ij,ikx,iky,iz) = TColl_
+ CALL LenardBernstein_e(ip,ij,iky,ikx,iz,TColl_)
+ TColl_e(ip,ij,iky,ikx,iz) = TColl_
ENDDO;ENDDO;ENDDO
ENDDO;ENDDO
ENDIF
DO ip = ips_i,ipe_i;DO ij = ijs_i,ije_i
DO ikx = ikxs, ikxe;DO iky = ikys, ikye; DO iz = izs,ize
- CALL LenardBernstein_i(ip,ij,ikx,iky,iz,TColl_)
- TColl_i(ip,ij,ikx,iky,iz) = TColl_
+ CALL LenardBernstein_i(ip,ij,iky,ikx,iz,TColl_)
+ TColl_i(ip,ij,iky,ikx,iz) = TColl_
ENDDO;ENDDO;ENDDO
ENDDO;ENDDO
END SUBROUTINE compute_lenard_bernstein
@@ -133,9 +133,9 @@ CONTAINS
!******************************************************************************!
!! for electrons
!******************************************************************************!
- SUBROUTINE LenardBernstein_e(ip_,ij_,ikx_,iky_,iz_,TColl_)
+ SUBROUTINE LenardBernstein_e(ip_,ij_,iky_,ikx_,iz_,TColl_)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_
+ INTEGER, INTENT(IN) :: ip_,ij_,iky_,ikx_,iz_
COMPLEX(dp), INTENT(OUT) :: TColl_
REAL(dp) :: j_dp, p_dp, be_2, kp
@@ -145,29 +145,29 @@ CONTAINS
eo_ = MODULO(parray_e(ip_),2)
p_dp = REAL(parray_e(ip_),dp)
j_dp = REAL(jarray_e(ij_),dp)
- kp = kparray(ikx_,iky_,iz_,eo_)
+ kp = kparray(iky_,ikx_,iz_,eo_)
be_2 = kp**2 * sigmae2_taue_o2 ! this is (be/2)^2
eo_ = MODULO(parray_e(ip_),2)
!** Assembling collison operator **
! -nuee (p + 2j) Nepj
- TColl_ = -nu_ee * (p_dp + 2._dp*j_dp)*nadiab_moments_e(ip_,ij_,ikx_,iky_,iz_)
+ TColl_ = -nu_ee * (p_dp + 2._dp*j_dp)*nadiab_moments_e(ip_,ij_,iky_,ikx_,iz_)
IF(gyrokin_CO) THEN
- TColl_ = TColl_ - nu_ee *2._dp*be_2*nadiab_moments_e(ip_,ij_,ikx_,iky_,iz_)
+ TColl_ = TColl_ - nu_ee *2._dp*be_2*nadiab_moments_e(ip_,ij_,iky_,ikx_,iz_)
ENDIF
END SUBROUTINE LenardBernstein_e
!******************************************************************************!
!! for ions
!******************************************************************************!
- SUBROUTINE LenardBernstein_i(ip_,ij_,ikx_,iky_,iz_,TColl_)
+ SUBROUTINE LenardBernstein_i(ip_,ij_,iky_,ikx_,iz_,TColl_)
USE fields, ONLY: moments_i
USE grid, ONLY: parray_i, jarray_i, kxarray, kyarray
USE basic
USE model, ONLY: sigmai2_taui_o2, nu_i
USE time_integration, ONLY : updatetlevel
IMPLICIT NONE
- INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_
+ INTEGER, INTENT(IN) :: ip_,ij_,iky_,ikx_,iz_
COMPLEX(dp), INTENT(OUT) :: TColl_
REAL(dp) :: j_dp, p_dp, kp, bi_2
@@ -177,14 +177,14 @@ CONTAINS
eo_ = MODULO(parray_i(ip_),2)
p_dp = REAL(parray_i(ip_),dp)
j_dp = REAL(jarray_i(ij_),dp)
- kp = kparray(ikx_,iky_,iz_,eo_)
+ kp = kparray(iky_,ikx_,iz_,eo_)
bi_2 = kp**2 * sigmai2_taui_o2 ! this is (be/2)^2
!** Assembling collison operator **
! -nuii (p + 2j) Nipj
- TColl_ = -nu_i * (p_dp + 2._dp*j_dp)*nadiab_moments_i(ip_,ij_,ikx_,iky_,iz_)
+ TColl_ = -nu_i * (p_dp + 2._dp*j_dp)*nadiab_moments_i(ip_,ij_,iky_,ikx_,iz_)
IF(gyrokin_CO) THEN
- TColl_ = TColl_ - nu_i *2._dp*bi_2*nadiab_moments_i(ip_,ij_,ikx_,iky_,iz_)
+ TColl_ = TColl_ - nu_i *2._dp*bi_2*nadiab_moments_i(ip_,ij_,iky_,ikx_,iz_)
ENDIF
END SUBROUTINE LenardBernstein_i
@@ -195,34 +195,40 @@ CONTAINS
IMPLICIT NONE
COMPLEX(dp) :: TColl_
IF (KIN_E) THEN
- DO ip = ips_e,ipe_e;DO ij = ijs_e,ije_e
- DO ikx = ikxs, ikxe;DO iky = ikys, ikye; DO iz = izs,ize
+ DO iz = izs,ize
+ DO iky = ikys, ikye;
+ DO ikx = ikxs, ikxe;
+ DO ij = ijs_e,ije_e
+ DO ip = ips_e,ipe_e;
IF(gyrokin_CO) THEN
- CALL DoughertyGK_ee(ip,ij,ikx,iky,iz,TColl_)
+ CALL DoughertyGK_ee(ip,ij,iky,ikx,iz,TColl_)
ELSE
- CALL DoughertyDK_ee(ip,ij,ikx,iky,iz,TColl_)
+ CALL DoughertyDK_ee(ip,ij,iky,ikx,iz,TColl_)
ENDIF
- TColl_e(ip,ij,ikx,iky,iz) = TColl_
+ TColl_e(ip,ij,iky,ikx,iz) = TColl_
ENDDO;ENDDO;ENDDO
ENDDO;ENDDO
ENDIF
- DO ip = ips_i,ipe_i;DO ij = ijs_i,ije_i
- DO ikx = ikxs, ikxe;DO iky = ikys, ikye; DO iz = izs,ize
+ DO iz = izs,ize
+ DO iky = ikys, ikye;
+ DO ikx = ikxs, ikxe;
+ DO ij = ijs_i,ije_i
+ DO ip = ips_i,ipe_i;
IF(gyrokin_CO) THEN
- CALL DoughertyGK_ii(ip,ij,ikx,iky,iz,TColl_)
+ CALL DoughertyGK_ii(ip,ij,iky,ikx,iz,TColl_)
ELSE
- CALL DoughertyDK_ii(ip,ij,ikx,iky,iz,TColl_)
+ CALL DoughertyDK_ii(ip,ij,iky,ikx,iz,TColl_)
ENDIF
- TColl_i(ip,ij,ikx,iky,iz) = TColl_
+ TColl_i(ip,ij,iky,ikx,iz) = TColl_
ENDDO;ENDDO;ENDDO
ENDDO;ENDDO
END SUBROUTINE compute_dougherty
!******************************************************************************!
!! Doughtery driftkinetic collision operator for electrons
!******************************************************************************!
- SUBROUTINE DoughertyDK_ee(ip_,ij_,ikx_,iky_,iz_,TColl_)
+ SUBROUTINE DoughertyDK_ee(ip_,ij_,iky_,ikx_,iz_,TColl_)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_
+ INTEGER, INTENT(IN) :: ip_,ij_,iky_,ikx_,iz_
COMPLEX(dp), INTENT(OUT) :: TColl_
COMPLEX(dp) :: upar, Tpar, Tperp
REAL(dp) :: j_dp, p_dp
@@ -230,24 +236,24 @@ CONTAINS
p_dp = REAL(parray_e(ip_),dp)
j_dp = REAL(jarray_e(ij_),dp)
!** Assembling collison operator **
- TColl_ = -(p_dp + 2._dp*j_dp)*nadiab_moments_e(ip_,ij_,ikx_,iky_,iz_)
+ TColl_ = -(p_dp + 2._dp*j_dp)*nadiab_moments_e(ip_,ij_,iky_,ikx_,iz_)
IF( (p_dp .EQ. 1._dp) .AND. (j_dp .EQ. 0._dp)) THEN !Ce10
- TColl_ = TColl_ + nadiab_moments_e(ip1_e,1,ikx_,iky_,iz_)
+ TColl_ = TColl_ + nadiab_moments_e(ip1_e,1,iky_,ikx_,iz_)
ELSEIF( (p_dp .EQ. 2._dp) .AND. (j_dp .EQ. 0._dp)) THEN ! Ce20
- TColl_ = TColl_ + twothird*nadiab_moments_e(ip2_e,1,ikx_,iky_,iz_) &
- - SQRT2*twothird*nadiab_moments_e(ip0_e,2,ikx_,iky_,iz_)
+ TColl_ = TColl_ + twothird*nadiab_moments_e(ip2_e,1,iky_,ikx_,iz_) &
+ - SQRT2*twothird*nadiab_moments_e(ip0_e,2,iky_,ikx_,iz_)
ELSEIF( (p_dp .EQ. 0._dp) .AND. (j_dp .EQ. 1._dp)) THEN ! Ce01
- TColl_ = TColl_ + 2._dp*twothird*nadiab_moments_e(ip0_e,2,ikx_,iky_,iz_) &
- - SQRT2*twothird*nadiab_moments_e(ip2_e,1,ikx_,iky_,iz_)
+ TColl_ = TColl_ + 2._dp*twothird*nadiab_moments_e(ip0_e,2,iky_,ikx_,iz_) &
+ - SQRT2*twothird*nadiab_moments_e(ip2_e,1,iky_,ikx_,iz_)
ENDIF
TColl_ = nu_ee * TColl_
END SUBROUTINE DoughertyDK_ee
!******************************************************************************!
!! Doughtery driftkinetic collision operator for ions
!******************************************************************************!
- SUBROUTINE DoughertyDK_ii(ip_,ij_,ikx_,iky_,iz_,TColl_)
+ SUBROUTINE DoughertyDK_ii(ip_,ij_,iky_,ikx_,iz_,TColl_)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_
+ INTEGER, INTENT(IN) :: ip_,ij_,iky_,ikx_,iz_
COMPLEX(dp), INTENT(OUT) :: TColl_
COMPLEX(dp) :: upar, Tpar, Tperp
REAL(dp) :: j_dp, p_dp
@@ -255,24 +261,24 @@ CONTAINS
p_dp = REAL(parray_i(ip_),dp)
j_dp = REAL(jarray_i(ij_),dp)
!** Assembling collison operator **
- TColl_ = -(p_dp + 2._dp*j_dp)*nadiab_moments_i(ip_,ij_,ikx_,iky_,iz_)
+ TColl_ = -(p_dp + 2._dp*j_dp)*nadiab_moments_i(ip_,ij_,iky_,ikx_,iz_)
IF( (p_dp .EQ. 1._dp) .AND. (j_dp .EQ. 0._dp)) THEN ! kronecker p1j0
- TColl_ = TColl_ + nadiab_moments_i(ip1_i,1,ikx_,iky_,iz_)
+ TColl_ = TColl_ + nadiab_moments_i(ip1_i,1,iky_,ikx_,iz_)
ELSEIF( (p_dp .EQ. 2._dp) .AND. (j_dp .EQ. 0._dp)) THEN ! kronecker p2j0
- TColl_ = TColl_ + twothird*nadiab_moments_i(ip2_i,1,ikx_,iky_,iz_) &
- - SQRT2*twothird*nadiab_moments_i(ip0_i,2,ikx_,iky_,iz_)
+ TColl_ = TColl_ + twothird*nadiab_moments_i(ip2_i,1,iky_,ikx_,iz_) &
+ - SQRT2*twothird*nadiab_moments_i(ip0_i,2,iky_,ikx_,iz_)
ELSEIF( (p_dp .EQ. 0._dp) .AND. (j_dp .EQ. 1._dp)) THEN ! kronecker p0j1
- TColl_ = TColl_ + 2._dp*twothird*nadiab_moments_i(ip0_i,2,ikx_,iky_,iz_) &
- - SQRT2*twothird*nadiab_moments_i(ip2_i,1,ikx_,iky_,iz_)
+ TColl_ = TColl_ + 2._dp*twothird*nadiab_moments_i(ip0_i,2,iky_,ikx_,iz_) &
+ - SQRT2*twothird*nadiab_moments_i(ip2_i,1,iky_,ikx_,iz_)
ENDIF
TColl_ = nu_i * TColl_
END SUBROUTINE DoughertyDK_ii
!******************************************************************************!
!! Doughtery gyrokinetic collision operator for electrons
!******************************************************************************!
- SUBROUTINE DoughertyGK_ee(ip_,ij_,ikx_,iky_,iz_,TColl_)
+ SUBROUTINE DoughertyGK_ee(ip_,ij_,iky_,ikx_,iz_,TColl_)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_
+ INTEGER, INTENT(IN) :: ip_,ij_,iky_,ikx_,iz_
COMPLEX(dp), INTENT(OUT) :: TColl_
COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_, T_
@@ -285,14 +291,14 @@ CONTAINS
p_dp = REAL(parray_e(ip_),dp)
eo_ = MODULO(parray_e(ip_),2)
j_dp = REAL(jarray_e(ij_),dp)
- kp = kparray(ikx_,iky_,iz_,eo_)
+ kp = kparray(iky_,ikx_,iz_,eo_)
be_2 = kp**2 * sigmae2_taue_o2 ! this is (be/2)^2
be_ = 2_dp*kp * sqrt_sigmae2_taue_o2 ! this is be
!** Assembling collison operator **
! Velocity-space diffusion (similar to Lenard Bernstein)
! -nuee (p + 2j + b^2/2) Nepj
- TColl_ = -(p_dp + 2._dp*j_dp + 2._dp*be_2)*nadiab_moments_e(ip_,ij_,ikx_,iky_,iz_)
+ TColl_ = -(p_dp + 2._dp*j_dp + 2._dp*be_2)*nadiab_moments_e(ip_,ij_,iky_,ikx_,iz_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
@@ -303,26 +309,26 @@ CONTAINS
DO in_ = 1,jmaxe+1
n_dp = REAL(in_-1,dp)
! Store the kernels for sparing readings
- Knp0 = Kernel_e(in_ ,ikx_,iky_,iz_,eo_)
- Knp1 = Kernel_e(in_+1,ikx_,iky_,iz_,eo_)
- Knm1 = Kernel_e(in_-1,ikx_,iky_,iz_,eo_)
+ Knp0 = Kernel_e(in_ ,iky_,ikx_,iz_,eo_)
+ Knp1 = Kernel_e(in_+1,iky_,ikx_,iz_,eo_)
+ Knm1 = Kernel_e(in_-1,iky_,ikx_,iz_,eo_)
! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = nadiab_moments_e(ip0_e,in_,ikx_,iky_,iz_)
+ nadiab_moment_0j = nadiab_moments_e(ip0_e,in_,iky_,ikx_,iz_)
! Density
n_ = n_ + Knp0 * nadiab_moment_0j
! Perpendicular velocity
uperp_ = uperp_ + be_*0.5_dp*(Knp0 - Knm1) * nadiab_moment_0j
! Parallel temperature
- Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_e(ip2_e,in_,ikx_,iky_,iz_) + nadiab_moment_0j)
+ Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_e(ip2_e,in_,iky_,ikx_,iz_) + nadiab_moment_0j)
! Perpendicular temperature
Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
ENDDO
T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
! Add energy restoring term
- TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_e(ij_ ,ikx_,iky_,iz_,eo_)
- TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_e(ij_+1,ikx_,iky_,iz_,eo_)
- TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_e(ij_-1,ikx_,iky_,iz_,eo_)
- TColl_ = TColl_ + uperp_*be_* (Kernel_e(ij_,ikx_,iky_,iz_,eo_) - Kernel_e(ij_-1,ikx_,iky_,iz_,eo_))
+ TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_e(ij_ ,iky_,ikx_,iz_,eo_)
+ TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_e(ij_+1,iky_,ikx_,iz_,eo_)
+ TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_e(ij_-1,iky_,ikx_,iz_,eo_)
+ TColl_ = TColl_ + uperp_*be_* (Kernel_e(ij_,iky_,ikx_,iz_,eo_) - Kernel_e(ij_-1,iky_,ikx_,iz_,eo_))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -331,9 +337,9 @@ CONTAINS
upar_ = 0._dp
DO in_ = 1,jmaxe+1
! Parallel velocity
- upar_ = upar_ + Kernel_e(in_,ikx_,iky_,iz_,eo_) * nadiab_moments_e(ip1_e,in_,ikx_,iky_,iz_)
+ upar_ = upar_ + Kernel_e(in_,iky_,ikx_,iz_,eo_) * nadiab_moments_e(ip1_e,in_,iky_,ikx_,iz_)
ENDDO
- TColl_ = TColl_ + upar_*Kernel_e(ij_,ikx_,iky_,iz_,eo_)
+ TColl_ = TColl_ + upar_*Kernel_e(ij_,iky_,ikx_,iz_,eo_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -345,20 +351,20 @@ CONTAINS
DO in_ = 1,jmaxe+1
n_dp = REAL(in_-1,dp)
! Store the kernels for sparing readings
- Knp0 = Kernel_e(in_ ,ikx_,iky_,iz_,eo_)
- Knp1 = Kernel_e(in_+1,ikx_,iky_,iz_,eo_)
- Knm1 = Kernel_e(in_-1,ikx_,iky_,iz_,eo_)
+ Knp0 = Kernel_e(in_ ,iky_,ikx_,iz_,eo_)
+ Knp1 = Kernel_e(in_+1,iky_,ikx_,iz_,eo_)
+ Knm1 = Kernel_e(in_-1,iky_,ikx_,iz_,eo_)
! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = nadiab_moments_e(ip0_e,in_,ikx_,iky_,iz_)
+ nadiab_moment_0j = nadiab_moments_e(ip0_e,in_,iky_,ikx_,iz_)
! Density
n_ = n_ + Knp0 * nadiab_moment_0j
! Parallel temperature
- Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_e(ip2_e,in_,ikx_,iky_,iz_) + nadiab_moment_0j)
+ Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_e(ip2_e,in_,iky_,ikx_,iz_) + nadiab_moment_0j)
! Perpendicular temperature
Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
ENDDO
T_ = (Tpar_ + 2._dp*Tperp_)/3._dp - n_
- TColl_ = TColl_ + T_*SQRT2*Kernel_e(ij_,ikx_,iky_,iz_,eo_)
+ TColl_ = TColl_ + T_*SQRT2*Kernel_e(ij_,iky_,ikx_,iz_,eo_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ENDIF
@@ -369,9 +375,9 @@ CONTAINS
!******************************************************************************!
!! Doughtery gyrokinetic collision operator for ions
!******************************************************************************!
- SUBROUTINE DoughertyGK_ii(ip_,ij_,ikx_,iky_,iz_,TColl_)
+ SUBROUTINE DoughertyGK_ii(ip_,ij_,iky_,ikx_,iz_,TColl_)
IMPLICIT NONE
- INTEGER, INTENT(IN) :: ip_,ij_,ikx_,iky_,iz_
+ INTEGER, INTENT(IN) :: ip_,ij_,iky_,ikx_,iz_
COMPLEX(dp), INTENT(OUT) :: TColl_
COMPLEX(dp) :: n_,upar_,uperp_,Tpar_, Tperp_, T_
@@ -385,14 +391,14 @@ CONTAINS
p_dp = REAL(parray_i(ip_),dp)
eo_ = MODULO(parray_i(ip_),2)
j_dp = REAL(jarray_i(ij_),dp)
- kp = kparray(ikx_,iky_,iz_,eo_)
+ kp = kparray(iky_,ikx_,iz_,eo_)
bi_2 = kp**2 *sigmai2_taui_o2 ! this is (bi/2)^2
bi_ = 2_dp*kp*sqrt_sigmai2_taui_o2 ! this is bi
!** Assembling collison operator **
! Velocity-space diffusion (similar to Lenard Bernstein)
! -nui (p + 2j + b^2/2) Nipj
- TColl_ = -(p_dp + 2._dp*j_dp + 2._dp*bi_2)*nadiab_moments_i(ip_,ij_,ikx_,iky_,iz_)
+ TColl_ = -(p_dp + 2._dp*j_dp + 2._dp*bi_2)*nadiab_moments_i(ip_,ij_,iky_,ikx_,iz_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( p_dp .EQ. 0 ) THEN ! Kronecker p0
@@ -403,26 +409,26 @@ CONTAINS
DO in_ = 1,jmaxi+1
n_dp = REAL(in_-1,dp)
! Store the kernels for sparing readings
- Knp0 = Kernel_i(in_ ,ikx_,iky_,iz_,eo_)
- Knp1 = Kernel_i(in_+1,ikx_,iky_,iz_,eo_)
- Knm1 = Kernel_i(in_-1,ikx_,iky_,iz_,eo_)
+ Knp0 = Kernel_i(in_ ,iky_,ikx_,iz_,eo_)
+ Knp1 = Kernel_i(in_+1,iky_,ikx_,iz_,eo_)
+ Knm1 = Kernel_i(in_-1,iky_,ikx_,iz_,eo_)
! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = nadiab_moments_i(ip0_i,in_,ikx_,iky_,iz_)
+ nadiab_moment_0j = nadiab_moments_i(ip0_i,in_,iky_,ikx_,iz_)
! Density
n_ = n_ + Knp0 * nadiab_moment_0j
! Perpendicular velocity
uperp_ = uperp_ + bi_*0.5_dp*(Knp0 - Knm1) * nadiab_moment_0j
! Parallel temperature
- Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_i(ip2_i,in_,ikx_,iky_,iz_) + nadiab_moment_0j)
+ Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_i(ip2_i,in_,iky_,ikx_,iz_) + nadiab_moment_0j)
! Perpendicular temperature
Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
ENDDO
T_ = (Tpar_ + 2._dp*Tperp_)*onethird - n_
! Add energy restoring term
- TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_i(ij_ ,ikx_,iky_,iz_,eo_)
- TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_i(ij_+1,ikx_,iky_,iz_,eo_)
- TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_i(ij_-1,ikx_,iky_,iz_,eo_)
- TColl_ = TColl_ + uperp_*bi_* (Kernel_i(ij_,ikx_,iky_,iz_,eo_) - Kernel_i(ij_-1,ikx_,iky_,iz_,eo_))
+ TColl_ = TColl_ + T_* 4._dp * j_dp * Kernel_i(ij_ ,iky_,ikx_,iz_,eo_)
+ TColl_ = TColl_ - T_* 2._dp * (j_dp + 1._dp) * Kernel_i(ij_+1,iky_,ikx_,iz_,eo_)
+ TColl_ = TColl_ - T_* 2._dp * j_dp * Kernel_i(ij_-1,iky_,ikx_,iz_,eo_)
+ TColl_ = TColl_ + uperp_*bi_* (Kernel_i(ij_,iky_,ikx_,iz_,eo_) - Kernel_i(ij_-1,iky_,ikx_,iz_,eo_))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -431,9 +437,9 @@ CONTAINS
upar_ = 0._dp
DO in_ = 1,jmaxi+1
! Parallel velocity
- upar_ = upar_ + Kernel_i(in_,ikx_,iky_,iz_,eo_) * nadiab_moments_i(ip1_i,in_,ikx_,iky_,iz_)
+ upar_ = upar_ + Kernel_i(in_,iky_,ikx_,iz_,eo_) * nadiab_moments_i(ip1_i,in_,iky_,ikx_,iz_)
ENDDO
- TColl_ = TColl_ + upar_*Kernel_i(ij_,ikx_,iky_,iz_,eo_)
+ TColl_ = TColl_ + upar_*Kernel_i(ij_,iky_,ikx_,iz_,eo_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Non zero term for p = 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -445,20 +451,20 @@ CONTAINS
DO in_ = 1,jmaxi+1
n_dp = REAL(in_-1,dp)
! Store the kernels for sparing readings
- Knp0 = Kernel_i(in_ ,ikx_,iky_,iz_,eo_)
- Knp1 = Kernel_i(in_+1,ikx_,iky_,iz_,eo_)
- Knm1 = Kernel_i(in_-1,ikx_,iky_,iz_,eo_)
+ Knp0 = Kernel_i(in_ ,iky_,ikx_,iz_,eo_)
+ Knp1 = Kernel_i(in_+1,iky_,ikx_,iz_,eo_)
+ Knm1 = Kernel_i(in_-1,iky_,ikx_,iz_,eo_)
! Nonadiabatic moments (only different from moments when p=0)
- nadiab_moment_0j = nadiab_moments_i(ip0_i,in_,ikx_,iky_,iz_)
+ nadiab_moment_0j = nadiab_moments_i(ip0_i,in_,iky_,ikx_,iz_)
! Density
n_ = n_ + Knp0 * nadiab_moment_0j
! Parallel temperature
- Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_i(ip2_i,in_,ikx_,iky_,iz_) + nadiab_moment_0j)
+ Tpar_ = Tpar_ + Knp0 * (SQRT2*nadiab_moments_i(ip2_i,in_,iky_,ikx_,iz_) + nadiab_moment_0j)
! Perpendicular temperature
Tperp_ = Tperp_ + ((2._dp*n_dp+1._dp)*Knp0 - (n_dp+1._dp)*Knp1 - n_dp*Knm1)*nadiab_moment_0j
ENDDO
T_ = (Tpar_ + 2._dp*Tperp_)*onethird - n_
- TColl_ = TColl_ + T_*SQRT2*Kernel_i(ij_,ikx_,iky_,iz_,eo_)
+ TColl_ = TColl_ + T_*SQRT2*Kernel_i(ij_,iky_,ikx_,iz_,eo_)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ENDIF
@@ -479,14 +485,14 @@ CONTAINS
COMPLEX(dp) :: TColl
INTEGER :: ikxs_C, ikxe_C, ikys_C, ikye_C
DO iz = izs,ize
- DO iky = ikys,ikye
- DO ikx = ikxs,ikxe
+ DO ikx = ikxs,ikxe
+ DO iky = ikys,ikye
IF(KIN_E) THEN
DO ij = 1,Jmaxe+1
! Electrons
! Loop over all p to compute sub collision term
DO ip = 1,total_np_e
- CALL apply_COSOlver_mat_e(ip,ij,ikx,iky,iz,TColl)
+ CALL apply_COSOlver_mat_e(ip,ij,iky,ikx,iz,TColl)
local_sum_e(ip) = TColl
ENDDO
IF (num_procs_p .GT. 1) THEN
@@ -501,14 +507,14 @@ CONTAINS
ENDIF
! Write in output variable
DO ip = ips_e,ipe_e
- TColl_e(ip,ij,ikx,iky,iz) = TColl_distr_e(ip)
+ TColl_e(ip,ij,iky,ikx,iz) = TColl_distr_e(ip)
ENDDO
ENDDO
ENDIF
! Ions
DO ij = 1,Jmaxi+1
DO ip = 1,total_np_i
- CALL apply_COSOlver_mat_i(ip,ij,ikx,iky,iz,TColl)
+ CALL apply_COSOlver_mat_i(ip,ij,iky,ikx,iz,TColl)
local_sum_i(ip) = TColl
ENDDO
IF (num_procs_p .GT. 1) THEN
@@ -524,7 +530,7 @@ CONTAINS
ENDIF
! Write in output variable
DO ip = ips_i,ipe_i
- TColl_i(ip,ij,ikx,iky,iz) = TColl_distr_i(ip)
+ TColl_i(ip,ij,iky,ikx,iz) = TColl_distr_i(ip)
ENDDO
ENDDO
ENDDO
@@ -535,13 +541,13 @@ CONTAINS
!******************************************************************************!
!!!!!!! Compute electron collision term
!******************************************************************************!
- SUBROUTINE apply_COSOlver_mat_e(ip_,ij_,ikx_,iky_,iz_,TColl_)
+ SUBROUTINE apply_COSOlver_mat_e(ip_,ij_,iky_,ikx_,iz_,TColl_)
IMPLICIT NONE
INTEGER, INTENT(IN) :: ip_, ij_ ,ikx_, iky_, iz_
COMPLEX(dp), INTENT(OUT) :: TColl_
- INTEGER :: ip2,ij2, p_int,j_int, p2_int,j2_int, ikx_C, iky_C, iz_C
+ INTEGER :: ip2,ij2, p_int,j_int, p2_int,j2_int, iky_C, ikx_C, iz_C
p_int = parray_e_full(ip_); j_int = jarray_e_full(ij_);
IF (gyrokin_CO) THEN ! GK operator (k-dependant)
@@ -558,9 +564,9 @@ CONTAINS
jloopee: DO ij2 = ijs_e,ije_e
j2_int = jarray_e(ij2)
IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxe))&
- TColl_ = TColl_ + nadiab_moments_e(ip2,ij2,ikx_,iky_,iz_) &
- *( nu_e * CeipjT(bare(p_int,j_int), bare(p2_int,j2_int),ikx_C, iky_C, iz_C) &
- +nu_ee * Ceepj (bare(p_int,j_int), bare(p2_int,j2_int),ikx_C, iky_C, iz_C))
+ TColl_ = TColl_ + nadiab_moments_e(ip2,ij2,iky_,ikx_,iz_) &
+ *( nu_e * CeipjT(bare(p_int,j_int), bare(p2_int,j2_int),iky_C, ikx_C, iz_C) &
+ +nu_ee * Ceepj (bare(p_int,j_int), bare(p2_int,j2_int),iky_C, ikx_C, iz_C))
ENDDO jloopee
ENDDO ploopee
@@ -570,8 +576,8 @@ CONTAINS
jloopei: DO ij2 = ijs_i,ije_i
j2_int = jarray_i(ij2)
IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxi))&
- TColl_ = TColl_ + nadiab_moments_i(ip2,ij2,ikx_,iky_,iz_) &
- *(nu_e * CeipjF(bare(p_int,j_int), bari(p2_int,j2_int),ikx_C, iky_C, iz_C))
+ TColl_ = TColl_ + nadiab_moments_i(ip2,ij2,iky_,ikx_,iz_) &
+ *(nu_e * CeipjF(bare(p_int,j_int), bari(p2_int,j2_int),iky_C, ikx_C, iz_C))
END DO jloopei
ENDDO ploopei
@@ -580,12 +586,12 @@ CONTAINS
!******************************************************************************!
!!!!!!! Compute ion collision term
!******************************************************************************!
- SUBROUTINE apply_COSOlver_mat_i(ip_,ij_,ikx_,iky_,iz_,TColl_)
+ SUBROUTINE apply_COSOlver_mat_i(ip_,ij_,iky_,ikx_,iz_,TColl_)
IMPLICIT NONE
INTEGER, INTENT(IN) :: ip_, ij_ ,ikx_, iky_, iz_
COMPLEX(dp), INTENT(OUT) :: TColl_
- INTEGER :: ip2,ij2, p_int,j_int, p2_int,j2_int, ikx_C, iky_C, iz_C
+ INTEGER :: ip2,ij2, p_int,j_int, p2_int,j2_int, iky_C, ikx_C, iz_C
p_int = parray_i_full(ip_); j_int = jarray_i_full(ij_);
IF (gyrokin_CO) THEN ! GK operator (k-dependant)
@@ -602,11 +608,11 @@ CONTAINS
j2_int = jarray_i(ij2)
IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxi))&
! Ion-ion collision
- TColl_ = TColl_ + nadiab_moments_i(ip2,ij2,ikx_,iky_,iz_) &
- * nu_i * Ciipj (bari(p_int,j_int), bari(p2_int,j2_int), ikx_C, iky_C, iz_C)
+ TColl_ = TColl_ + nadiab_moments_i(ip2,ij2,iky_,ikx_,iz_) &
+ * nu_i * Ciipj (bari(p_int,j_int), bari(p2_int,j2_int), iky_C, ikx_C, iz_C)
IF(KIN_E) & ! Ion-electron collision test
- TColl_ = TColl_ + nadiab_moments_i(ip2,ij2,ikx_,iky_,iz_) &
- * nu_ie * CiepjT(bari(p_int,j_int), bari(p2_int,j2_int), ikx_C, iky_C, iz_C)
+ TColl_ = TColl_ + nadiab_moments_i(ip2,ij2,iky_,ikx_,iz_) &
+ * nu_ie * CiepjT(bari(p_int,j_int), bari(p2_int,j2_int), iky_C, ikx_C, iz_C)
ENDDO jloopii
ENDDO ploopii
@@ -616,8 +622,8 @@ CONTAINS
jloopie: DO ij2 = ijs_e,ije_e
j2_int = jarray_e(ij2)
IF((CLOS .NE. 1) .OR. (p2_int+2*j2_int .LE. dmaxe))&
- TColl_ = TColl_ + nadiab_moments_e(ip2,ij2,ikx_,iky_,iz_) &
- *(nu_ie * CiepjF(bari(p_int,j_int), bare(p2_int,j2_int), ikx_C, iky_C, iz_C))
+ TColl_ = TColl_ + nadiab_moments_e(ip2,ij2,iky_,ikx_,iz_) &
+ *(nu_ie * CiepjF(bari(p_int,j_int), bare(p2_int,j2_int), iky_C, ikx_C, iz_C))
ENDDO jloopie
ENDDO ploopie
ENDIF
@@ -834,7 +840,7 @@ CONTAINS
DO iky = ikys,ikye
DO iz = izs,ize
! Check for nonlinear case if we are in the anti aliased domain or the filtered one
- kperp_sim = kparray(ikx,iky,iz,0) ! current simulation kperp
+ kperp_sim = kparray(iky,ikx,iz,0) ! current simulation kperp
! Find the interval in kp grid mat where kperp_sim is contained
! Loop over the whole kp mat grid to find the smallest kperp that is
! larger than the current kperp_sim (brute force...)
@@ -854,18 +860,18 @@ CONTAINS
! 0->1 variable for linear interp, i.e. zero2one = (k-k0)/(k1-k0)
zerotoone = (kperp_sim - kp_grid_mat(ikp_prev))/(kp_grid_mat(ikp_next) - kp_grid_mat(ikp_prev))
! Linear interpolation between previous and next kperp matrix values
- Ceepj (:,:,ikx,iky,iz) = (Ceepj__kp(:,:,ikp_next) - Ceepj__kp(:,:,ikp_prev))*zerotoone + Ceepj__kp(:,:,ikp_prev)
- Ciipj (:,:,ikx,iky,iz) = (Ciipj__kp(:,:,ikp_next) - Ciipj__kp(:,:,ikp_prev))*zerotoone + Ciipj__kp(:,:,ikp_prev)
+ Ceepj (:,:,iky,ikx,iz) = (Ceepj__kp(:,:,ikp_next) - Ceepj__kp(:,:,ikp_prev))*zerotoone + Ceepj__kp(:,:,ikp_prev)
+ Ciipj (:,:,iky,ikx,iz) = (Ciipj__kp(:,:,ikp_next) - Ciipj__kp(:,:,ikp_prev))*zerotoone + Ciipj__kp(:,:,ikp_prev)
IF(interspecies) THEN
- CeipjT(:,:,ikx,iky,iz) = (CeipjT_kp(:,:,ikp_next) - CeipjT_kp(:,:,ikp_prev))*zerotoone + CeipjT_kp(:,:,ikp_prev)
- CeipjF(:,:,ikx,iky,iz) = (CeipjF_kp(:,:,ikp_next) - CeipjF_kp(:,:,ikp_prev))*zerotoone + CeipjF_kp(:,:,ikp_prev)
- CiepjT(:,:,ikx,iky,iz) = (CiepjT_kp(:,:,ikp_next) - CiepjT_kp(:,:,ikp_prev))*zerotoone + CiepjT_kp(:,:,ikp_prev)
- CiepjF(:,:,ikx,iky,iz) = (CiepjF_kp(:,:,ikp_next) - CiepjF_kp(:,:,ikp_prev))*zerotoone + CiepjF_kp(:,:,ikp_prev)
+ CeipjT(:,:,iky,ikx,iz) = (CeipjT_kp(:,:,ikp_next) - CeipjT_kp(:,:,ikp_prev))*zerotoone + CeipjT_kp(:,:,ikp_prev)
+ CeipjF(:,:,iky,ikx,iz) = (CeipjF_kp(:,:,ikp_next) - CeipjF_kp(:,:,ikp_prev))*zerotoone + CeipjF_kp(:,:,ikp_prev)
+ CiepjT(:,:,iky,ikx,iz) = (CiepjT_kp(:,:,ikp_next) - CiepjT_kp(:,:,ikp_prev))*zerotoone + CiepjT_kp(:,:,ikp_prev)
+ CiepjF(:,:,iky,ikx,iz) = (CiepjF_kp(:,:,ikp_next) - CiepjF_kp(:,:,ikp_prev))*zerotoone + CiepjF_kp(:,:,ikp_prev)
ELSE
- CeipjT(:,:,ikx,iky,iz) = 0._dp
- CeipjF(:,:,ikx,iky,iz) = 0._dp
- CiepjT(:,:,ikx,iky,iz) = 0._dp
- CiepjF(:,:,ikx,iky,iz) = 0._dp
+ CeipjT(:,:,iky,ikx,iz) = 0._dp
+ CeipjF(:,:,iky,ikx,iz) = 0._dp
+ CiepjT(:,:,iky,ikx,iz) = 0._dp
+ CiepjF(:,:,iky,ikx,iz) = 0._dp
ENDIF
ENDDO
ENDDO
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index c5b45f0dd97a95be2ae6a2eaf559840ca7fad87b..925f515b039d666129b75cd6ccf1106e7b00aea5 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -88,8 +88,8 @@ SUBROUTINE diagnose(kstep)
CALL putarrnd(fidres, "/data/metric/gradzB", gradzB(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/Jacobian", Jacobian(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gradz_coeff", gradz_coeff(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/Ckxky", Ckxky(ikxs:ikxe,ikys:ikye,izs:ize,0:1), (/1, 1, 3/))
- CALL putarrnd(fidres, "/data/metric/kernel_i", kernel_i(ijs_i:ije_i,ikxs:ikxe,ikys:ikye,izs:ize,0:1), (/ 1, 2, 4/))
+ CALL putarrnd(fidres, "/data/metric/Ckxky", Ckxky(ikys:ikye,ikxs:ikxe,izs:ize,0:1), (/1, 1, 3/))
+ CALL putarrnd(fidres, "/data/metric/kernel_i", kernel_i(ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize,0:1), (/ 1, 2, 4/))
! var0d group (gyro transport)
IF (nsave_0d .GT. 0) THEN
@@ -217,7 +217,7 @@ SUBROUTINE diagnose(kstep)
CALL attach(fidres, TRIM(str), "cpu_time", -1)
CALL attach(fidres, TRIM(str), "Nproc", num_procs)
CALL attach(fidres, TRIM(str), "Np_p" , num_procs_p)
- CALL attach(fidres, TRIM(str), "Np_kx",num_procs_kx)
+ CALL attach(fidres, TRIM(str), "Np_kx",num_procs_ky)
CALL attach(fidres, TRIM(str), "write_gamma", write_gamma)
CALL attach(fidres, TRIM(str), "write_hf", write_hf)
CALL attach(fidres, TRIM(str), "write_phi", write_phi)
@@ -383,7 +383,7 @@ SUBROUTINE diagnose_2d
IMPLICIT NONE
- COMPLEX(dp) :: buffer(ikxs:ikxe,ikys:ikye)
+ COMPLEX(dp) :: buffer(ikys:ikye,ikxs:ikxe)
INTEGER :: i_, root, world_rank, world_size
CALL append(fidres, "/data/var2d/time", time,ionode=0)
@@ -398,24 +398,24 @@ CONTAINS
USE futils, ONLY: attach, putarr
USE grid, ONLY: ikxs,ikxe, ikys,ikye, Nkx, Nky, local_nkx
USE prec_const
- USE basic, ONLY : comm_kx, num_procs_p, rank_p
+ USE basic, ONLY : comm_ky, num_procs_p, rank_p
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(ikxs:ikxe, ikys:ikye), INTENT(IN) :: field
+ COMPLEX(dp), DIMENSION(ikys:ikye,ikxs:ikxe), INTENT(IN) :: field
CHARACTER(*), INTENT(IN) :: text
- COMPLEX(dp) :: buffer_dist(ikxs:ikxe,ikys:ikye)
- COMPLEX(dp) :: buffer_full(1:Nkx,1:Nky)
+ COMPLEX(dp) :: buffer_dist(ikys:ikye,ikxs:ikxe)
+ COMPLEX(dp) :: buffer_full(1:Nky,1:Nkx)
INTEGER :: scount, rcount
CHARACTER(LEN=50) :: dset_name
- scount = (ikxe-ikxs+1) * (ikye-ikys+1)
+ scount = (ikye-ikys+1) * (ikxe-ikxs+1)
rcount = scount
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var2d", TRIM(text), iframe2d
IF (num_procs .EQ. 1) THEN ! no data distribution
- CALL putarr(fidres, dset_name, field(ikxs:ikxe, ikys:ikye), ionode=0)
+ CALL putarr(fidres, dset_name, field(ikys:ikye,ikxs:ikxe), ionode=0)
ELSE
- CALL putarrnd(fidres, dset_name, field(ikxs:ikxe, ikys:ikye), (/1, 1/))
+ CALL putarrnd(fidres, dset_name, field(ikys:ikye,ikxs:ikxe), (/1, 1/))
ENDIF
CALL attach(fidres, dset_name, "time", time)
@@ -446,21 +446,21 @@ SUBROUTINE diagnose_3d
iframe3d=iframe3d+1
CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
- IF (write_phi) CALL write_field3d(phi (ikxs:ikxe,ikys:ikye,izs:ize), 'phi')
+ IF (write_phi) CALL write_field3d(phi (ikys:ikye,ikxs:ikxe,izs:ize), 'phi')
IF (write_Na00) THEN
IF(KIN_E)THEN
IF (ips_e .EQ. 1) THEN
- Ne00(ikxs:ikxe,ikys:ikye,izs:ize) = moments_e(ips_e,1,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel)
+ Ne00(ikys:ikye,ikxs:ikxe,izs:ize) = moments_e(ips_e,1,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel)
ENDIF
- ! CALL manual_3D_bcast(Ne00(ikxs:ikxe,ikys:ikye,izs:ize))
- CALL write_field3d(Ne00(ikxs:ikxe,ikys:ikye,izs:ize), 'Ne00')
+ ! CALL manual_3D_bcast(Ne00(ikys:ikye,ikxs:ikxe,izs:ize))
+ CALL write_field3d(Ne00(ikys:ikye,ikxs:ikxe,izs:ize), 'Ne00')
ENDIF
IF (ips_i .EQ. 1) THEN
- Ni00(ikxs:ikxe,ikys:ikye,izs:ize) = moments_i(ips_e,1,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel)
+ Ni00(ikys:ikye,ikxs:ikxe,izs:ize) = moments_i(ips_e,1,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel)
ENDIF
- ! CALL manual_3D_bcast(Ni00(ikxs:ikxe,ikys:ikye,izs:ize))
- CALL write_field3d(Ni00(ikxs:ikxe,ikys:ikye,izs:ize), 'Ni00')
+ ! CALL manual_3D_bcast(Ni00(ikys:ikye,ikxs:ikxe,izs:ize))
+ CALL write_field3d(Ni00(ikys:ikye,ikxs:ikxe,izs:ize), 'Ni00')
ENDIF
!! Fuid moments
@@ -469,29 +469,29 @@ SUBROUTINE diagnose_3d
IF (write_dens) THEN
IF(KIN_E)&
- CALL write_field3d(dens_e(ikxs:ikxe,ikys:ikye,izs:ize), 'dens_e')
- CALL write_field3d(dens_i(ikxs:ikxe,ikys:ikye,izs:ize), 'dens_i')
+ CALL write_field3d(dens_e(ikys:ikye,ikxs:ikxe,izs:ize), 'dens_e')
+ CALL write_field3d(dens_i(ikys:ikye,ikxs:ikxe,izs:ize), 'dens_i')
ENDIF
IF (write_fvel) THEN
IF(KIN_E)&
- CALL write_field3d(upar_e(ikxs:ikxe,ikys:ikye,izs:ize), 'upar_e')
- CALL write_field3d(upar_i(ikxs:ikxe,ikys:ikye,izs:ize), 'upar_i')
+ CALL write_field3d(upar_e(ikys:ikye,ikxs:ikxe,izs:ize), 'upar_e')
+ CALL write_field3d(upar_i(ikys:ikye,ikxs:ikxe,izs:ize), 'upar_i')
IF(KIN_E)&
- CALL write_field3d(uper_e(ikxs:ikxe,ikys:ikye,izs:ize), 'uper_e')
- CALL write_field3d(uper_i(ikxs:ikxe,ikys:ikye,izs:ize), 'uper_i')
+ CALL write_field3d(uper_e(ikys:ikye,ikxs:ikxe,izs:ize), 'uper_e')
+ CALL write_field3d(uper_i(ikys:ikye,ikxs:ikxe,izs:ize), 'uper_i')
ENDIF
IF (write_temp) THEN
IF(KIN_E)&
- CALL write_field3d(Tpar_e(ikxs:ikxe,ikys:ikye,izs:ize), 'Tpar_e')
- CALL write_field3d(Tpar_i(ikxs:ikxe,ikys:ikye,izs:ize), 'Tpar_i')
+ CALL write_field3d(Tpar_e(ikys:ikye,ikxs:ikxe,izs:ize), 'Tpar_e')
+ CALL write_field3d(Tpar_i(ikys:ikye,ikxs:ikxe,izs:ize), 'Tpar_i')
IF(KIN_E)&
- CALL write_field3d(Tper_e(ikxs:ikxe,ikys:ikye,izs:ize), 'Tper_e')
- CALL write_field3d(Tper_i(ikxs:ikxe,ikys:ikye,izs:ize), 'Tper_i')
+ CALL write_field3d(Tper_e(ikys:ikye,ikxs:ikxe,izs:ize), 'Tper_e')
+ CALL write_field3d(Tper_i(ikys:ikye,ikxs:ikxe,izs:ize), 'Tper_i')
IF(KIN_E)&
- CALL write_field3d(temp_e(ikxs:ikxe,ikys:ikye,izs:ize), 'temp_e')
- CALL write_field3d(temp_i(ikxs:ikxe,ikys:ikye,izs:ize), 'temp_i')
+ CALL write_field3d(temp_e(ikys:ikye,ikxs:ikxe,izs:ize), 'temp_e')
+ CALL write_field3d(temp_i(ikys:ikye,ikxs:ikxe,izs:ize), 'temp_i')
ENDIF
CONTAINS
@@ -500,18 +500,18 @@ CONTAINS
USE futils, ONLY: attach, putarr
USE grid, ONLY: ikxs,ikxe, ikys,ikye, Nkx, Nky, local_nkx
USE prec_const
- USE basic, ONLY : comm_kx, num_procs_p, rank_p
+ USE basic, ONLY : comm_ky, num_procs_p, rank_p
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(ikxs:ikxe, ikys:ikye, izs:ize), INTENT(IN) :: field
+ COMPLEX(dp), DIMENSION(ikys:ikye,ikxs:ikxe, izs:ize), INTENT(IN) :: field
CHARACTER(*), INTENT(IN) :: text
CHARACTER(LEN=50) :: dset_name
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var3d", TRIM(text), iframe3d
IF (num_procs .EQ. 1) THEN ! no data distribution
- CALL putarr(fidres, dset_name, field(ikxs:ikxe, ikys:ikye, izs:ize), ionode=0)
+ CALL putarr(fidres, dset_name, field(ikys:ikye,ikxs:ikxe, izs:ize), ionode=0)
ELSE
- CALL putarrnd(fidres, dset_name, field(ikxs:ikxe, ikys:ikye, izs:ize), (/1, 1, 3/))
+ CALL putarrnd(fidres, dset_name, field(ikys:ikye,ikxs:ikxe, izs:ize), (/1, 1, 3/))
ENDIF
CALL attach(fidres, dset_name, "time", time)
@@ -526,7 +526,8 @@ SUBROUTINE diagnose_5d
USE fields
USE array!, ONLY: Sepj, Sipj
USE grid, ONLY: ips_e,ipe_e, ips_i, ipe_i, &
- ijs_e,ije_e, ijs_i, ije_i
+ ijs_e,ije_e, ijs_i, ije_i, &
+ ikxs,ikxe,ikys,ikye,izs,ize
USE time_integration
USE diagnostics_par
USE prec_const
@@ -541,14 +542,14 @@ SUBROUTINE diagnose_5d
IF (write_Napj) THEN
IF(KIN_E)&
- CALL write_field5d_e(moments_e(ips_e:ipe_e,ijs_e:ije_e,:,:,:,updatetlevel), 'moments_e')
- CALL write_field5d_i(moments_i(ips_i:ipe_i,ijs_i:ije_i,:,:,:,updatetlevel), 'moments_i')
+ CALL write_field5d_e(moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel), 'moments_e')
+ CALL write_field5d_i(moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize,updatetlevel), 'moments_i')
ENDIF
IF (write_Sapj) THEN
IF(KIN_E)&
- CALL write_field5d_e(Sepj(ips_e:ipe_e,ijs_e:ije_e,:,:,:), 'Sepj')
- CALL write_field5d_i(Sipj(ips_i:ipe_i,ijs_i:ije_i,:,:,:), 'Sipj')
+ CALL write_field5d_e(Sepj(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), 'Sepj')
+ CALL write_field5d_i(Sipj(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), 'Sipj')
ENDIF
CONTAINS
@@ -559,16 +560,16 @@ SUBROUTINE diagnose_5d
USE prec_const
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,ikys:ikye,izs:ize), INTENT(IN) :: field
+ COMPLEX(dp), DIMENSION(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), INTENT(IN) :: field
CHARACTER(*), INTENT(IN) :: text
CHARACTER(LEN=50) :: dset_name
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
IF (num_procs .EQ. 1) THEN
- CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,ikys:ikye,izs:ize), ionode=0)
+ CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), ionode=0)
ELSE
- CALL putarrnd(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,ikys:ikye,izs:ize), (/1,3,5/))
+ CALL putarrnd(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikxs:ikxe,izs:ize), (/1,3,5/))
ENDIF
CALL attach(fidres, dset_name, 'cstep', cstep)
CALL attach(fidres, dset_name, 'time', time)
@@ -585,16 +586,16 @@ SUBROUTINE diagnose_5d
USE prec_const
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,ikys:ikye,izs:ize), INTENT(IN) :: field
+ COMPLEX(dp), DIMENSION(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), INTENT(IN) :: field
CHARACTER(*), INTENT(IN) :: text
CHARACTER(LEN=50) :: dset_name
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
IF (num_procs .EQ. 1) THEN
- CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,ikys:ikye,izs:ize), ionode=0)
+ CALL putarr(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), ionode=0)
ELSE
- CALL putarrnd(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,ikys:ikye,izs:ize), (/1,3,5/))
+ CALL putarrnd(fidres, dset_name, field(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikxs:ikxe,izs:ize), (/1,3,5/))
ENDIF
CALL attach(fidres, dset_name, 'cstep', cstep)
CALL attach(fidres, dset_name, 'time', time)
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index 5b947db46e95fbfcc79853df21afed50341d2d15..cdc2fefa87e990ffae94a1d50baaf3460b622e7c 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -20,7 +20,7 @@ MODULE fourier
type(C_PTR) , PUBLIC :: planf, planb
integer(C_INTPTR_T) :: i, ix, iy
integer(C_INTPTR_T), PUBLIC :: alloc_local_1, alloc_local_2
- integer(C_INTPTR_T) :: NX_, NY_, NX_halved
+ integer(C_INTPTR_T) :: NX_, NY_, NY_halved
integer :: communicator
! many plan data variables
integer(C_INTPTR_T) :: howmany=9 ! numer of eleemnt of the tensor
@@ -34,38 +34,38 @@ MODULE fourier
INTEGER, INTENT(IN) :: Nx,Ny
NX_ = Nx; NY_ = Ny
- NX_halved = NX_/2 + 1
+ NY_halved = NY_/2 + 1
! communicator = MPI_COMM_WORLD
- communicator = comm_kx
+ communicator = comm_ky
!! Complex arrays F, G
! Compute the room to allocate
- alloc_local_1 = fftw_mpi_local_size_2d(NX_halved, NY_, communicator, local_nkx, local_nkx_offset)
+ alloc_local_1 = fftw_mpi_local_size_2d(NY_halved, NX_, communicator, local_nky, local_nky_offset)
! Initalize pointers to this room
cdatac_f = fftw_alloc_complex(alloc_local_1)
cdatac_g = fftw_alloc_complex(alloc_local_1)
cdatac_c = fftw_alloc_complex(alloc_local_1)
! Initalize the arrays with the rooms pointed
- call c_f_pointer(cdatac_f, cmpx_data_f, [NY_ ,local_nkx])
- call c_f_pointer(cdatac_g, cmpx_data_g, [NY_ ,local_nkx])
- call c_f_pointer(cdatac_c, bracket_sum_c, [NY_ ,local_nkx])
+ call c_f_pointer(cdatac_f, cmpx_data_f, [NX_ ,local_nky])
+ call c_f_pointer(cdatac_g, cmpx_data_g, [NX_ ,local_nky])
+ call c_f_pointer(cdatac_c, bracket_sum_c, [NX_ ,local_nky])
!! Real arrays iFFT(F), iFFT(G)
! Compute the room to allocate
- alloc_local_2 = fftw_mpi_local_size_2d(NY_, NX_halved, communicator, local_nky, local_nky_offset)
+ alloc_local_2 = fftw_mpi_local_size_2d(NX_, NY_halved, communicator, local_nkx, local_nkx_offset)
! Initalize pointers to this room
cdatar_f = fftw_alloc_real(2*alloc_local_2)
cdatar_g = fftw_alloc_real(2*alloc_local_2)
cdatar_c = fftw_alloc_real(2*alloc_local_2)
! Initalize the arrays with the rooms pointed
- call c_f_pointer(cdatar_f, real_data_f, [2*(NX_/2 + 1),local_nky])
- call c_f_pointer(cdatar_g, real_data_g, [2*(NX_/2 + 1),local_nky])
- call c_f_pointer(cdatar_c, bracket_sum_r, [2*(NX_/2 + 1),local_nky])
+ call c_f_pointer(cdatar_f, real_data_f, [2*(NY_/2 + 1),local_nkx])
+ call c_f_pointer(cdatar_g, real_data_g, [2*(NY_/2 + 1),local_nkx])
+ call c_f_pointer(cdatar_c, bracket_sum_r, [2*(NY_/2 + 1),local_nkx])
! Plan Creation (out-of-place forward and backward FFT)
- planf = fftw_mpi_plan_dft_r2c_2D(NY_, NX_, real_data_f, cmpx_data_f, communicator, ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_OUT))
- planb = fftw_mpi_plan_dft_c2r_2D(NY_, NX_, cmpx_data_f, real_data_f, communicator, ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_IN))
+ planf = fftw_mpi_plan_dft_r2c_2D(NX_, NY_, real_data_f, cmpx_data_f, communicator, ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_OUT))
+ planb = fftw_mpi_plan_dft_c2r_2D(NX_, NY_, cmpx_data_f, real_data_f, communicator, ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_IN))
if ((.not. c_associated(planf)) .OR. (.not. c_associated(planb))) then
write(*,*) "plan creation error!!"
@@ -79,15 +79,15 @@ MODULE fourier
! - Compute the convolution using the convolution theorem
SUBROUTINE poisson_bracket_and_sum( F_, G_)
IMPLICIT NONE
- COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(ikxs:ikxe,ikys:ikye),&
+ COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(ikys:ikye,ikxs:ikxe),&
INTENT(IN) :: F_, G_ ! input fields
! First term df/dx x dg/dy
DO ikx = ikxs, ikxe
DO iky = ikys, ikye
- cmpx_data_f(iky,ikx-local_nkx_offset) = &
- imagu*kxarray(ikx)*F_(ikx,iky)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
- cmpx_data_g(iky,ikx-local_nkx_offset) = &
- imagu*kyarray(iky)*G_(ikx,iky)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ cmpx_data_f(ikx,iky-local_nky_offset) = &
+ imagu*kxarray(ikx)*F_(iky,ikx)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ cmpx_data_g(ikx,iky-local_nky_offset) = &
+ imagu*kyarray(iky)*G_(iky,ikx)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
ENDDO
ENDDO
call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
@@ -96,10 +96,10 @@ MODULE fourier
! Second term -df/dy x dg/dx
DO ikx = ikxs, ikxe
DO iky = ikys, ikye
- cmpx_data_f(iky,ikx-local_nkx_offset) = &
- imagu*kyarray(iky)*F_(ikx,iky)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
- cmpx_data_g(iky,ikx-local_nkx_offset) = &
- imagu*kxarray(ikx)*G_(ikx,iky)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ cmpx_data_f(ikx,iky-local_nky_offset) = &
+ imagu*kyarray(iky)*F_(iky,ikx)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ cmpx_data_g(ikx,iky-local_nky_offset) = &
+ imagu*kxarray(ikx)*G_(iky,ikx)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
ENDDO
ENDDO
call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
@@ -111,13 +111,13 @@ END SUBROUTINE poisson_bracket_and_sum
! - Compute the convolution using the convolution theorem
SUBROUTINE convolve_and_add( F_, G_)
IMPLICIT NONE
- COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(ikxs:ikxe,ikys:ikye),&
+ COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(ikys:ikye,ikxs:ikxe),&
INTENT(IN) :: F_, G_ ! input fields
! First term df/dx x dg/dy
DO ikx = ikxs, ikxe
DO iky = ikys, ikye
- cmpx_data_f(iky,ikx-local_nkx_offset) = F_(ikx,iky)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
- cmpx_data_g(iky,ikx-local_nkx_offset) = G_(ikx,iky)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ cmpx_data_f(ikx,iky-local_nky_offset) = F_(iky,ikx)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ cmpx_data_g(ikx,iky-local_nky_offset) = G_(iky,ikx)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
ENDDO
ENDDO
call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 9c8ff8097e06448e69e5835116a215f32fc2361e..f57942cc4c2deeb5125c95c420aa424b9a6b9b4f 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -86,7 +86,7 @@ CONTAINS
DO ikx = ikxs, ikxe
kx = kxarray(ikx)
DO iz = izgs,izge
- kparray(ikx, iky, iz, eo) = &
+ kparray(iky, ikx, iz, eo) = &
SQRT( gxx(iz,eo)*kx**2 + 2._dp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/hatB(iz,eo)
! there is a factor 1/B from the normalization; important to match GENE
ENDDO
@@ -209,7 +209,7 @@ CONTAINS
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
- Ckxky(ikx, iky, iz,eo) = - ky * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
+ Ckxky(iky, ikx, iz,eo) = - ky * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
! coefficient in the front of parallel derivative
@@ -266,7 +266,7 @@ CONTAINS
SUBROUTINE geometry_allocate_mem
! Curvature and geometry
- CALL allocate_array( Ckxky, ikxs,ikxe, ikys,ikye,izgs,izge,0,1)
+ CALL allocate_array( Ckxky, ikys,ikye, ikxs,ikxe,izgs,izge,0,1)
CALL allocate_array( Jacobian,izgs,izge, 0,1)
CALL allocate_array( gxx,izgs,izge, 0,1)
CALL allocate_array( gxy,izgs,izge, 0,1)
@@ -286,7 +286,7 @@ CONTAINS
CALL allocate_array( dxdR,izgs,izge, 0,1)
CALL allocate_array( dxdZ,izgs,izge, 0,1)
call allocate_array(gradz_coeff,izgs,izge, 0,1)
- CALL allocate_array( kparray, ikxs,ikxe, ikys,ikye,izgs,izge,0,1)
+ CALL allocate_array( kparray, ikys,ikye, ikxs,ikxe,izgs,izge,0,1)
END SUBROUTINE geometry_allocate_mem
diff --git a/src/ghosts_mod.F90 b/src/ghosts_mod.F90
index fee016af6f8b112176188be96c381bb3a580c8b9..59cf8ab72e3341fc606ab98160fe7c935449ba57 100644
--- a/src/ghosts_mod.F90
+++ b/src/ghosts_mod.F90
@@ -45,7 +45,7 @@ SUBROUTINE update_ghosts_p_e
IMPLICIT NONE
- count = (ijge_e-ijgs_e+1)*(ikxe-ikxs+1)*(ikye-ikys+1)*(izge-izgs+1)
+ count = (ijge_e-ijgs_e+1)*(ikye-ikys+1)*(ikxe-ikxs+1)*(izge-izgs+1)
!!!!!!!!!!! Send ghost to right neighbour !!!!!!!!!!!!!!!!!!!!!!
! Send the last local moment to fill the -1 neighbour ghost
@@ -73,7 +73,7 @@ SUBROUTINE update_ghosts_p_i
IMPLICIT NONE
- count = (ijge_i-ijgs_i+1)*(ikxe-ikxs+1)*(ikye-ikys+1)*(izge-izgs+1) ! Number of elements sent
+ count = (ijge_i-ijgs_i+1)*(ikye-ikys+1)*(ikxe-ikxs+1)*(izge-izgs+1) ! Number of elements sent
!!!!!!!!!!! Send ghost to right neighbour !!!!!!!!!!!!!!!!!!!!!!
CALL mpi_sendrecv(moments_i(ipe_i ,:,:,:,:,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 14, &
@@ -126,7 +126,7 @@ SUBROUTINE update_ghosts_z_e
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
IF (num_procs_z .GT. 1) THEN
- count = (ipge_e-ipgs_e+1)*(ijge_e-ijgs_e+1)*(ikxe-ikxs+1)*(ikye-ikys+1)
+ count = (ipge_e-ipgs_e+1)*(ijge_e-ijgs_e+1)*(ikye-ikys+1)*(ikxe-ikxs+1)
!!!!!!!!!!! Send ghost to up neighbour !!!!!!!!!!!!!!!!!!!!!!
! Send the last local moment to fill the -1 neighbour ghost
@@ -163,7 +163,7 @@ SUBROUTINE update_ghosts_z_i
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
IF (num_procs_z .GT. 1) THEN
- count = (ipge_i-ipgs_i+1)*(ijge_i-ijgs_i+1)*(ikxe-ikxs+1)*(ikye-ikys+1)
+ count = (ipge_i-ipgs_i+1)*(ijge_i-ijgs_i+1)*(ikye-ikys+1)*(ikxe-ikxs+1)
!!!!!!!!!!! Send ghost to up neighbour !!!!!!!!!!!!!!!!!!!!!!
! Send the last local moment to fill the -1 neighbour ghost
@@ -199,7 +199,7 @@ SUBROUTINE update_ghosts_z_phi
IF(Nz .GT. 1) THEN
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
IF (num_procs_z .GT. 1) THEN
- count = (ikxe-ikxs+1) * (ikye-ikys+1)
+ count = (ikye-ikys+1) * (ikxe-ikxs+1)
!!!!!!!!!!! Send ghost to up neighbour !!!!!!!!!!!!!!!!!!!!!!
! Send the last local moment to fill the -1 neighbour ghost
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 52889c105ddc87785ee789e476cae2d1db99acd5..50b9e5c76958ae06f4385e554ee3e2d17fdb9385 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -78,7 +78,7 @@ MODULE grid
INTEGER, PUBLIC :: ikx, iky, ip, ij, ikp, pp2, eo ! counters
LOGICAL, PUBLIC, PROTECTED :: contains_kx0 = .false. ! flag if the proc contains kx=0 index
LOGICAL, PUBLIC, PROTECTED :: contains_ky0 = .false. ! flag if the proc contains ky=0 index
- LOGICAL, PUBLIC, PROTECTED :: contains_kxmax = .false. ! flag if the proc contains kx=kxmax index
+ LOGICAL, PUBLIC, PROTECTED :: contains_kymax = .false. ! flag if the proc contains kx=kxmax index
! Grid containing the polynomials degrees
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: parray_e, parray_e_full
@@ -147,10 +147,10 @@ CONTAINS
SUBROUTINE init_1Dgrid_distr
! write(*,*) Nx
- local_nkx = (Nx/2+1)/num_procs_kx
+ local_nky = (Ny/2+1)/num_procs_ky
! write(*,*) local_nkx
- local_nkx_offset = rank_kx*local_nkx
- if (rank_kx .EQ. num_procs_kx-1) local_nkx = (Nx/2+1)-local_nkx_offset
+ local_nky_offset = rank_ky*local_nky
+ if (rank_ky .EQ. num_procs_ky-1) local_nky = (Ny/2+1)-local_nky_offset
END SUBROUTINE init_1Dgrid_distr
SUBROUTINE set_pgrid
@@ -271,136 +271,126 @@ CONTAINS
END SUBROUTINE set_jgrid
- SUBROUTINE set_kxgrid
+ SUBROUTINE set_kygrid
USE prec_const
USE model, ONLY: LINEARITY
IMPLICIT NONE
INTEGER :: i_
- Nkx = Nx/2+1 ! Defined only on positive kx since fields are real
+ Nky = Ny/2+1 ! Defined only on positive kx since fields are real
! Grid spacings
- IF (Nx .EQ. 1) THEN
- deltakx = 0._dp
- kx_max = 0._dp
- kx_min = 0._dp
+ IF (Ny .EQ. 1) THEN
+ deltaky = 0._dp
+ ky_max = 0._dp
+ ky_min = 0._dp
ELSE
- deltakx = 2._dp*PI/Lx
- kx_max = Nkx*deltakx
- kx_min = deltakx
+ deltaky = 2._dp*PI/Ly
+ ky_max = Nky*deltaky
+ ky_min = deltaky
ENDIF
! Build the full grids on process 0 to diagnose it without comm
- ALLOCATE(kxarray_full(1:Nkx))
- DO ikx = 1,Nkx
- kxarray_full(ikx) = REAL(ikx-1,dp) * deltakx
+ ALLOCATE(kyarray_full(1:Nky))
+ DO iky = 1,Nky
+ kyarray_full(iky) = REAL(iky-1,dp) * deltaky
END DO
!! Parallel distribution
- ikxs = local_nkx_offset + 1
- ikxe = ikxs + local_nkx - 1
- ALLOCATE(kxarray(ikxs:ikxe))
- local_kxmax = 0._dp
+ ikys = local_nky_offset + 1
+ ikye = ikys + local_nky - 1
+ ALLOCATE(kyarray(ikys:ikye))
+ local_kymax = 0._dp
! Creating a grid ordered as dk*(0 1 2 3)
- DO ikx = ikxs,ikxe
- kxarray(ikx) = REAL(ikx-1,dp) * deltakx
+ DO iky = ikys,ikye
+ kyarray(iky) = REAL(iky-1,dp) * deltaky
! Finding kx=0
- IF (kxarray(ikx) .EQ. 0) THEN
- ikx_0 = ikx
- contains_kx0 = .true.
+ IF (kyarray(iky) .EQ. 0) THEN
+ iky_0 = iky
+ contains_ky0 = .true.
ENDIF
! Finding local kxmax value
- IF (ABS(kxarray(ikx)) .GT. local_kxmax) THEN
- local_kxmax = ABS(kxarray(ikx))
+ IF (ABS(kyarray(iky)) .GT. local_kymax) THEN
+ local_kymax = ABS(kyarray(iky))
ENDIF
! Finding kxmax idx
- IF (kxarray(ikx) .EQ. kx_max) THEN
- ikx_max = ikx
- contains_kxmax = .true.
+ IF (kyarray(iky) .EQ. ky_max) THEN
+ iky_max = iky
+ contains_kymax = .true.
ENDIF
END DO
! Orszag 2/3 filter
- two_third_kxmax = 2._dp/3._dp*deltakx*(Nkx-1)
- ALLOCATE(AA_x(ikxs:ikxe))
- DO ikx = ikxs,ikxe
- IF ( (kxarray(ikx) .LT. two_third_kxmax) .OR. (LINEARITY .EQ. 'linear')) THEN
- AA_x(ikx) = 1._dp;
+ two_third_kymax = 2._dp/3._dp*deltaky*(Nky-1)
+ ALLOCATE(AA_y(ikys:ikye))
+ DO iky = ikys,ikye
+ IF ( (kyarray(iky) .LT. two_third_kymax) .OR. (LINEARITY .EQ. 'linear')) THEN
+ AA_y(iky) = 1._dp;
ELSE
- AA_x(ikx) = 0._dp;
+ AA_y(iky) = 0._dp;
ENDIF
END DO
- END SUBROUTINE set_kxgrid
+ END SUBROUTINE set_kygrid
- SUBROUTINE set_kygrid
+ SUBROUTINE set_kxgrid
USE prec_const
USE model, ONLY: LINEARITY
IMPLICIT NONE
INTEGER :: i_, counter
- Nky = Ny;
+ Nkx = Nx;
! Local data
! Start and END indices of grid
- ikys = 1
- ikye = Nky
- local_nky = ikye - ikys + 1
- ALLOCATE(kyarray(ikys:ikye))
- ALLOCATE(kyarray_full(1:Nky))
- IF (Ny .EQ. 1) THEN ! "cancel" y dimension
- deltaky = 1._dp
- kyarray(1) = 0._dp
- iky_0 = 1
- contains_ky0 = .true.
- ky_max = 0._dp
- iky_max = 1
- ky_min = 0._dp
- kyarray_full(1) = 0._dp
- local_kymax = 0._dp
+ ikxs = 1
+ ikxe = Nkx
+ local_nkx = ikxe - ikxs + 1
+ ALLOCATE(kxarray(ikxs:ikxe))
+ ALLOCATE(kxarray_full(1:Nkx))
+ IF (Nx .EQ. 1) THEN ! "cancel" y dimension
+ deltakx = 1._dp
+ kxarray(1) = 0._dp
+ ikx_0 = 1
+ contains_kx0 = .true.
+ kx_max = 0._dp
+ ikx_max = 1
+ kx_min = 0._dp
+ kxarray_full(1) = 0._dp
+ local_kxmax = 0._dp
ELSE ! Build apprpopriate grid
- deltaky = 2._dp*PI/Ly
- ky_max = (Ny/2)*deltakx
- ky_min = deltaky
+ deltakx = 2._dp*PI/Lx
+ kx_max = (Ny/2)*deltakx
+ kx_min = deltakx
! Creating a grid ordered as dk*(0 1 2 3 -2 -1)
- local_kymax = 0._dp
- DO iky = ikys,ikye
- SELECT CASE (LINEARITY)
- CASE ('linear') ! only positive freq for linear runs dk*(0 1 2 3 4 5)
- kyarray(iky) = deltaky*REAL(iky-1,dp)
- CASE DEFAULT !
- kyarray(iky) = deltaky*(MODULO(iky-1,Nky/2)-Nky/2*FLOOR(2.*real(iky-1)/real(Nky)))
- END SELECT
- if (iky .EQ. Ny/2+1) kyarray(iky) = -kyarray(iky)
- ! Finding ky=0
- IF (kyarray(iky) .EQ. 0) THEN
- iky_0 = iky
- contains_ky0 = .true.
+ local_kxmax = 0._dp
+ DO ikx = ikxs,ikxe
+ kxarray(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1)/real(Nkx)))
+ if (ikx .EQ. Nx/2+1) kxarray(ikx) = -kxarray(ikx)
+ ! Finding kx=0
+ IF (kxarray(ikx) .EQ. 0) THEN
+ ikx_0 = ikx
+ contains_kx0 = .true.
ENDIF
- ! Finding local kymax
- IF (ABS(kyarray(iky)) .GT. local_kymax) THEN
- local_kymax = ABS(kyarray(iky))
+ ! Finding local kxmax
+ IF (ABS(kxarray(ikx)) .GT. local_kxmax) THEN
+ local_kxmax = ABS(kxarray(ikx))
ENDIF
- ! Finding kymax
- IF (kyarray(iky) .EQ. ky_max) ikx_max = ikx
+ ! Finding kxmax
+ IF (kxarray(ikx) .EQ. kx_max) ikx_max = ikx
END DO
! Build the full grids on process 0 to diagnose it without comm
- ! ky
- DO iky = 1,Nky
- SELECT CASE (LINEARITY)
- CASE ('linear') ! only positive freq for linear runs dk*(0 1 2 3 4 5)
- kyarray_full(iky) = deltaky*REAL(iky-1,dp)
- CASE DEFAULT !
- kyarray_full(iky) = deltaky*(MODULO(iky-1,Nky/2)-Nky/2*FLOOR(2.*real(iky-1)/real(Nky)))
- IF (iky .EQ. Ny/2+1) kyarray_full(iky) = -kyarray_full(iky)
- END SELECT
+ ! kx
+ DO ikx = 1,Nkx
+ kxarray_full(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1)/real(Nkx)))
+ IF (ikx .EQ. Nx/2+1) kxarray_full(ikx) = -kxarray_full(ikx)
END DO
ENDIF
! Orszag 2/3 filter
- two_third_kymax = 2._dp/3._dp*deltaky*(Nky/2-1);
- ALLOCATE(AA_y(ikys:ikye))
- DO iky = ikys,ikye
- IF ( ((kyarray(iky) .GT. -two_third_kymax) .AND. &
- (kyarray(iky) .LT. two_third_kymax)) .OR. (LINEARITY .EQ. 'linear')) THEN
- AA_y(iky) = 1._dp;
+ two_third_kxmax = 2._dp/3._dp*deltakx*(Nkx/2-1);
+ ALLOCATE(AA_x(ikxs:ikxe))
+ DO ikx = ikxs,ikxe
+ IF ( ((kxarray(ikx) .GT. -two_third_kxmax) .AND. &
+ (kxarray(ikx) .LT. two_third_kxmax)) .OR. (LINEARITY .EQ. 'linear')) THEN
+ AA_x(ikx) = 1._dp;
ELSE
- AA_y(iky) = 0._dp;
+ AA_x(ikx) = 0._dp;
ENDIF
END DO
- END SUBROUTINE set_kygrid
+ END SUBROUTINE set_kxgrid
SUBROUTINE set_zgrid
diff --git a/src/inital.F90 b/src/inital.F90
index ae42fbfbdef65dfb21dce0f7359576a5d6cda15e..aacfd20a0586eecf71f3db0131749ef1545f2c0c 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -124,13 +124,13 @@ SUBROUTINE init_moments
DO iky=ikys,ikye
DO iz=izs,ize
CALL RANDOM_NUMBER(noise)
- moments_e(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
+ moments_e(ip,ij,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
END DO
END DO
END DO
- IF ( contains_kx0 ) THEN
- DO iky=2,Nky/2 !symmetry at kx = 0 for all z
- moments_e(ip,ij,ikx_0,iky,:,:) = moments_e( ip,ij,ikx_0,Nky+2-iky,:, :)
+ IF ( contains_ky0 ) THEN
+ DO iky=2,Nky/2 !symmetry at ky = 0 for all z
+ moments_e(ip,ij,iky_0,ikx,:,:) = moments_e( ip,ij,iky_0,Nkx+2-ikx,:, :)
END DO
ENDIF
END DO
@@ -144,14 +144,14 @@ SUBROUTINE init_moments
DO iky=ikys,ikye
DO iz=izs,ize
CALL RANDOM_NUMBER(noise)
- moments_i(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
+ moments_i(ip,ij,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
END DO
END DO
END DO
- IF ( contains_kx0 ) THEN
- DO iky=2,Nky/2 !symmetry at kx = 0 for all z
- moments_i( ip,ij,ikx_0,iky,:,:) = moments_i( ip,ij,ikx_0,Nky+2-iky,:,:)
+ IF ( contains_ky0 ) THEN
+ DO ikx=2,Nkx/2 !symmetry at ky = 0 for all z
+ moments_i( ip,ij,iky_0,ikx,:,:) = moments_i( ip,ij,iky_0,Nkx+2-ikx,:,:)
END DO
ENDIF
@@ -166,13 +166,13 @@ SUBROUTINE init_moments
IF(KIN_E) THEN
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
- moments_e( ip,ij,ikx,iky,iz, :) = moments_e( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
+ moments_e( ip,ij,iky,ikx,iz, :) = moments_e( ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
ENDIF
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
- moments_i( ip,ij,ikx,iky,iz, :) = moments_i( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
+ moments_i( ip,ij,iky,ikx,iz, :) = moments_i( ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
ENDDO
@@ -212,16 +212,16 @@ SUBROUTINE init_gyrodens
DO iz=izs,ize
CALL RANDOM_NUMBER(noise)
IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
- moments_e(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
+ moments_e(ip,ij,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
ELSE
- moments_e(ip,ij,ikx,iky,iz,:) = 0._dp
+ moments_e(ip,ij,iky,ikx,iz,:) = 0._dp
ENDIF
END DO
END DO
END DO
- IF ( contains_kx0 ) THEN
- DO iky=2,Nky/2 !symmetry at kx = 0 for all z
- moments_e(ip,ij,ikx_0,iky,:,:) = moments_e( ip,ij,ikx_0,Nky+2-iky,:, :)
+ IF ( contains_ky0 ) THEN
+ DO ikx=2,Nkx/2 !symmetry at ky = 0 for all z
+ moments_e(ip,ij,iky_0,ikx,:,:) = moments_e( ip,ij,iky_0,Nkx+2-ikx,:, :)
END DO
ENDIF
END DO
@@ -235,16 +235,16 @@ SUBROUTINE init_gyrodens
DO iz=izs,ize
CALL RANDOM_NUMBER(noise)
IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
- moments_i(ip,ij,ikx,iky,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
+ moments_i(ip,ij,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
ELSE
- moments_i(ip,ij,ikx,iky,iz,:) = 0._dp
+ moments_i(ip,ij,iky,ikx,iz,:) = 0._dp
ENDIF
END DO
END DO
END DO
- IF ( contains_kx0 ) THEN
- DO iky=2,Nky/2 !symmetry at kx = 0 for all z
- moments_i( ip,ij,ikx_0,iky,:,:) = moments_i( ip,ij,ikx_0,Nky+2-iky,:,:)
+ IF ( contains_ky0 ) THEN
+ DO iky=2,Nky/2 !symmetry at ky = 0 for all z
+ moments_i( ip,ij,iky_0,ikx,:,:) = moments_i( ip,ij,iky_0,Nkx+2-ikx,:,:)
END DO
ENDIF
END DO
@@ -258,13 +258,13 @@ SUBROUTINE init_gyrodens
IF(KIN_E) THEN
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
- moments_e( ip,ij,ikx,iky,iz, :) = moments_e( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
+ moments_e( ip,ij,iky,ikx,iz, :) = moments_e( ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
ENDIF
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
- moments_i( ip,ij,ikx,iky,iz, :) = moments_i( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
+ moments_i( ip,ij,iky,ikx,iz, :) = moments_i( ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
ENDDO
@@ -300,17 +300,17 @@ SUBROUTINE init_phi
DO iky=ikys,ikye
DO iz=izs,ize
CALL RANDOM_NUMBER(noise)
- phi(ikx,iky,iz) = (init_background + init_noiselvl*(noise-0.5_dp))!*AA_x(ikx)*AA_y(iky)
+ phi(iky,ikx,iz) = (init_background + init_noiselvl*(noise-0.5_dp))!*AA_x(ikx)*AA_y(iky)
ENDDO
END DO
END DO
- !symmetry at kx = 0 to keep real inverse transform
- IF ( contains_kx0 ) THEN
- DO iky=2,Nky/2
- phi(ikx_0,iky,izs:ize) = phi(ikx_0,Nky+2-iky,izs:ize)
+ !symmetry at ky = 0 to keep real inverse transform
+ IF ( contains_ky0 ) THEN
+ DO ikx=2,Nkx/2
+ phi(iky_0,ikx,izs:ize) = phi(iky_0,Nkx+2-ikx,izs:ize)
END DO
- phi(ikx_0,Ny/2,izs:ize) = REAL(phi(ikx_0,Ny/2,izs:ize)) !origin must be real
+ phi(iky_0,ikx_0,izs:ize) = REAL(phi(iky_0,ikx_0,izs:ize)) !origin must be real
ENDIF
!**** ensure no previous moments initialization
@@ -323,9 +323,9 @@ SUBROUTINE init_phi
IF (my_id .EQ. 0) WRITE(*,*) 'Init ZF phi'
IF( (INIT_ZF+1 .GT. ikxs) .AND. (INIT_ZF+1 .LT. ikxe) ) THEN
DO iz = izs,ize
- phi(INIT_ZF+1,iky_0,iz) = ZF_AMP*(2._dp*PI)**2/deltakx/deltaky/2._dp * COS((iz-1)/Nz*2._dp*PI)
- moments_i(1,1,INIT_ZF+1,iky_0,iz,:) = kxarray(INIT_ZF+1)**2*phi(INIT_ZF+1,iky_0,iz)* COS((iz-1)/Nz*2._dp*PI)
- IF(KIN_E) moments_e(1,1,INIT_ZF+1,iky_0,iz,:) = 0._dp
+ phi(iky_0,INIT_ZF+1,iz) = ZF_AMP*(2._dp*PI)**2/deltakx/deltaky/2._dp * COS((iz-1)/Nz*2._dp*PI)
+ moments_i(1,1,iky_0,INIT_ZF+1,iz,:) = kxarray(INIT_ZF+1)**2*phi(iky_0,INIT_ZF+1,iz)* COS((iz-1)/Nz*2._dp*PI)
+ IF(KIN_E) moments_e(1,1,iky_0,INIT_ZF+1,iz,:) = 0._dp
ENDDO
ENDIF
ENDIF
@@ -354,7 +354,7 @@ SUBROUTINE initialize_blob
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
IF( (iky .NE. iky_0) .AND. (ip .EQ. 1) .AND. (ij .EQ. 1)) THEN
- moments_i( ip,ij,ikx,iky,iz, :) = moments_i( ip,ij,ikx,iky,iz, :) &
+ moments_i( ip,ij,iky,ikx,iz, :) = moments_i( ip,ij,iky,ikx,iz, :) &
+ gain*sigma/SQRT2 * exp(-(kx**2+ky**2)*sigma**2/4._dp) &
* AA_x(ikx)*AA_y(iky)!&
! * exp(sigmai2_taui_o2*(kx**2+ky**2))
@@ -365,7 +365,7 @@ SUBROUTINE initialize_blob
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
IF( (iky .NE. iky_0) .AND. (ip .EQ. 1) .AND. (ij .EQ. 1)) THEN
- moments_e( ip,ij,ikx,iky,iz, :) = moments_e( ip,ij,ikx,iky,iz, :) &
+ moments_e( ip,ij,iky,ikx,iz, :) = moments_e( ip,ij,iky,ikx,iz, :) &
+ gain*sigma/SQRT2 * exp(-(kx**2+ky**2)*sigma**2/4._dp) &
* AA_x(ikx)*AA_y(iky)!&
! * exp(sigmai2_taui_o2*(kx**2+ky**2))
@@ -416,28 +416,28 @@ SUBROUTINE init_ppj
z = zarray(iz,0)
IF (kx .EQ. 0) THEN
IF(ky .EQ. 0) THEN
- moments_e(ip,ij,ikx,iky,iz,:) = 0._dp
+ moments_e(ip,ij,iky,ikx,iz,:) = 0._dp
ELSE
- moments_e(ip,ij,ikx,iky,iz,:) = 0.5_dp * ky_min/(ABS(ky)+ky_min)
+ moments_e(ip,ij,iky,ikx,iz,:) = 0.5_dp * ky_min/(ABS(ky)+ky_min)
ENDIF
ELSE
IF(ky .GT. 0) THEN
- moments_e(ip,ij,ikx,iky,iz,:) = (kx_min/(ABS(kx)+kx_min))*(ky_min/(ABS(ky)+ky_min))
+ moments_e(ip,ij,iky,ikx,iz,:) = (kx_min/(ABS(kx)+kx_min))*(ky_min/(ABS(ky)+ky_min))
ELSE
- moments_e(ip,ij,ikx,iky,iz,:) = 0.5_dp*(kx_min/(ABS(kx)+kx_min))
+ moments_e(ip,ij,iky,ikx,iz,:) = 0.5_dp*(kx_min/(ABS(kx)+kx_min))
ENDIF
ENDIF
! z-dep
- moments_e(ip,ij,ikx,iky,iz,:) = moments_e(ip,ij,ikx,iky,iz,:) * &
+ moments_e(ip,ij,iky,ikx,iz,:) = moments_e(ip,ij,iky,ikx,iz,:) * &
! (1 + exp(-(z/sigma_z)**2/2.0)*sqrt(2.0*sqrt(pi)/sigma_z))
(Jacobian(iz,0)*iInt_Jacobian)**2
END DO
END DO
END DO
- IF ( contains_kx0 ) THEN
- DO iky=2,Nky/2 !symmetry at kx = 0 for all z
- moments_e(ip,ij,ikx_0,iky,:,:) = moments_e( ip,ij,ikx_0,Nky+2-iky,:, :)
+ IF ( contains_ky0 ) THEN
+ DO ikx=2,Nkx/2 !symmetry at kx = 0 for all z
+ moments_e(ip,ij,iky_0,ikx,:,:) = moments_e( ip,ij,iky_0,Nkx+2-ikx,:, :)
END DO
ENDIF
ELSE
@@ -459,28 +459,28 @@ SUBROUTINE init_ppj
z = zarray(iz,0)
IF (kx .EQ. 0) THEN
IF(ky .EQ. 0) THEN
- moments_i(ip,ij,ikx,iky,iz,:) = 0._dp
+ moments_i(ip,ij,iky,ikx,iz,:) = 0._dp
ELSE
- moments_i(ip,ij,ikx,iky,iz,:) = 0.5_dp * ky_min/(ABS(ky)+ky_min)
+ moments_i(ip,ij,iky,ikx,iz,:) = 0.5_dp * ky_min/(ABS(ky)+ky_min)
ENDIF
ELSE
IF(ky .GT. 0) THEN
- moments_i(ip,ij,ikx,iky,iz,:) = (kx_min/(ABS(kx)+kx_min))*(ky_min/(ABS(ky)+ky_min))
+ moments_i(ip,ij,iky,ikx,iz,:) = (kx_min/(ABS(kx)+kx_min))*(ky_min/(ABS(ky)+ky_min))
ELSE
- moments_i(ip,ij,ikx,iky,iz,:) = 0.5_dp*(kx_min/(ABS(kx)+kx_min))
+ moments_i(ip,ij,iky,ikx,iz,:) = 0.5_dp*(kx_min/(ABS(kx)+kx_min))
ENDIF
ENDIF
! z-dep
- moments_i(ip,ij,ikx,iky,iz,:) = moments_i(ip,ij,ikx,iky,iz,:) * &
+ moments_i(ip,ij,iky,ikx,iz,:) = moments_i(ip,ij,iky,ikx,iz,:) * &
! (1 + exp(-(z/sigma_z)**2/2.0)*sqrt(2.0*sqrt(pi)/sigma_z))
(Jacobian(iz,0)*iInt_Jacobian)**2
END DO
END DO
END DO
- IF ( contains_kx0 ) THEN
- DO iky=2,Nky/2 !symmetry at kx = 0 for all z
- moments_i( ip,ij,ikx_0,iky,:,:) = moments_i( ip,ij,ikx_0,Nky+2-iky,:,:)
+ IF ( contains_ky0 ) THEN
+ DO ikx=2,Nkx/2 !symmetry at kx = 0 for all z
+ moments_i(ip,ij,iky_0,ikx,:,:) = moments_i( ip,ij,iky_0,Nkx+2-ikx,:, :)
END DO
ENDIF
ELSE
@@ -496,12 +496,12 @@ SUBROUTINE init_ppj
DO iz=izs,ize
DO ip=ips_e,ipe_e
DO ij=ijs_e,ije_e
- moments_e( ip,ij,ikx,iky,iz, :) = moments_e( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
+ moments_e( ip,ij,iky,ikx,iz, :) = moments_e( ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
DO ip=ips_i,ipe_i
DO ij=ijs_i,ije_i
- moments_i( ip,ij,ikx,iky,iz, :) = moments_i( ip,ij,ikx,iky,iz, :)*AA_x(ikx)*AA_y(iky)
+ moments_i( ip,ij,iky,ikx,iz, :) = moments_i( ip,ij,iky,ikx,iz, :)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
ENDDO
diff --git a/src/memory.F90 b/src/memory.F90
index a90fec8c24722e6bfb299bb2909f1e9ee04a08ee..a1ebae20539ab897cb0e0887645321db8ac0b196 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -13,31 +13,31 @@ SUBROUTINE memory
IMPLICIT NONE
! Electrostatic potential
- CALL allocate_array( phi, ikxs,ikxe, ikys,ikye, izgs,izge)
+ CALL allocate_array( phi, ikys,ikye, ikxs,ikxe, izgs,izge)
CALL allocate_array( phi_ZF, ikxs,ikxe, izs,ize)
CALL allocate_array( phi_EM, ikys,ikye, izs,ize)
- CALL allocate_array(inv_poisson_op, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array(inv_poisson_op, ikys,ikye, ikxs,ikxe, izs,ize)
!Electrons arrays
IF(KIN_E) THEN
- CALL allocate_array( Ne00, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( dens_e, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( upar_e, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( uper_e, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( Tpar_e, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( Tper_e, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( temp_e, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( Kernel_e, ijgs_e,ijge_e, ikxs,ikxe, ikys,ikye, izgs,izge, 0,1)
- CALL allocate_array( moments_e, ipgs_e,ipge_e, ijgs_e,ijge_e, ikxs,ikxe, ikys,ikye, izgs,izge, 1,ntimelevel )
- CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
- CALL allocate_array( nadiab_moments_e, ipgs_e,ipge_e, ijgs_e,ijge_e, ikxs,ikxe, ikys,ikye, izgs,izge)
- CALL allocate_array( ddz_nepj, ipgs_e,ipge_e, ijgs_e,ijge_e, ikxs,ikxe, ikys,ikye, izgs,izge)
- CALL allocate_array( ddz2_Nepj, ipgs_e,ipge_e, ijgs_e,ijge_e, ikxs,ikxe, ikys,ikye, izgs,izge)
- CALL allocate_array( interp_nepj, ipgs_e,ipge_e, ijgs_e,ijge_e, ikxs,ikxe, ikys,ikye, izgs,izge)
+ CALL allocate_array( Ne00, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( dens_e, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( upar_e, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( uper_e, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( Tpar_e, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( Tper_e, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( temp_e, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( Kernel_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge, 0,1)
+ CALL allocate_array( moments_e, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge, 1,ntimelevel )
+ CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikys,ikye, ikxs,ikxe, izs,ize, 1,ntimelevel )
+ CALL allocate_array( nadiab_moments_e, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge)
+ CALL allocate_array( ddz_nepj, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge)
+ CALL allocate_array( ddz2_Nepj, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge)
+ CALL allocate_array( interp_nepj, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, ikxs,ikxe, izgs,izge)
CALL allocate_array( moments_e_ZF, ipgs_e,ipge_e, ijgs_e,ijge_e, ikxs,ikxe, izs,ize)
CALL allocate_array( moments_e_EM, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, izs,ize)
- CALL allocate_array( TColl_e, ips_e,ipe_e, ijs_e,ije_e , ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( Sepj, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( TColl_e, ips_e,ipe_e, ijs_e,ije_e , ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( Sepj, ips_e,ipe_e, ijs_e,ije_e, ikys,ikye, ikxs,ikxe, izs,ize)
CALL allocate_array( xnepj, ips_e,ipe_e, ijs_e,ije_e)
CALL allocate_array( xnepp2j, ips_e,ipe_e)
CALL allocate_array( xnepp1j, ips_e,ipe_e)
@@ -55,24 +55,24 @@ SUBROUTINE memory
ENDIF
!Ions arrays
- CALL allocate_array( Ni00, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( dens_i, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( upar_i, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( uper_i, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( Tpar_i, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( Tper_i, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( temp_i, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( Kernel_i, ijgs_i,ijge_i, ikxs,ikxe, ikys,ikye, izgs,izge, 0,1)
- CALL allocate_array( moments_i, ipgs_i,ipge_i, ijgs_i,ijge_i, ikxs,ikxe, ikys,ikye, izgs,izge, 1,ntimelevel )
- CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
- CALL allocate_array( nadiab_moments_i, ipgs_i,ipge_i, ijgs_i,ijge_i, ikxs,ikxe, ikys,ikye, izgs,izge)
- CALL allocate_array( ddz_nipj, ipgs_i,ipge_i, ijgs_i,ijge_i, ikxs,ikxe, ikys,ikye, izgs,izge)
- CALL allocate_array( ddz2_Nipj, ipgs_i,ipge_i, ijgs_i,ijge_i, ikxs,ikxe, ikys,ikye, izgs,izge)
- CALL allocate_array( interp_nipj, ipgs_i,ipge_i, ijgs_i,ijge_i, ikxs,ikxe, ikys,ikye, izgs,izge)
+ CALL allocate_array( Ni00, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( dens_i, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( upar_i, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( uper_i, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( Tpar_i, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( Tper_i, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( temp_i, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( Kernel_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge, 0,1)
+ CALL allocate_array( moments_i, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge, 1,ntimelevel )
+ CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikys,ikye, ikxs,ikxe, izs,ize, 1,ntimelevel )
+ CALL allocate_array( nadiab_moments_i, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge)
+ CALL allocate_array( ddz_nipj, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge)
+ CALL allocate_array( ddz2_Nipj, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge)
+ CALL allocate_array( interp_nipj, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, ikxs,ikxe, izgs,izge)
CALL allocate_array( moments_i_ZF, ipgs_i,ipge_i, ijgs_i,ijge_i, ikxs,ikxe, izs,ize)
CALL allocate_array( moments_i_EM, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, izs,ize)
- CALL allocate_array( TColl_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( Sipj, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( TColl_i, ips_i,ipe_i, ijs_i,ije_i, ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( Sipj, ips_i,ipe_i, ijs_i,ije_i, ikys,ikye, ikxs,ikxe, izs,ize)
CALL allocate_array( xnipj, ips_i,ipe_i, ijs_i,ije_i)
CALL allocate_array( xnipp2j, ips_i,ipe_i)
CALL allocate_array( xnipp1j, ips_i,ipe_i)
@@ -105,13 +105,13 @@ SUBROUTINE memory
!! Collision matrices
IF (gyrokin_CO) THEN !GK collision matrices (one for each kperp)
IF (KIN_E) THEN
- CALL allocate_array( Ceepj, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( CeipjT, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( CeipjF, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxi+1)*(jmaxi+1), ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( CiepjT, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( CiepjF, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxe+1)*(jmaxe+1), ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( Ceepj, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( CeipjT, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( CeipjF, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxi+1)*(jmaxi+1), ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( CiepjT, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), ikys,ikye, ikxs,ikxe, izs,ize)
+ CALL allocate_array( CiepjF, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxe+1)*(jmaxe+1), ikys,ikye, ikxs,ikxe, izs,ize)
ENDIF
- CALL allocate_array( Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( Ciipj, 1,(pmaxi+1)*(jmaxi+1), 1,(pmaxi+1)*(jmaxi+1), ikys,ikye, ikxs,ikxe, izs,ize)
ELSE !DK collision matrix (same for every k)
IF (KIN_E) THEN
CALL allocate_array( Ceepj, 1,(pmaxe+1)*(jmaxe+1), 1,(pmaxe+1)*(jmaxe+1), 1,1, 1,1, 1,1)
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 27012952444ee2b504fea4dddaef45a38222e89d..9d4dc7a12113430be824b38167a06342f9819f8f 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -35,20 +35,20 @@ SUBROUTINE moments_eq_rhs_e
COMPLEX(dp) :: Tnepp1j, Tnepm1j, Tnepp1jm1, Tnepm1jm1 ! Terms from mirror force with non adiab moments
COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi
COMPLEX(dp) :: Unepm1j, Unepm1jp1, Unepm1jm1 ! Terms from mirror force with adiab moments
- COMPLEX(dp) :: i_ky,phikxkyz
+ COMPLEX(dp) :: i_ky,phikykxz
! Measuring execution time
CALL cpu_time(t0_rhs)
! Spatial loops
zloope : DO iz = izs,ize
- kyloope : DO iky = ikys,ikye
- ky = kyarray(iky) ! toroidal wavevector
- i_ky = imagu * ky ! toroidal derivative
+ kxloope : DO ikx = ikxs,ikxe
+ kx = kxarray(ikx) ! radial wavevector
- kxloope : DO ikx = ikxs,ikxe
- kx = kxarray(ikx) ! radial wavevector
- phikxkyz = phi(ikx,iky,iz)! tmp phi value
+ kyloope : DO iky = ikys,ikye
+ ky = kyarray(iky) ! toroidal wavevector
+ i_ky = imagu * ky ! toroidal derivative
+ phikykxz = phi(iky,ikx,iz)! tmp phi value
! Kinetic loops
jloope : DO ij = ijs_e, ije_e ! This loop is from 1 to jmaxi+1
@@ -57,50 +57,50 @@ SUBROUTINE moments_eq_rhs_e
ploope : DO ip = ips_e, ipe_e ! Hermite loop
p_int = parray_e(ip) ! Hermite degree
eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid
- kperp2= kparray(ikx,iky,iz,eo)**2
+ kperp2= kparray(iky,ikx,iz,eo)**2
!! Compute moments mixing terms
Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp
! Perpendicular dynamic
! term propto n_e^{p,j}
- Tnepj = xnepj(ip,ij)* nadiab_moments_e(ip,ij,ikx,iky,iz)
+ Tnepj = xnepj(ip,ij)* nadiab_moments_e(ip,ij,iky,ikx,iz)
! term propto n_e^{p+2,j}
- Tnepp2j = xnepp2j(ip) * nadiab_moments_e(ip+pp2,ij,ikx,iky,iz)
+ Tnepp2j = xnepp2j(ip) * nadiab_moments_e(ip+pp2,ij,iky,ikx,iz)
! term propto n_e^{p-2,j}
- Tnepm2j = xnepm2j(ip) * nadiab_moments_e(ip-pp2,ij,ikx,iky,iz)
+ Tnepm2j = xnepm2j(ip) * nadiab_moments_e(ip-pp2,ij,iky,ikx,iz)
! term propto n_e^{p,j+1}
- Tnepjp1 = xnepjp1(ij) * nadiab_moments_e(ip,ij+1,ikx,iky,iz)
+ Tnepjp1 = xnepjp1(ij) * nadiab_moments_e(ip,ij+1,iky,ikx,iz)
! term propto n_e^{p,j-1}
- Tnepjm1 = xnepjm1(ij) * nadiab_moments_e(ip,ij-1,ikx,iky,iz)
+ Tnepjm1 = xnepjm1(ij) * nadiab_moments_e(ip,ij-1,iky,ikx,iz)
! Parallel dynamic
! ddz derivative for Landau damping term
- Tpar = xnepp1j(ip) * ddz_nepj(ip+1,ij,ikx,iky,iz) &
- + xnepm1j(ip) * ddz_nepj(ip-1,ij,ikx,iky,iz)
+ Tpar = xnepp1j(ip) * ddz_nepj(ip+1,ij,iky,ikx,iz) &
+ + xnepm1j(ip) * ddz_nepj(ip-1,ij,iky,ikx,iz)
! Mirror terms
- Tnepp1j = ynepp1j (ip,ij) * interp_nepj(ip+1,ij ,ikx,iky,iz)
- Tnepp1jm1 = ynepp1jm1(ip,ij) * interp_nepj(ip+1,ij-1,ikx,iky,iz)
- Tnepm1j = ynepm1j (ip,ij) * interp_nepj(ip-1,ij ,ikx,iky,iz)
- Tnepm1jm1 = ynepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,ikx,iky,iz)
+ Tnepp1j = ynepp1j (ip,ij) * interp_nepj(ip+1,ij ,iky,ikx,iz)
+ Tnepp1jm1 = ynepp1jm1(ip,ij) * interp_nepj(ip+1,ij-1,iky,ikx,iz)
+ Tnepm1j = ynepm1j (ip,ij) * interp_nepj(ip-1,ij ,iky,ikx,iz)
+ Tnepm1jm1 = ynepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,iky,ikx,iz)
! Trapping terms
- Unepm1j = znepm1j (ip,ij) * interp_nepj(ip-1,ij ,ikx,iky,iz)
- Unepm1jp1 = znepm1jp1(ip,ij) * interp_nepj(ip-1,ij+1,ikx,iky,iz)
- Unepm1jm1 = znepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,ikx,iky,iz)
+ Unepm1j = znepm1j (ip,ij) * interp_nepj(ip-1,ij ,iky,ikx,iz)
+ Unepm1jp1 = znepm1jp1(ip,ij) * interp_nepj(ip-1,ij+1,iky,ikx,iz)
+ Unepm1jm1 = znepm1jm1(ip,ij) * interp_nepj(ip-1,ij-1,iky,ikx,iz)
Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + Unepm1j + Unepm1jp1 + Unepm1jm1
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
- Tphi = (xphij_i (ip,ij)*kernel_e(ij ,ikx,iky,iz,eo) &
- + xphijp1_i(ip,ij)*kernel_e(ij+1,ikx,iky,iz,eo) &
- + xphijm1_i(ip,ij)*kernel_e(ij-1,ikx,iky,iz,eo))*phikxkyz
+ Tphi = (xphij_i (ip,ij)*kernel_e(ij ,iky,ikx,iz,eo) &
+ + xphijp1_i(ip,ij)*kernel_e(ij+1,iky,ikx,iz,eo) &
+ + xphijm1_i(ip,ij)*kernel_e(ij-1,iky,ikx,iz,eo))*phikykxz
ELSE
Tphi = 0._dp
ENDIF
!! Sum of all RHS terms
- moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) = &
+ moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
- - imagu*Ckxky(ikx,iky,iz,eo)*hatR(iz,eo)* (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)&
+ - imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo)* (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)&
! Parallel coupling (Landau Damping)
- Tpar*gradz_coeff(iz,eo) &
! Mirror term (parallel magnetic gradient)
@@ -108,18 +108,18 @@ SUBROUTINE moments_eq_rhs_e
! Drives (density + temperature gradients)
- i_ky * Tphi &
! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose)
- - (mu_x*kx**4 + mu_y*ky**4)*moments_e(ip,ij,ikx,iky,iz,updatetlevel) &
+ - (mu_x*kx**4 + mu_y*ky**4)*moments_e(ip,ij,iky,ikx,iz,updatetlevel) &
! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)"
- + mu_z * diff_dz_coeff * ddz2_Nepj(ip,ij,ikx,iky,iz) &
+ + mu_z * diff_dz_coeff * ddz2_Nepj(ip,ij,iky,ikx,iz) &
! Collision term
- + TColl_e(ip,ij,ikx,iky,iz) &
+ + TColl_e(ip,ij,iky,ikx,iz) &
! Nonlinear term
- - Sepj(ip,ij,ikx,iky,iz)
+ - Sepj(ip,ij,iky,ikx,iz)
END DO ploope
END DO jloope
- END DO kxloope
- END DO kyloope
+ END DO kyloope
+ END DO kxloope
END DO zloope
! Execution time end
@@ -153,7 +153,7 @@ SUBROUTINE moments_eq_rhs_i
COMPLEX(dp) :: Tnipp1j, Tnipm1j, Tnipp1jm1, Tnipm1jm1 ! Terms from mirror force with non adiab moments
COMPLEX(dp) :: Unipm1j, Unipm1jp1, Unipm1jm1 ! Terms from mirror force with adiab moments
COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi
- COMPLEX(dp) :: i_ky, phikxkyz
+ COMPLEX(dp) :: i_ky, phikykxz
! Measuring execution time
CALL cpu_time(t0_rhs)
@@ -161,13 +161,13 @@ SUBROUTINE moments_eq_rhs_i
! Spatial loops
zloopi : DO iz = izs,ize
- kyloopi : DO iky = ikys,ikye
- ky = kyarray(iky) ! toroidal wavevector
- i_ky = imagu * ky ! toroidal derivative
+ kxloopi : DO ikx = ikxs,ikxe
+ kx = kxarray(ikx) ! radial wavevector
- kxloopi : DO ikx = ikxs,ikxe
- kx = kxarray(ikx) ! radial wavevector
- phikxkyz = phi(ikx,iky,iz)! tmp phi value
+ kyloopi : DO iky = ikys,ikye
+ ky = kyarray(iky) ! toroidal wavevector
+ i_ky = imagu * ky ! toroidal derivative
+ phikykxz = phi(iky,ikx,iz)! tmp phi value
! Kinetic loops
jloopi : DO ij = ijs_i, ije_i ! This loop is from 1 to jmaxi+1
@@ -176,52 +176,52 @@ SUBROUTINE moments_eq_rhs_i
ploopi : DO ip = ips_i, ipe_i ! Hermite loop
p_int = parray_i(ip) ! Hermite degree
eo = MODULO(p_int,2) ! Indicates if we are on odd or even z grid
- kperp2= kparray(ikx,iky,iz,eo)**2
+ kperp2= kparray(iky,ikx,iz,eo)**2
!! Compute moments mixing terms
Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp
! Perpendicular dynamic
! term propto n_i^{p,j}
- Tnipj = xnipj(ip,ij) * nadiab_moments_i(ip ,ij ,ikx,iky,iz)
+ Tnipj = xnipj(ip,ij) * nadiab_moments_i(ip ,ij ,iky,ikx,iz)
! term propto n_i^{p+2,j}
- Tnipp2j = xnipp2j(ip) * nadiab_moments_i(ip+pp2,ij ,ikx,iky,iz)
+ Tnipp2j = xnipp2j(ip) * nadiab_moments_i(ip+pp2,ij ,iky,ikx,iz)
! term propto n_i^{p-2,j}
- Tnipm2j = xnipm2j(ip) * nadiab_moments_i(ip-pp2,ij ,ikx,iky,iz)
+ Tnipm2j = xnipm2j(ip) * nadiab_moments_i(ip-pp2,ij ,iky,ikx,iz)
! term propto n_i^{p,j+1}
- Tnipjp1 = xnipjp1(ij) * nadiab_moments_i(ip ,ij+1,ikx,iky,iz)
+ Tnipjp1 = xnipjp1(ij) * nadiab_moments_i(ip ,ij+1,iky,ikx,iz)
! term propto n_i^{p,j-1}
- Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip ,ij-1,ikx,iky,iz)
+ Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip ,ij-1,iky,ikx,iz)
! Tperp
Tperp = Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1
! Parallel dynamic
! ddz derivative for Landau damping term
- Tpar = xnipp1j(ip) * ddz_nipj(ip+1,ij,ikx,iky,iz) &
- + xnipm1j(ip) * ddz_nipj(ip-1,ij,ikx,iky,iz)
+ Tpar = xnipp1j(ip) * ddz_nipj(ip+1,ij,iky,ikx,iz) &
+ + xnipm1j(ip) * ddz_nipj(ip-1,ij,iky,ikx,iz)
! Mirror terms
- Tnipp1j = ynipp1j (ip,ij) * interp_nipj(ip+1,ij ,ikx,iky,iz)
- Tnipp1jm1 = ynipp1jm1(ip,ij) * interp_nipj(ip+1,ij-1,ikx,iky,iz)
- Tnipm1j = ynipm1j (ip,ij) * interp_nipj(ip-1,ij ,ikx,iky,iz)
- Tnipm1jm1 = ynipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,ikx,iky,iz)
+ Tnipp1j = ynipp1j (ip,ij) * interp_nipj(ip+1,ij ,iky,ikx,iz)
+ Tnipp1jm1 = ynipp1jm1(ip,ij) * interp_nipj(ip+1,ij-1,iky,ikx,iz)
+ Tnipm1j = ynipm1j (ip,ij) * interp_nipj(ip-1,ij ,iky,ikx,iz)
+ Tnipm1jm1 = ynipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,iky,ikx,iz)
! Trapping terms
- Unipm1j = znipm1j (ip,ij) * interp_nipj(ip-1,ij ,ikx,iky,iz)
- Unipm1jp1 = znipm1jp1(ip,ij) * interp_nipj(ip-1,ij+1,ikx,iky,iz)
- Unipm1jm1 = znipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,ikx,iky,iz)
+ Unipm1j = znipm1j (ip,ij) * interp_nipj(ip-1,ij ,iky,ikx,iz)
+ Unipm1jp1 = znipm1jp1(ip,ij) * interp_nipj(ip-1,ij+1,iky,ikx,iz)
+ Unipm1jm1 = znipm1jm1(ip,ij) * interp_nipj(ip-1,ij-1,iky,ikx,iz)
Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + Unipm1j + Unipm1jp1 + Unipm1jm1
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
- Tphi = (xphij_i (ip,ij)*kernel_i(ij ,ikx,iky,iz,eo) &
- + xphijp1_i(ip,ij)*kernel_i(ij+1,ikx,iky,iz,eo) &
- + xphijm1_i(ip,ij)*kernel_i(ij-1,ikx,iky,iz,eo))*phikxkyz
+ Tphi = (xphij_i (ip,ij)*kernel_i(ij ,iky,ikx,iz,eo) &
+ + xphijp1_i(ip,ij)*kernel_i(ij+1,iky,ikx,iz,eo) &
+ + xphijm1_i(ip,ij)*kernel_i(ij-1,iky,ikx,iz,eo))*phikykxz
ELSE
Tphi = 0._dp
ENDIF
!! Sum of all RHS terms
- moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) = &
+ moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
- - imagu*Ckxky(ikx,iky,iz,eo)*hatR(iz,eo) * Tperp &
+ - imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp &
! Parallel coupling (Landau damping)
- gradz_coeff(iz,eo) * Tpar &
! Mirror term (parallel magnetic gradient)
@@ -229,18 +229,18 @@ SUBROUTINE moments_eq_rhs_i
! Drives (density + temperature gradients)
- i_ky * Tphi &
! Numerical hyperdiffusion (totally artificial, for stability purpose)
- - (mu_x*kx**4 + mu_y*ky**4)*moments_i(ip,ij,ikx,iky,iz,updatetlevel) &
+ - (mu_x*kx**4 + mu_y*ky**4)*moments_i(ip,ij,iky,ikx,iz,updatetlevel) &
! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)"
- + mu_z * diff_dz_coeff * ddz2_Nipj(ip,ij,ikx,iky,iz) &
+ + mu_z * diff_dz_coeff * ddz2_Nipj(ip,ij,iky,ikx,iz) &
! Collision term
- + TColl_i(ip,ij,ikx,iky,iz)&
+ + TColl_i(ip,ij,iky,ikx,iz)&
! Nonlinear term
- - Sipj(ip,ij,ikx,iky,iz)
+ - Sipj(ip,ij,iky,ikx,iz)
END DO ploopi
END DO jloopi
- END DO kxloopi
- END DO kyloopi
+ END DO kyloopi
+ END DO kxloopi
END DO zloopi
! Execution time end
diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90
index c70de47a619e10428f2fe0740c66df9a3415eea5..78ffce96817fe8a68febcd436eb73455525ac3c0 100644
--- a/src/nonlinear_mod.F90
+++ b/src/nonlinear_mod.F90
@@ -24,15 +24,15 @@ MODULE nonlinear
CONTAINS
SUBROUTINE nonlinear_init
- ALLOCATE( F_cmpx(ikxs:ikxe,ikys:ikye))
- ALLOCATE( G_cmpx(ikxs:ikxe,ikys:ikye))
+ ALLOCATE( F_cmpx(ikys:ikye,ikxs:ikxe))
+ ALLOCATE( G_cmpx(ikys:ikye,ikxs:ikxe))
- ALLOCATE(Fx_cmpx(ikxs:ikxe,ikys:ikye))
- ALLOCATE(Gy_cmpx(ikxs:ikxe,ikys:ikye))
- ALLOCATE(Fy_cmpx(ikxs:ikxe,ikys:ikye))
- ALLOCATE(Gx_cmpx(ikxs:ikxe,ikys:ikye))
+ ALLOCATE(Fx_cmpx(ikys:ikye,ikxs:ikxe))
+ ALLOCATE(Gy_cmpx(ikys:ikye,ikxs:ikxe))
+ ALLOCATE(Fy_cmpx(ikys:ikye,ikxs:ikxe))
+ ALLOCATE(Gx_cmpx(ikys:ikye,ikxs:ikxe))
- ALLOCATE(F_conv_G(ikxs:ikxe,ikys:ikye))
+ ALLOCATE(F_conv_G(ikys:ikye,ikxs:ikxe))
END SUBROUTINE nonlinear_init
SUBROUTINE compute_Sapj
@@ -101,7 +101,7 @@ SUBROUTINE compute_nonlinear
! Retrieve convolution in input format
DO ikx = ikxs, ikxe
DO iky = ikys, ikye
- Sepj(ip,ij,ikx,iky,iz) = bracket_sum_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ Sepj(ip,ij,iky,ikx,iz) = bracket_sum_c(ikx,iky-local_nky_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
ENDDO
ENDDO
ENDDO jloope
@@ -143,7 +143,7 @@ zloopi: DO iz = izs,ize
! Retrieve convolution in input format
DO ikx = ikxs, ikxe
DO iky = ikys, ikye
- Sipj(ip,ij,ikx,iky,iz) = bracket_sum_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky)
+ Sipj(ip,ij,iky,ikx,iz) = bracket_sum_c(ikx,iky-local_nky_offset)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
ENDDO jloopi
@@ -184,26 +184,26 @@ SUBROUTINE compute_semi_linear_ZF
ky = kyarray(iky)
kerneln = kernel_e(in, ikx, iky, iz, eo)
! Zonal terms (=0 for all ky not 0)
- Fx_cmpx(ikx,iky) = 0._dp
- Gx_cmpx(ikx,iky) = 0._dp
+ Fx_cmpx(iky,ikx) = 0._dp
+ Gx_cmpx(iky,ikx) = 0._dp
IF(iky .EQ. iky_0) THEN
- Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,iz) * kerneln
+ Fx_cmpx(iky,ikx) = imagu*kx* phi(iky,ikx,iz) * kerneln
smax = MIN( (in-1)+(ij-1), jmaxe );
DO is = 1, smax+1 ! sum truncation on number of moments
- Gx_cmpx(ikx,iky) = Gx_cmpx(ikx,iky) + &
- dnjs(in,ij,is) * moments_e(ip,is,ikx,iky,iz,updatetlevel)
+ Gx_cmpx(iky,ikx) = Gx_cmpx(iky,ikx) + &
+ dnjs(in,ij,is) * moments_e(ip,is,iky,ikx,iz,updatetlevel)
ENDDO
- Gx_cmpx(ikx,iky) = imagu*kx*Gx_cmpx(ikx,iky)
+ Gx_cmpx(iky,ikx) = imagu*kx*Gx_cmpx(iky,ikx)
ENDIF
! NZ terms
- Fy_cmpx(ikx,iky) = imagu*ky* phi(ikx,iky,iz) * kerneln
- Gy_cmpx(ikx,iky) = 0._dp ! initialization of the sum
+ Fy_cmpx(iky,ikx) = imagu*ky* phi(iky,ikx,iz) * kerneln
+ Gy_cmpx(iky,ikx) = 0._dp ! initialization of the sum
smax = MIN( (in-1)+(ij-1), jmaxe );
DO is = 1, smax+1 ! sum truncation on number of moments
- Gy_cmpx(ikx,iky) = Gy_cmpx(ikx,iky) + &
- dnjs(in,ij,is) * moments_e(ip,is,ikx,iky,iz,updatetlevel)
+ Gy_cmpx(iky,ikx) = Gy_cmpx(iky,ikx) + &
+ dnjs(in,ij,is) * moments_e(ip,is,iky,ikx,iz,updatetlevel)
ENDDO
- Gy_cmpx(ikx,iky) = imagu*ky*Gy_cmpx(ikx,iky)
+ Gy_cmpx(iky,ikx) = imagu*ky*Gy_cmpx(iky,ikx)
ENDDO kyloope
ENDDO kxloope
! First term df/dx x dg/dy
@@ -216,7 +216,7 @@ SUBROUTINE compute_semi_linear_ZF
! Retrieve convolution in input format
DO ikx = ikxs, ikxe
DO iky = ikys, ikye
- Sepj(ip,ij,ikx,iky,iz) = bracket_sum_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ Sepj(ip,ij,iky,ikx,iz) = bracket_sum_c(ikx,iky-local_nky_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
ENDDO
ENDDO
ENDDO jloope
@@ -244,27 +244,27 @@ zloopi: DO iz = izs,ize
kxloopi: DO ikx = ikxs,ikxe ! Loop over kx
kyloopi: DO iky = ikys,ikye ! Loop over ky
! Zonal terms (=0 for all ky not 0)
- Fx_cmpx(ikx,iky) = 0._dp
- Gx_cmpx(ikx,iky) = 0._dp
+ Fx_cmpx(iky,ikx) = 0._dp
+ Gx_cmpx(iky,ikx) = 0._dp
IF(iky .EQ. iky_0) THEN
- Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,iz) * kerneln
+ Fx_cmpx(iky,ikx) = imagu*kx* phi(iky,ikx,iz) * kerneln
smax = MIN( (in-1)+(ij-1), jmaxi );
DO is = 1, smax+1 ! sum truncation on number of moments
- Gx_cmpx(ikx,iky) = Gx_cmpx(ikx,iky) + &
- dnjs(in,ij,is) * moments_i(ip,is,ikx,iky,iz,updatetlevel)
+ Gx_cmpx(iky,ikx) = Gx_cmpx(iky,ikx) + &
+ dnjs(in,ij,is) * moments_i(ip,is,iky,ikx,iz,updatetlevel)
ENDDO
- Gx_cmpx(ikx,iky) = imagu*kx*Gx_cmpx(ikx,iky)
+ Gx_cmpx(iky,ikx) = imagu*kx*Gx_cmpx(iky,ikx)
ENDIF
! NZ terms
- Fy_cmpx(ikx,iky) = imagu*ky* phi(ikx,iky,iz) * kerneln
- Gy_cmpx(ikx,iky) = 0._dp ! initialization of the sum
+ Fy_cmpx(iky,ikx) = imagu*ky* phi(iky,ikx,iz) * kerneln
+ Gy_cmpx(iky,ikx) = 0._dp ! initialization of the sum
smax = MIN( (in-1)+(ij-1), jmaxi );
DO is = 1, smax+1 ! sum truncation on number of moments
- Gy_cmpx(ikx,iky) = Gy_cmpx(ikx,iky) + &
- dnjs(in,ij,is) * moments_i(ip,is,ikx,iky,iz,updatetlevel)
+ Gy_cmpx(iky,ikx) = Gy_cmpx(iky,ikx) + &
+ dnjs(in,ij,is) * moments_i(ip,is,iky,ikx,iz,updatetlevel)
ENDDO
- Gy_cmpx(ikx,iky) = imagu*ky*Gy_cmpx(ikx,iky)
+ Gy_cmpx(iky,ikx) = imagu*ky*Gy_cmpx(iky,ikx)
ENDDO kyloopi
ENDDO kxloopi
! First term drphi x dzf
@@ -277,7 +277,7 @@ zloopi: DO iz = izs,ize
! Retrieve convolution in input format
DO ikx = ikxs, ikxe
DO iky = ikys, ikye
- Sipj(ip,ij,ikx,iky,iz) = bracket_sum_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky)
+ Sipj(ip,ij,iky,ikx,iz) = bracket_sum_c(ikx,iky-local_nky_offset)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
ENDDO jloopi
@@ -316,25 +316,25 @@ SUBROUTINE compute_semi_linear_NZ
ky = kyarray(iky)
kerneln = kernel_e(in, ikx, iky, iz, eo)
! All terms
- Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,iz) * kerneln
+ Fx_cmpx(iky,ikx) = imagu*kx* phi(iky,ikx,iz) * kerneln
smax = MIN( (in-1)+(ij-1), jmaxe );
DO is = 1, smax+1 ! sum truncation on number of moments
- Gx_cmpx(ikx,iky) = Gx_cmpx(ikx,iky) + &
- dnjs(in,ij,is) * moments_e(ip,is,ikx,iky,iz,updatetlevel)
+ Gx_cmpx(iky,ikx) = Gx_cmpx(iky,ikx) + &
+ dnjs(in,ij,is) * moments_e(ip,is,iky,ikx,iz,updatetlevel)
ENDDO
- Gx_cmpx(ikx,iky) = imagu*kx*Gx_cmpx(ikx,iky)
+ Gx_cmpx(iky,ikx) = imagu*kx*Gx_cmpx(iky,ikx)
! Kx = 0 terms
- Fy_cmpx(ikx,iky) = 0._dp
- Gy_cmpx(ikx,iky) = 0._dp
+ Fy_cmpx(iky,ikx) = 0._dp
+ Gy_cmpx(iky,ikx) = 0._dp
IF (ikx .EQ. ikx_0) THEN
- Fy_cmpx(ikx,iky) = imagu*ky* phi(ikx,iky,iz) * kerneln
- Gy_cmpx(ikx,iky) = 0._dp ! initialization of the sum
+ Fy_cmpx(iky,ikx) = imagu*ky* phi(iky,ikx,iz) * kerneln
+ Gy_cmpx(iky,ikx) = 0._dp ! initialization of the sum
smax = MIN( (in-1)+(ij-1), jmaxe );
DO is = 1, smax+1 ! sum truncation on number of moments
- Gy_cmpx(ikx,iky) = Gy_cmpx(ikx,iky) + &
- dnjs(in,ij,is) * moments_e(ip,is,ikx,iky,iz,updatetlevel)
+ Gy_cmpx(iky,ikx) = Gy_cmpx(iky,ikx) + &
+ dnjs(in,ij,is) * moments_e(ip,is,iky,ikx,iz,updatetlevel)
ENDDO
- Gy_cmpx(ikx,iky) = imagu*ky*Gy_cmpx(ikx,iky)
+ Gy_cmpx(iky,ikx) = imagu*ky*Gy_cmpx(iky,ikx)
ENDIF
ENDDO kyloope
ENDDO kxloope
@@ -348,7 +348,7 @@ SUBROUTINE compute_semi_linear_NZ
! Retrieve convolution in input format
DO ikx = ikxs, ikxe
DO iky = ikys, ikye
- Sepj(ip,ij,ikx,iky,iz) = bracket_sum_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ Sepj(ip,ij,iky,ikx,iz) = bracket_sum_c(ikx,iky-local_nky_offset)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
ENDDO
ENDDO
ENDDO jloope
@@ -376,26 +376,26 @@ zloopi: DO iz = izs,ize
kxloopi: DO ikx = ikxs,ikxe ! Loop over kx
kyloopi: DO iky = ikys,ikye ! Loop over ky
! Zonal terms (=0 for all ky not 0)
- Fx_cmpx(ikx,iky) = imagu*kx* phi(ikx,iky,iz) * kerneln
+ Fx_cmpx(iky,ikx) = imagu*kx* phi(iky,ikx,iz) * kerneln
smax = MIN( (in-1)+(ij-1), jmaxi );
DO is = 1, smax+1 ! sum truncation on number of moments
- Gx_cmpx(ikx,iky) = Gx_cmpx(ikx,iky) + &
- dnjs(in,ij,is) * moments_i(ip,is,ikx,iky,iz,updatetlevel)
+ Gx_cmpx(iky,ikx) = Gx_cmpx(iky,ikx) + &
+ dnjs(in,ij,is) * moments_i(ip,is,iky,ikx,iz,updatetlevel)
ENDDO
- Gx_cmpx(ikx,iky) = imagu*kx*Gx_cmpx(ikx,iky)
+ Gx_cmpx(iky,ikx) = imagu*kx*Gx_cmpx(iky,ikx)
! Kx = 0 terms
- Fy_cmpx(ikx,iky) = 0._dp
- Gy_cmpx(ikx,iky) = 0._dp
+ Fy_cmpx(iky,ikx) = 0._dp
+ Gy_cmpx(iky,ikx) = 0._dp
IF (ikx .EQ. ikx_0) THEN
- Fy_cmpx(ikx,iky) = imagu*ky* phi(ikx,iky,iz) * kerneln
- Gy_cmpx(ikx,iky) = 0._dp ! initialization of the sum
+ Fy_cmpx(iky,ikx) = imagu*ky* phi(iky,ikx,iz) * kerneln
+ Gy_cmpx(iky,ikx) = 0._dp ! initialization of the sum
smax = MIN( (in-1)+(ij-1), jmaxi );
DO is = 1, smax+1 ! sum truncation on number of moments
- Gy_cmpx(ikx,iky) = Gy_cmpx(ikx,iky) + &
- dnjs(in,ij,is) * moments_i(ip,is,ikx,iky,iz,updatetlevel)
+ Gy_cmpx(iky,ikx) = Gy_cmpx(iky,ikx) + &
+ dnjs(in,ij,is) * moments_i(ip,is,iky,ikx,iz,updatetlevel)
ENDDO
- Gy_cmpx(ikx,iky) = imagu*ky*Gy_cmpx(ikx,iky)
+ Gy_cmpx(iky,ikx) = imagu*ky*Gy_cmpx(iky,ikx)
ENDIF
ENDDO kyloopi
ENDDO kxloopi
@@ -409,7 +409,7 @@ zloopi: DO iz = izs,ize
! Retrieve convolution in input format
DO ikx = ikxs, ikxe
DO iky = ikys, ikye
- Sipj(ip,ij,ikx,iky,iz) = bracket_sum_c(iky,ikx-local_nkx_offset)*AA_x(ikx)*AA_y(iky)
+ Sipj(ip,ij,iky,ikx,iz) = bracket_sum_c(ikx,iky-local_nky_offset)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
ENDDO jloopi
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index ba8d0653e120c44d9ba44c48e120751cf89b3616..77d0d35749e853075968b23cdf9ac5d5f74ff357 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -55,7 +55,7 @@ SUBROUTINE evaluate_kernels
lambdaD, CLOS, sigmae2_taue_o2, sigmai2_taui_o2, KIN_E
IMPLICIT NONE
INTEGER :: j_int
- REAL(dp) :: j_dp, y_, kp2_, kx_, ky_
+ REAL(dp) :: j_dp, y_
DO eo = 0,1
DO ikx = ikxs,ikxe
@@ -66,21 +66,21 @@ DO iz = izgs,izge
DO ij = ijgs_e, ijge_e
j_int = jarray_e(ij)
j_dp = REAL(j_int,dp)
- y_ = sigmae2_taue_o2 * kparray(ikx,iky,iz,eo)**2
- kernel_e(ij,ikx,iky,iz,eo) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
+ y_ = sigmae2_taue_o2 * kparray(iky,ikx,iz,eo)**2
+ kernel_e(ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
ENDDO
IF (ijs_e .EQ. 1) &
- kernel_e(ijgs_e,ikx,iky,iz,eo) = 0._dp
+ kernel_e(ijgs_e,iky,ikx,iz,eo) = 0._dp
ENDIF
!!!!! Ion kernels !!!!!
DO ij = ijgs_i, ijge_i
j_int = jarray_i(ij)
j_dp = REAL(j_int,dp)
- y_ = sigmai2_taui_o2 * kparray(ikx,iky,iz,eo)**2
- kernel_i(ij,ikx,iky,iz,eo) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
+ y_ = sigmai2_taui_o2 * kparray(iky,ikx,iz,eo)**2
+ kernel_i(ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
ENDDO
IF (ijs_i .EQ. 1) &
- kernel_i(ijgs_i,ikx,iky,iz,eo) = 0._dp
+ kernel_i(ijgs_i,iky,ikx,iz,eo) = 0._dp
ENDDO
ENDDO
ENDDO
@@ -107,13 +107,13 @@ SUBROUTINE evaluate_poisson_op
kyloop: DO iky = ikys,ikye
zloop: DO iz = izs,ize
IF( (kxarray(ikx).EQ.0._dp) .AND. (kyarray(iky).EQ.0._dp) ) THEN
- inv_poisson_op(ikx, iky, iz) = 0._dp
+ inv_poisson_op(iky, ikx, iz) = 0._dp
ELSE
!!!!!!!!!!!!!!!!! Ion contribution
! loop over n only if the max polynomial degree
pol_i = 0._dp
DO ini=1,jmaxi+1
- pol_i = pol_i + qi2_taui*kernel_i(ini,ikx,iky,iz,0)**2 ! ... sum recursively ...
+ pol_i = pol_i + qi2_taui*kernel_i(ini,iky,ikx,iz,0)**2 ! ... sum recursively ...
END DO
!!!!!!!!!!!!! Electron contribution
pol_e = 0._dp
@@ -121,13 +121,13 @@ SUBROUTINE evaluate_poisson_op
IF (KIN_E) THEN
! loop over n only if the max polynomial degree
DO ine=1,jmaxe+1 ! ine = n+1
- pol_e = pol_e + qe2_taue*kernel_e(ine,ikx,iky,iz,0)**2 ! ... sum recursively ...
+ pol_e = pol_e + qe2_taue*kernel_e(ine,iky,ikx,iz,0)**2 ! ... sum recursively ...
END DO
! Adiabatic model
ELSE
pol_e = 1._dp - qe2_taue
ENDIF
- inv_poisson_op(ikx, iky, iz) = 1._dp/(qe2_taue + qi2_taui - pol_i - pol_e)
+ inv_poisson_op(iky, ikx, iz) = 1._dp/(qe2_taue + qi2_taui - pol_i - pol_e)
ENDIF
END DO zloop
END DO kyloop
@@ -281,15 +281,22 @@ SUBROUTINE save_EM_ZF_modes
USE model, ONLY: KIN_E
IMPLICIT NONE
! Store Zonal and entropy modes
+ IF(contains_ky0) THEN
IF(KIN_E) &
- moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize) = moments_e(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,iky_0,izs:ize,updatetlevel)
- moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize) = moments_i(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,iky_0,izs:ize,updatetlevel)
- phi_ZF(ikxs:ikxe,izs:ize) = phi(ikxs:ikxe,iky_0,izs:ize)
+ moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize) = moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel)
+ moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize) = moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel)
+ phi_ZF(ikxs:ikxe,izs:ize) = phi(iky_0,ikxs:ikxe,izs:ize)
+ ELSE
+ IF(KIN_E) &
+ moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize) = 0._dp
+ moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize) = 0._dp
+ phi_ZF(ikxs:ikxe,izs:ize) = 0._dp
+ ENDIF
IF(contains_kx0) THEN
IF(KIN_E) &
- moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize) = moments_e(ips_e:ipe_e,ijs_e:ije_e,ikx_0,ikys:ikye,izs:ize,updatetlevel)
- moments_i_EM(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,izs:ize) = moments_i(ips_i:ipe_i,ijs_i:ije_i,ikx_0,ikys:ikye,izs:ize,updatetlevel)
- phi_EM(ikys:ikye,izs:ize) = phi(ikx_0,ikys:ikye,izs:ize)
+ moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize) = moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikx_0,izs:ize,updatetlevel)
+ moments_i_EM(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,izs:ize) = moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikx_0,izs:ize,updatetlevel)
+ phi_EM(ikys:ikye,izs:ize) = phi(ikys:ikye,ikx_0,izs:ize)
ELSE
IF(KIN_E) &
moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize) = 0._dp
@@ -311,22 +318,22 @@ SUBROUTINE play_with_modes
SELECT CASE(ACT_ON_MODES)
CASE('wipe_zonal') ! Errase the zonal flow
IF(KIN_E) &
- moments_e(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,iky_0,izs:ize,updatetlevel) = 0._dp
- moments_i(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,iky_0,izs:ize,updatetlevel) = 0._dp
- phi(ikxs:ikxe,iky_0,izs:ize) = 0._dp
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = 0._dp
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = 0._dp
+ phi(iky_0,ikxs:ikxe,izs:ize) = 0._dp
CASE('wipe_entropymode')
IF(KIN_E) &
- moments_e(ips_e:ipe_e,ijs_e:ije_e,ikx_0,ikys:ikye,izs:ize,updatetlevel) = 0._dp
- moments_i(ips_i:ipe_i,ijs_i:ije_i,ikx_0,ikys:ikye,izs:ize,updatetlevel) = 0._dp
- phi(ikx_0,ikys:ikye,izs:ize) = 0._dp
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikx_0,izs:ize,updatetlevel) = 0._dp
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikx_0,izs:ize,updatetlevel) = 0._dp
+ phi(ikys:ikye,ikx_0,izs:ize) = 0._dp
CASE('wipe_turbulence')
DO ikx = ikxs,ikxe
DO iky = ikys, ikye
IF ( (ikx .NE. ikx_0) .AND. (iky .NE. iky_0) ) THEN
IF(KIN_E) &
- moments_e(ips_e:ipe_e,ijs_e:ije_e,ikx,iky,izs:ize,updatetlevel) = 0._dp
- moments_i(ips_i:ipe_i,ijs_i:ije_i,ikx,iky,izs:ize,updatetlevel) = 0._dp
- phi(ikx,iky,izs:ize) = 0._dp
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel) = 0._dp
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,iky,ikx,izs:ize,updatetlevel) = 0._dp
+ phi(iky,ikx,izs:ize) = 0._dp
ENDIF
ENDDO
ENDDO
@@ -335,29 +342,29 @@ SUBROUTINE play_with_modes
DO iky = ikys, ikye
IF ( (ikx .NE. ikx_0) ) THEN
IF(KIN_E) &
- moments_e(ips_e:ipe_e,ijs_e:ije_e,ikx,iky,izs:ize,updatetlevel) = 0._dp
- moments_i(ips_i:ipe_i,ijs_i:ije_i,ikx,iky,izs:ize,updatetlevel) = 0._dp
- phi(ikx,iky,izs:ize) = 0._dp
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,iky,ikx,izs:ize,updatetlevel) = 0._dp
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,iky,ikx,izs:ize,updatetlevel) = 0._dp
+ phi(iky,ikx,izs:ize) = 0._dp
ENDIF
ENDDO
ENDDO
CASE('freeze_zonal')
IF(KIN_E) &
- moments_e(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,iky_0,izs:ize,updatetlevel) = moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize)
- moments_i(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,iky_0,izs:ize,updatetlevel) = moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize)
- phi(ikxs:ikxe,iky_0,izs:ize) = phi_ZF(ikxs:ikxe,izs:ize)
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize)
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize)
+ phi(iky_0,ikxs:ikxe,izs:ize) = phi_ZF(ikxs:ikxe,izs:ize)
CASE('freeze_entropymode')
IF(contains_kx0) THEN
IF(KIN_E) &
- moments_e(ips_e:ipe_e,ijs_e:ije_e,ikx_0,ikys:ikye,izs:ize,updatetlevel) = moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize)
- moments_i(ips_i:ipe_i,ijs_i:ije_i,ikx_0,ikys:ikye,izs:ize,updatetlevel) = moments_i_EM(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,izs:ize)
- phi(ikx_0,ikys:ikye,izs:ize) = phi_EM(ikys:ikye,izs:ize)
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,ikx_0,izs:ize,updatetlevel) = moments_e_EM(ips_e:ipe_e,ijs_e:ije_e,ikys:ikye,izs:ize)
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,ikx_0,izs:ize,updatetlevel) = moments_i_EM(ips_i:ipe_i,ijs_i:ije_i,ikys:ikye,izs:ize)
+ phi(ikys:ikye,ikx_0,izs:ize) = phi_EM(ikys:ikye,izs:ize)
ENDIF
CASE('amplify_zonal')
IF(KIN_E) &
- moments_e(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,iky_0,izs:ize,updatetlevel) = AMP*moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize)
- moments_i(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,iky_0,izs:ize,updatetlevel) = AMP*moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize)
- phi(ikxs:ikxe,iky_0,izs:ize) = AMP*phi_ZF(ikxs:ikxe,izs:ize)
+ moments_e(ips_e:ipe_e,ijs_e:ije_e,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = AMP*moments_e_ZF(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,izs:ize)
+ moments_i(ips_i:ipe_i,ijs_i:ije_i,iky_0,ikxs:ikxe,izs:ize,updatetlevel) = AMP*moments_i_ZF(ips_i:ipe_i,ijs_i:ije_i,ikxs:ikxe,izs:ize)
+ phi(iky_0,ikxs:ikxe,izs:ize) = AMP*phi_ZF(ikxs:ikxe,izs:ize)
END SELECT
END SUBROUTINE
diff --git a/src/poisson.F90 b/src/poisson.F90
index c30d9d1298e016c8de0beabb5312df95bd06d44a..3936324eee0339353c29520b0486c7fd79ef04ee 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -17,7 +17,6 @@ SUBROUTINE poisson
INTEGER :: ini,ine, i_, root_bcast
COMPLEX(dp) :: fa_phi, intf_ ! current flux averaged phi
INTEGER :: count !! mpi integer to broadcast the electric potential at the end
- COMPLEX(dp) :: buffer(ikxs:ikxe,ikys:ikye)
COMPLEX(dp), DIMENSION(izs:ize) :: rho_i, rho_e ! charge density q_a n_a
! Execution time start
@@ -33,7 +32,7 @@ SUBROUTINE poisson
DO ini=1,jmaxi+1
DO iz = izs,ize
rho_i(iz) = rho_i(iz) &
- +q_i*kernel_i(ini,ikx,iky,iz,0)*moments_i(ip0_i,ini,ikx,iky,iz,updatetlevel)
+ +q_i*kernel_i(ini,iky,ikx,iz,0)*moments_i(ip0_i,ini,iky,ikx,iz,updatetlevel)
ENDDO
END DO
@@ -43,7 +42,7 @@ SUBROUTINE poisson
DO ine=1,jmaxe+1
DO iz = izs,ize
rho_e(iz) = rho_e(iz) &
- +q_e*kernel_e(ine,ikx,iky,iz,0)*moments_e(ip0_e,ine,ikx,iky,iz,updatetlevel)
+ +q_e*kernel_e(ine,iky,ikx,iz,0)*moments_e(ip0_e,ine,iky,ikx,iz,updatetlevel)
ENDDO
END DO
ELSE ! Adiabatic electrons
@@ -51,8 +50,8 @@ SUBROUTINE poisson
fa_phi = 0._dp
IF(kyarray(iky).EQ.0._dp) THEN ! take ky=0 mode (average)
DO ini=1,jmaxi+1 ! some over FLR contributions
- rho_e(izs:ize) = Jacobian(izs:ize,0)*moments_i(ip0_i,ini,ikx,iky,izs:ize,updatetlevel)&
- *kernel_i(ini,ikx,iky,izs:ize,0)*(inv_poisson_op(ikx,iky,izs:ize)-1._dp)
+ rho_e(izs:ize) = Jacobian(izs:ize,0)*moments_i(ip0_i,ini,iky,ikx,izs:ize,updatetlevel)&
+ *kernel_i(ini,iky,ikx,izs:ize,0)*(inv_poisson_op(iky,ikx,izs:ize)-1._dp)
call simpson_rule_z(rho_e(izs:ize),intf_) ! integrate over z
fa_phi = fa_phi + intf_
ENDDO
@@ -62,18 +61,18 @@ SUBROUTINE poisson
!!!!!!!!!!!!!!! Inverting the poisson equation !!!!!!!!!!!!!!!!!!!!!!!!!!
DO iz = izs,ize
- phi(ikx, iky, iz) = (rho_e(iz) + rho_i(iz))*inv_poisson_op(ikx,iky,iz)
+ phi(iky, ikx, iz) = (rho_e(iz) + rho_i(iz))*inv_poisson_op(iky,ikx,iz)
ENDDO
END DO kyloop
END DO kxloop
! Cancel origin singularity
- IF (contains_kx0 .AND. contains_ky0) phi(ikx_0,iky_0,:) = 0._dp
+ IF (contains_kx0 .AND. contains_ky0) phi(iky_0,ikx_0,:) = 0._dp
ENDIF
! Transfer phi to all the others process along p
- CALL manual_3D_bcast(phi(ikxs:ikxe,ikys:ikye,izs:ize))
+ CALL manual_3D_bcast(phi(ikys:ikye,ikxs:ikxe,izs:ize))
! Execution time end
CALL cpu_time(t1_poisson)
diff --git a/src/ppinit.F90 b/src/ppinit.F90
index b403fdbad15ad14d3bbd5c2bbcaffb1aa455e573..1a662babc6bc41ab9188b17ded43759e49d0bffc 100644
--- a/src/ppinit.F90
+++ b/src/ppinit.F90
@@ -37,13 +37,13 @@ SUBROUTINE ppinit
END IF
num_procs_p = dims(1) ! Number of processes along p
- num_procs_kx = dims(2) ! Number of processes along kx
+ num_procs_ky = dims(2) ! Number of processes along kx
num_procs_z = dims(3) ! Number of processes along z
!
!periodicity in p
periods(1)=.FALSE.
- !periodicity in kx
+ !periodicity in ky
periods(2)=.FALSE.
!periodicity in z
periods(3)=.TRUE.
@@ -58,11 +58,11 @@ SUBROUTINE ppinit
! Partitions 3-dim cartesian topology of comm0 into 1-dim cartesian subgrids
!
CALL MPI_CART_SUB (comm0, (/.TRUE.,.FALSE.,.FALSE./), comm_p, ierr)
- CALL MPI_CART_SUB (comm0, (/.FALSE.,.TRUE.,.FALSE./), comm_kx, ierr)
+ CALL MPI_CART_SUB (comm0, (/.FALSE.,.TRUE.,.FALSE./), comm_ky, ierr)
CALL MPI_CART_SUB (comm0, (/.FALSE.,.FALSE.,.TRUE./), comm_z, ierr)
! Find id inside the 1d-sub communicators
CALL MPI_COMM_RANK(comm_p, rank_p, ierr)
- CALL MPI_COMM_RANK(comm_kx, rank_kx, ierr)
+ CALL MPI_COMM_RANK(comm_ky, rank_ky, ierr)
CALL MPI_COMM_RANK(comm_z, rank_z, ierr)
! Find neighbours
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index abcd31dbd4e5d46aa9adfa81edf30df8f9b03c98..84ca47df51f8d639c4d174b3998076a006e3c003 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -38,7 +38,7 @@ SUBROUTINE compute_radial_ion_transport
ky_ = kyarray(iky)
DO ikx = ikxs,ikxe
DO iz = izgs,izge
- integrant(iz) = Jacobian(iz,0)*imagu * ky_ * moments_i(ip0_i,1,ikx,iky,iz,updatetlevel) * CONJG(phi(ikx,iky,iz))
+ integrant(iz) = Jacobian(iz,0)*imagu * ky_ * moments_i(ip0_i,1,iky,ikx,iz,updatetlevel) * CONJG(phi(iky,ikx,iz))
ENDDO
! Integrate over z
call simpson_rule_z(integrant,integral)
@@ -54,8 +54,8 @@ SUBROUTINE compute_radial_ion_transport
DO ij = ijs_i, ije_i
DO iz = izgs,izge
integrant(iz) = integrant(iz) + &
- Jacobian(iz,0)*imagu * ky_ * kernel_i(ij,ikx,iky,iz,0) &
- *moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel) * CONJG(phi(ikx,iky,iz))
+ Jacobian(iz,0)*imagu * ky_ * kernel_i(ij,iky,ikx,iz,0) &
+ *moments_i(ip0_i,ij,iky,ikx,iz,updatetlevel) * CONJG(phi(iky,ikx,iz))
ENDDO
ENDDO
! Integrate over z
@@ -71,15 +71,15 @@ SUBROUTINE compute_radial_ion_transport
!Gather manually among the rank_p=0 processes and perform the sum
gflux_ri = 0
pflux_ri = 0
- IF (num_procs_kx .GT. 1) THEN
+ IF (num_procs_ky .GT. 1) THEN
!! Everyone sends its local_sum to root = 0
- IF (rank_kx .NE. root) THEN
- CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_kx, ierr)
+ IF (rank_ky .NE. root) THEN
+ CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr)
ELSE
! Recieve from all the other processes
- DO i_ = 0,num_procs_kx-1
- IF (i_ .NE. rank_kx) &
- CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_kx, MPI_STATUS_IGNORE, ierr)
+ DO i_ = 0,num_procs_ky-1
+ IF (i_ .NE. rank_ky) &
+ CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
gflux_ri = gflux_ri + buffer(1)
pflux_ri = pflux_ri + buffer(2)
ENDDO
@@ -113,23 +113,23 @@ SUBROUTINE compute_radial_heatflux
integrant = 0._dp
DO ij = ijs_i, ije_i
DO iz = izgs,izge
- integrant(iz) = integrant(iz) + Jacobian(iz,0)*imagu*ky_*CONJG(phi(ikx,iky,iz))&
- *(twothird * ( 2._dp*j_dp * kernel_i(ij ,ikx,iky,iz,0) &
- - (j_dp+1) * kernel_i(ij+1,ikx,iky,iz,0) &
- - j_dp * kernel_i(ij-1,ikx,iky,iz,0))&
- * (moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel)+qi_taui*kernel_i(ij,ikx,iky,iz,0)*phi(ikx,iky,iz)) &
- + SQRT2*onethird * kernel_i(ij,ikx,iky,iz,0) * moments_i(ip2_i,ij,ikx,iky,iz,updatetlevel))
+ integrant(iz) = integrant(iz) + Jacobian(iz,0)*imagu*ky_*CONJG(phi(iky,ikx,iz))&
+ *(twothird * ( 2._dp*j_dp * kernel_i(ij ,iky,ikx,iz,0) &
+ - (j_dp+1) * kernel_i(ij+1,iky,ikx,iz,0) &
+ - j_dp * kernel_i(ij-1,iky,ikx,iz,0))&
+ * (moments_i(ip0_i,ij,iky,ikx,iz,updatetlevel)+qi_taui*kernel_i(ij,iky,ikx,iz,0)*phi(iky,ikx,iz)) &
+ + SQRT2*onethird * kernel_i(ij,iky,ikx,iz,0) * moments_i(ip2_i,ij,iky,ikx,iz,updatetlevel))
ENDDO
ENDDO
IF(KIN_E) THEN
DO ij = ijs_e, ije_e
DO iz = izgs,izge
- integrant(iz) = integrant(iz) + Jacobian(iz,0)*imagu*ky_*CONJG(phi(ikx,iky,iz))&
- *(twothird * ( 2._dp*j_dp * kernel_e(ij ,ikx,iky,iz,0) &
- - (j_dp+1) * kernel_e(ij+1,ikx,iky,iz,0) &
- - j_dp * kernel_e(ij-1,ikx,iky,iz,0))&
- * (moments_e(ip0_e,ij,ikx,iky,iz,updatetlevel)+qe_taue*kernel_e(ij,ikx,iky,iz,0)*phi(ikx,iky,iz)) &
- + SQRT2*onethird * kernel_e(ij,ikx,iky,iz,0) * moments_e(ip2_e,ij,ikx,iky,iz,updatetlevel))
+ integrant(iz) = integrant(iz) + Jacobian(iz,0)*imagu*ky_*CONJG(phi(iky,ikx,iz))&
+ *(twothird * ( 2._dp*j_dp * kernel_e(ij ,iky,ikx,iz,0) &
+ - (j_dp+1) * kernel_e(ij+1,iky,ikx,iz,0) &
+ - j_dp * kernel_e(ij-1,iky,ikx,iz,0))&
+ * (moments_e(ip0_e,ij,iky,ikx,iz,updatetlevel)+qe_taue*kernel_e(ij,iky,ikx,iz,0)*phi(iky,ikx,iz)) &
+ + SQRT2*onethird * kernel_e(ij,iky,ikx,iz,0) * moments_e(ip2_e,ij,iky,ikx,iz,updatetlevel))
ENDDO
ENDDO
ENDIF
@@ -142,15 +142,15 @@ SUBROUTINE compute_radial_heatflux
root = 0
!Gather manually among the rank_p=0 processes and perform the sum
hflux_x = 0
- IF (num_procs_kx .GT. 1) THEN
+ IF (num_procs_ky .GT. 1) THEN
!! Everyone sends its local_sum to root = 0
- IF (rank_kx .NE. root) THEN
- CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_kx, ierr)
+ IF (rank_ky .NE. root) THEN
+ CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_ky, ierr)
ELSE
! Recieve from all the other processes
- DO i_ = 0,num_procs_kx-1
- IF (i_ .NE. rank_kx) &
- CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_kx, MPI_STATUS_IGNORE, ierr)
+ DO i_ = 0,num_procs_ky-1
+ IF (i_ .NE. rank_ky) &
+ CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_ky, MPI_STATUS_IGNORE, ierr)
hflux_x = hflux_x + buffer(2)
ENDDO
ENDIF
@@ -209,36 +209,36 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
!------------- INTERP AND GRADIENTS ALONG Z ----------------------------------
IF (KIN_E) THEN
- DO iky = ikys,ikye
- DO ikx = ikxs,ikxe
+ DO ikx = ikxs,ikxe
+ DO iky = ikys,ikye
DO ij = ijgs_e,ijge_e
DO ip = ipgs_e,ipge_e
p_int = parray_e(ip)
eo = MODULO(p_int,2) ! Indicates if we are on even or odd z grid
! Compute z derivatives
- CALL grad_z(eo,nadiab_moments_e(ip,ij,ikx,iky,izgs:izge), ddz_nepj(ip,ij,ikx,iky,izs:ize))
+ CALL grad_z(eo,nadiab_moments_e(ip,ij,iky,ikx,izgs:izge), ddz_nepj(ip,ij,iky,ikx,izs:ize))
! Parallel hyperdiffusion
- CALL grad_z2(moments_e(ip,ij,ikx,iky,izgs:izge,updatetlevel),ddz2_Nepj(ip,ij,ikx,iky,izs:ize))
+ CALL grad_z2(moments_e(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddz2_Nepj(ip,ij,iky,ikx,izs:ize))
! Compute even odd grids interpolation
- CALL interp_z(eo,nadiab_moments_e(ip,ij,ikx,iky,izgs:izge), interp_nepj(ip,ij,ikx,iky,izs:ize))
+ CALL interp_z(eo,nadiab_moments_e(ip,ij,iky,ikx,izgs:izge), interp_nepj(ip,ij,iky,ikx,izs:ize))
ENDDO
ENDDO
ENDDO
ENDDO
ENDIF
- DO iky = ikys,ikye
- DO ikx = ikxs,ikxe
+ DO ikx = ikxs,ikxe
+ DO iky = ikys,ikye
DO ij = ijgs_i,ijge_i
DO ip = ipgs_i,ipge_i
p_int = parray_i(ip)
eo = MODULO(p_int,2) ! Indicates if we are on even or odd z grid
! Compute z first derivative
- CALL grad_z(eo,nadiab_moments_i(ip,ij,ikx,iky,izgs:izge), ddz_nipj(ip,ij,ikx,iky,izs:ize))
+ CALL grad_z(eo,nadiab_moments_i(ip,ij,iky,ikx,izgs:izge), ddz_nipj(ip,ij,iky,ikx,izs:ize))
! Parallel hyperdiffusion
- CALL grad_z2(moments_i(ip,ij,ikx,iky,izgs:izge,updatetlevel),ddz2_Nipj(ip,ij,ikx,iky,izs:ize))
+ CALL grad_z2(moments_i(ip,ij,iky,ikx,izgs:izge,updatetlevel),ddz2_Nipj(ip,ij,iky,ikx,izs:ize))
! Compute even odd grids interpolation
- CALL interp_z(eo,nadiab_moments_i(ip,ij,ikx,iky,izgs:izge), interp_nipj(ip,ij,ikx,iky,izs:ize))
+ CALL interp_z(eo,nadiab_moments_i(ip,ij,iky,ikx,izgs:izge), interp_nipj(ip,ij,iky,ikx,izs:ize))
ENDDO
ENDDO
ENDDO
@@ -268,16 +268,16 @@ SUBROUTINE compute_density
! electron density
dens = 0._dp
DO ij = ijs_e, ije_e
- dens = dens + kernel_e(ij,ikx,iky,iz,0) * nadiab_moments_e(ip0_e,ij,ikx,iky,iz)
+ dens = dens + kernel_e(ij,iky,ikx,iz,0) * nadiab_moments_e(ip0_e,ij,iky,ikx,iz)
ENDDO
- dens_e(ikx,iky,iz) = dens
+ dens_e(iky,ikx,iz) = dens
ENDIF
! ion density
dens = 0._dp
DO ij = ijs_i, ije_i
- dens = dens + kernel_i(ij,ikx,iky,iz,0) * nadiab_moments_i(ip0_e,ij,ikx,iky,iz)
+ dens = dens + kernel_i(ij,iky,ikx,iz,0) * nadiab_moments_i(ip0_e,ij,iky,ikx,iz)
ENDDO
- dens_i(ikx,iky,iz) = dens
+ dens_i(iky,ikx,iz) = dens
ENDDO
ENDDO
ENDDO
@@ -303,18 +303,18 @@ SUBROUTINE compute_uperp
! electron
uperp = 0._dp
DO ij = ijs_e, ije_e
- uperp = uperp + kernel_e(ij,ikx,iky,iz,0) *&
- 0.5_dp*(nadiab_moments_e(ip0_e,ij,ikx,iky,iz) - nadiab_moments_e(ip0_e,ij-1,ikx,iky,iz))
+ uperp = uperp + kernel_e(ij,iky,ikx,iz,0) *&
+ 0.5_dp*(nadiab_moments_e(ip0_e,ij,iky,ikx,iz) - nadiab_moments_e(ip0_e,ij-1,iky,ikx,iz))
ENDDO
- uper_e(ikx,iky,iz) = uperp
+ uper_e(iky,ikx,iz) = uperp
ENDIF
! ion
uperp = 0._dp
DO ij = ijs_i, ije_i
- uperp = uperp + kernel_i(ij,ikx,iky,iz,0) *&
- 0.5_dp*(nadiab_moments_i(ip0_i,ij,ikx,iky,iz) - nadiab_moments_i(ip0_i,ij-1,ikx,iky,iz))
+ uperp = uperp + kernel_i(ij,iky,ikx,iz,0) *&
+ 0.5_dp*(nadiab_moments_i(ip0_i,ij,iky,ikx,iz) - nadiab_moments_i(ip0_i,ij-1,iky,ikx,iz))
ENDDO
- uper_i(ikx,iky,iz) = uperp
+ uper_i(iky,ikx,iz) = uperp
ENDDO
ENDDO
ENDDO
@@ -339,16 +339,16 @@ SUBROUTINE compute_upar
! electron
upar = 0._dp
DO ij = ijs_e, ije_e
- upar = upar + kernel_e(ij,ikx,iky,iz,1)*nadiab_moments_e(ip1_e,ij,ikx,iky,iz)
+ upar = upar + kernel_e(ij,iky,ikx,iz,1)*nadiab_moments_e(ip1_e,ij,iky,ikx,iz)
ENDDO
- upar_e(ikx,iky,iz) = upar
+ upar_e(iky,ikx,iz) = upar
ENDIF
! ion
upar = 0._dp
DO ij = ijs_i, ije_i
- upar = upar + kernel_i(ij,ikx,iky,iz,1)*nadiab_moments_i(ip1_i,ij,ikx,iky,iz)
+ upar = upar + kernel_i(ij,iky,ikx,iz,1)*nadiab_moments_i(ip1_i,ij,iky,ikx,iz)
ENDDO
- upar_i(ikx,iky,iz) = upar
+ upar_i(iky,ikx,iz) = upar
ENDDO
ENDDO
ENDDO
@@ -381,23 +381,23 @@ SUBROUTINE compute_tperp
Tperp = 0._dp
DO ij = ijs_e, ije_e
j_dp = REAL(ij-1,dp)
- Tperp = Tperp + kernel_e(ij,ikx,iky,iz,0)*&
- ((2_dp*j_dp+1)*nadiab_moments_e(ip0_e,ij ,ikx,iky,iz)&
- -j_dp *nadiab_moments_e(ip0_e,ij-1,ikx,iky,iz)&
- -j_dp+1 *nadiab_moments_e(ip0_e,ij+1,ikx,iky,iz))
+ Tperp = Tperp + kernel_e(ij,iky,ikx,iz,0)*&
+ ((2_dp*j_dp+1)*nadiab_moments_e(ip0_e,ij ,iky,ikx,iz)&
+ -j_dp *nadiab_moments_e(ip0_e,ij-1,iky,ikx,iz)&
+ -j_dp+1 *nadiab_moments_e(ip0_e,ij+1,iky,ikx,iz))
ENDDO
- Tper_e(ikx,iky,iz) = Tperp
+ Tper_e(iky,ikx,iz) = Tperp
ENDIF
! ion temperature
Tperp = 0._dp
DO ij = ijs_i, ije_i
j_dp = REAL(ij-1,dp)
- Tperp = Tperp + kernel_i(ij,ikx,iky,iz,0)*&
- ((2_dp*j_dp+1)*nadiab_moments_i(ip0_i,ij ,ikx,iky,iz)&
- -j_dp *nadiab_moments_i(ip0_i,ij-1,ikx,iky,iz)&
- -j_dp+1 *nadiab_moments_i(ip0_i,ij+1,ikx,iky,iz))
+ Tperp = Tperp + kernel_i(ij,iky,ikx,iz,0)*&
+ ((2_dp*j_dp+1)*nadiab_moments_i(ip0_i,ij ,iky,ikx,iz)&
+ -j_dp *nadiab_moments_i(ip0_i,ij-1,iky,ikx,iz)&
+ -j_dp+1 *nadiab_moments_i(ip0_i,ij+1,iky,ikx,iz))
ENDDO
- Tper_i(ikx,iky,iz) = Tperp
+ Tper_i(iky,ikx,iz) = Tperp
ENDDO
ENDDO
ENDDO
@@ -426,21 +426,21 @@ SUBROUTINE compute_Tpar
Tpar = 0._dp
DO ij = ijs_e, ije_e
j_dp = REAL(ij-1,dp)
- Tpar = Tpar + kernel_e(ij,ikx,iky,iz,0)*&
- (SQRT2 * nadiab_moments_e(ip2_e,ij,ikx,iky,iz) &
- + nadiab_moments_e(ip0_e,ij,ikx,iky,iz))
+ Tpar = Tpar + kernel_e(ij,iky,ikx,iz,0)*&
+ (SQRT2 * nadiab_moments_e(ip2_e,ij,iky,ikx,iz) &
+ + nadiab_moments_e(ip0_e,ij,iky,ikx,iz))
ENDDO
- Tpar_e(ikx,iky,iz) = Tpar
+ Tpar_e(iky,ikx,iz) = Tpar
ENDIF
! ion temperature
Tpar = 0._dp
DO ij = ijs_i, ije_i
j_dp = REAL(ij-1,dp)
- Tpar = Tpar + kernel_i(ij,ikx,iky,iz,0)*&
- (SQRT2 * nadiab_moments_i(ip2_i,ij,ikx,iky,iz) &
- + nadiab_moments_i(ip0_i,ij,ikx,iky,iz))
+ Tpar = Tpar + kernel_i(ij,iky,ikx,iz,0)*&
+ (SQRT2 * nadiab_moments_i(ip2_i,ij,iky,ikx,iz) &
+ + nadiab_moments_i(ip0_i,ij,iky,ikx,iz))
ENDDO
- Tpar_i(ikx,iky,iz) = Tpar
+ Tpar_i(iky,ikx,iz) = Tpar
ENDDO
ENDDO
ENDDO
diff --git a/src/restarts_mod.F90 b/src/restarts_mod.F90
index 613cf37f8175cc1b1a118e7fb867b39ee48f11ec..ea1697c311180786225e5ce74880114b9108f669 100644
--- a/src/restarts_mod.F90
+++ b/src/restarts_mod.F90
@@ -72,10 +72,10 @@ CONTAINS
! Read state of system from checkpoint file
IF (KIN_E) THEN
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_
- CALL getarrnd(fidrst, dset_name, moments_e(ips_e:ipe_e, ijs_e:ije_e, ikxs:ikxe, ikys:ikye, izs:ize, 1),(/1,3,5/))
+ CALL getarrnd(fidrst, dset_name, moments_e(ips_e:ipe_e, ijs_e:ije_e, ikys:ikye, ikxs:ikxe, izs:ize, 1),(/1,3,5/))
ENDIF
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_
- CALL getarrnd(fidrst, dset_name, moments_i(ips_i:ipe_i, ijs_i:ije_i, ikxs:ikxe, ikys:ikye, izs:ize, 1),(/1,3,5/))
+ CALL getarrnd(fidrst, dset_name, moments_i(ips_i:ipe_i, ijs_i:ije_i, ikys:ikye, ikxs:ikxe, izs:ize, 1),(/1,3,5/))
CALL closef(fidrst)
@@ -113,7 +113,7 @@ CONTAINS
CALL getatt(fidrst,"/data/input/" , "start_iframe5d", n0)
! Allocate the required size to load checkpoints moments
! CALL allocate_array(moments_e_cp, 1,pmaxe_cp+1, 1,jmaxe_cp+1, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array(moments_e_cp, 1,pmaxe_cp+1, 1,jmaxe_cp+1, 1,Nkx_cp, 1,Nky_cp, 1,Nz_cp)
+ CALL allocate_array(moments_e_cp, 1,pmaxe_cp+1, 1,jmaxe_cp+1, 1,Nky_cp, 1,Nkx_cp, 1,Nz_cp)
! Find the last results of the checkpoint file by iteration
n_ = n0+1
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_ ! start with moments_e/000001
@@ -124,8 +124,8 @@ CONTAINS
n_ = n_ - 1 ! n_ is not a file so take the previous one n_-1
! Read state of system from checkpoint file and load every moment to change the distribution
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_
- ! CALL getarrnd(fidrst, dset_name, moments_e_cp(1:pmaxe_cp+1, 1:jmaxe_cp+1, ikxs:ikxe, ikys:ikye, izs:ize),(/1,3/))
- CALL getarr(fidrst, dset_name, moments_e_cp(1:pmaxe_cp+1, 1:jmaxe_cp+1, 1:Nkx_cp, 1:Nky_cp, 1:Nz_cp))
+ ! CALL getarrnd(fidrst, dset_name, moments_e_cp(1:pmaxe_cp+1, 1:jmaxe_cp+1, ikys:ikye, ikxs:ikxe, izs:ize),(/1,3/))
+ CALL getarr(fidrst, dset_name, moments_e_cp(1:pmaxe_cp+1, 1:jmaxe_cp+1, 1:Nky_cp, 1:Nkx_cp, 1:Nz_cp))
! Initialize simulation moments array with checkpoints ones
! (they may have a larger number of polynomials, set to 0 at the begining)
moments_e = 0._dp;
@@ -138,7 +138,7 @@ CONTAINS
DO ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz = izs,ize
- moments_e(ip,ij,ikx,iky,iz,:) = moments_e_cp(ip,ij,ikx,iky,iz)
+ moments_e(ip,ij,iky,ikx,iz,:) = moments_e_cp(ip,ij,iky,ikx,iz)
ENDDO
ENDDO
ENDDO
@@ -158,7 +158,7 @@ CONTAINS
CALL getatt(fidrst,"/data/input/" , "start_iframe5d", n0)
! Allocate the required size to load checkpoints moments
! CALL allocate_array(moments_i_cp, 1,pmaxi_cp+1, 1,jmaxi_cp+1, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array(moments_i_cp, 1,pmaxi_cp+1, 1,jmaxi_cp+1, 1,Nkx_cp, 1,Nky_cp, 1,Nz_cp)
+ CALL allocate_array(moments_i_cp, 1,pmaxi_cp+1, 1,jmaxi_cp+1, 1,Nky_cp, 1,Nkx_cp, 1,Nz_cp)
! Find the last results of the checkpoint file by iteration
n_ = n0+1
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_ ! start with moments_e/000001
@@ -170,8 +170,8 @@ CONTAINS
! Read state of system from checkpoint file and load every moment to change the distribution
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_
- ! CALL getarrnd(fidrst, dset_name, moments_i_cp(1:pmaxi_cp+1, 1:jmaxi_cp+1, ikxs:ikxe, ikys:ikye, izs:ize),(/1,3/))
- CALL getarr(fidrst, dset_name, moments_i_cp(1:pmaxi_cp+1, 1:jmaxi_cp+1, 1:Nkx_cp, 1:Nky_cp, 1:Nz_cp))
+ ! CALL getarrnd(fidrst, dset_name, moments_i_cp(1:pmaxi_cp+1, 1:jmaxi_cp+1, ikys:ikye, ikxs:ikxe, izs:ize),(/1,3/))
+ CALL getarr(fidrst, dset_name, moments_i_cp(1:pmaxi_cp+1, 1:jmaxi_cp+1, 1:Nky_cp, 1:Nkx_cp, 1:Nz_cp))
! Initialize simulation moments array with checkpoints ones
! (they may have a larger number of polynomials, set to 0 at the begining)
@@ -185,7 +185,7 @@ CONTAINS
DO ikx=ikxs,ikxe
DO iky=ikys,ikye
DO iz = izs,ize
- moments_i(ip,ij,ikx,iky,iz,:) = moments_i_cp(ip,ij,ikx,iky,iz)
+ moments_i(ip,ij,iky,ikx,iz,:) = moments_i_cp(ip,ij,iky,ikx,iz)
ENDDO
ENDDO
ENDDO
diff --git a/src/srcinfo.h b/src/srcinfo.h
index 37401921917253200fe476e0e73ef57cad15ad0f..1521a89638fe2b6222543ddba406baaf4ff274c5 100644
--- a/src/srcinfo.h
+++ b/src/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='7642d96-dirty')
-parameter (BRANCH='master')
+parameter (VERSION='6458026-dirty')
+parameter (BRANCH='kx_pos_plane')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Thu Apr 28 15:46:58 CEST 2022')
+parameter (EXECDATE='Tue May 3 10:57:28 CEST 2022')
parameter (HOST ='spcpc606')
diff --git a/src/srcinfo/srcinfo.h b/src/srcinfo/srcinfo.h
index 37401921917253200fe476e0e73ef57cad15ad0f..1521a89638fe2b6222543ddba406baaf4ff274c5 100644
--- a/src/srcinfo/srcinfo.h
+++ b/src/srcinfo/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='7642d96-dirty')
-parameter (BRANCH='master')
+parameter (VERSION='6458026-dirty')
+parameter (BRANCH='kx_pos_plane')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Thu Apr 28 15:46:58 CEST 2022')
+parameter (EXECDATE='Tue May 3 10:57:28 CEST 2022')
parameter (HOST ='spcpc606')
diff --git a/src/stepon.F90 b/src/stepon.F90
index d928a4b897f7e93a0eea422acc4fc40eba62e292..230f6bf0c47e6b95682c7174f0bb435e8a65a1bd 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -102,11 +102,11 @@ SUBROUTINE stepon
SUBROUTINE anti_aliasing
IF(KIN_E)THEN
DO iz=izgs,izge
- DO iky=ikys,ikye
- DO ikx=ikxs,ikxe
+ DO ikx=ikxs,ikxe
+ DO iky=ikys,ikye
DO ij=ijgs_e,ijge_e
DO ip=ipgs_e,ipge_e
- moments_e( ip,ij,ikx,iky,iz,:) = AA_x(ikx)* AA_y(iky) * moments_e( ip,ij,ikx,iky,iz,:)
+ moments_e( ip,ij,iky,ikx,iz,:) = AA_x(ikx)* AA_y(iky) * moments_e( ip,ij,iky,ikx,iz,:)
END DO
END DO
END DO
@@ -114,11 +114,11 @@ SUBROUTINE stepon
END DO
ENDIF
DO iz=izgs,izge
- DO iky=ikys,ikye
- DO ikx=ikxs,ikxe
+ DO ikx=ikxs,ikxe
+ DO iky=ikys,ikye
DO ij=ijs_i,ije_i
DO ip=ipgs_i,ipge_i
- moments_i( ip,ij,ikx,iky,iz,:) = AA_x(ikx)* AA_y(iky) * moments_i( ip,ij,ikx,iky,iz,:)
+ moments_i( ip,ij,iky,ikx,iz,:) = AA_x(ikx)* AA_y(iky) * moments_i( ip,ij,iky,ikx,iz,:)
END DO
END DO
END DO
@@ -127,17 +127,17 @@ SUBROUTINE stepon
END SUBROUTINE anti_aliasing
SUBROUTINE enforce_symmetry ! Force X(k) = X(N-k)* complex conjugate symmetry
- IF ( contains_kx0 ) THEN
+ IF ( contains_ky0 ) THEN
! Electron moments
IF(KIN_E) THEN
DO iz=izs,ize
DO ij=ijgs_e,ijge_e
DO ip=ipgs_e,ipge_e
- DO iky=2,Nky/2 !symmetry at kx = 0
- moments_e( ip,ij,ikx_0,iky,iz, :) = CONJG(moments_e( ip,ij,ikx_0,Nky+2-iky,iz, :))
+ DO ikx=2,Nkx/2 !symmetry at ky = 0
+ moments_e( ip,ij,iky_0,ikx,iz, :) = CONJG(moments_e( ip,ij,iky_0,Nkx+2-iky,iz, :))
END DO
! must be real at origin
- moments_e(ip,ij, ikx_0,iky_0,iz, :) = REAL(moments_e(ip,ij, ikx_0,iky_0,iz, :))
+ moments_e(ip,ij, iky_0,ikx_0,iz, :) = REAL(moments_e(ip,ij, iky_0,ikx_0,iz, :))
END DO
END DO
END DO
@@ -146,20 +146,20 @@ SUBROUTINE stepon
DO iz=izs,ize
DO ij=ijgs_i,ijge_i
DO ip=ipgs_i,ipge_i
- DO iky=2,Nky/2 !symmetry at kx = 0
- moments_i( ip,ij,ikx_0,iky,iz, :) = CONJG(moments_i( ip,ij,ikx_0,Nky+2-iky,iz, :))
+ DO ikx=2,Nkx/2 !symmetry at ky = 0
+ moments_i( ip,ij,iky_0,ikx,iz, :) = CONJG(moments_i( ip,ij,iky_0,Nkx+2-ikx,iz, :))
END DO
! must be real at origin and top right
- moments_i(ip,ij, ikx_0,iky_0,iz, :) = REAL(moments_i(ip,ij, ikx_0,iky_0,iz, :))
+ moments_i(ip,ij, iky_0,ikx_0,iz, :) = REAL(moments_i(ip,ij, iky_0,ikx_0,iz, :))
END DO
END DO
END DO
! Phi
- DO iky=2,Nky/2 !symmetry at kx = 0
- phi(ikx_0,iky,izs:ize) = phi(ikx_0,Nky+2-iky,izs:ize)
+ DO ikx=2,Nkx/2 !symmetry at ky = 0
+ phi(iky_0,ikx,izs:ize) = phi(iky_0,Nkx+2-ikx,izs:ize)
END DO
! must be real at origin
- phi(ikx_0,iky_0,izs:ize) = REAL(phi(ikx_0,iky_0,izs:ize))
+ phi(iky_0,ikx_0,izs:ize) = REAL(phi(iky_0,ikx_0,izs:ize))
ENDIF
END SUBROUTINE enforce_symmetry
diff --git a/src/utility_mod.F90 b/src/utility_mod.F90
index bd3ab13b13a5322f9d3987c41316f2baba98e5b0..b6e989a746f8dedcdf7f634ee7d769835060af71 100644
--- a/src/utility_mod.F90
+++ b/src/utility_mod.F90
@@ -59,7 +59,7 @@ CONTAINS
use prec_const
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(ikxs:ikxe,ikys:ikye), INTENT(IN) :: field
+ COMPLEX(dp), DIMENSION(ikys:ikye,ikxs:ikxe), INTENT(IN) :: field
CHARACTER(LEN=*), INTENT(IN) :: str
LOGICAL :: mlend
COMPLEX(dp) :: sumfield
@@ -75,8 +75,8 @@ CONTAINS
SUBROUTINE manual_2D_bcast(field_)
USE grid
IMPLICIT NONE
- COMPLEX(dp), INTENT(INOUT) :: field_(ikxs:ikxe,ikys:ikye)
- COMPLEX(dp) :: buffer(ikxs:ikxe,ikys:ikye)
+ COMPLEX(dp), INTENT(INOUT) :: field_(ikys:ikye,ikxs:ikxe)
+ COMPLEX(dp) :: buffer(ikys:ikye,ikxs:ikxe)
INTEGER :: i_, root, world_rank, world_size
root = 0;
CALL MPI_COMM_RANK(comm_p,world_rank,ierr)
@@ -87,21 +87,21 @@ CONTAINS
! Fill the buffer
DO ikx = ikxs,ikxe
DO iky = ikys,ikye
- buffer(ikx,iky) = field_(ikx,iky)
+ buffer(iky,ikx) = field_(iky,ikx)
ENDDO
ENDDO
! Send it to all the other processes
DO i_ = 0,num_procs_p-1
IF (i_ .NE. world_rank) &
- CALL MPI_SEND(buffer, local_nkx * Nky , MPI_DOUBLE_COMPLEX, i_, 0, comm_p, ierr)
+ CALL MPI_SEND(buffer, local_nky * Nkx , MPI_DOUBLE_COMPLEX, i_, 0, comm_p, ierr)
ENDDO
ELSE
! Recieve buffer from root
- CALL MPI_RECV(buffer, local_nkx * Nky , MPI_DOUBLE_COMPLEX, root, 0, comm_p, MPI_STATUS_IGNORE, ierr)
+ CALL MPI_RECV(buffer, local_nky * Nkx , MPI_DOUBLE_COMPLEX, root, 0, comm_p, MPI_STATUS_IGNORE, ierr)
! Write it in phi
DO ikx = ikxs,ikxe
DO iky = ikys,ikye
- field_(ikx,iky) = buffer(ikx,iky)
+ field_(iky,ikx) = buffer(iky,ikx)
ENDDO
ENDDO
ENDIF
@@ -112,11 +112,11 @@ CONTAINS
SUBROUTINE manual_3D_bcast(field_)
USE grid
IMPLICIT NONE
- COMPLEX(dp), INTENT(INOUT) :: field_(ikxs:ikxe,ikys:ikye,izs:ize)
- COMPLEX(dp) :: buffer(ikxs:ikxe,ikys:ikye,izs:ize)
+ COMPLEX(dp), INTENT(INOUT) :: field_(ikys:ikye,ikxs:ikxe,izs:ize)
+ COMPLEX(dp) :: buffer(ikys:ikye,ikxs:ikxe,izs:ize)
INTEGER :: i_, root, world_rank, world_size, count
root = 0;
- count = (ikxe-ikxs+1) * (ikye-ikys+1) * (ize-izs+1);
+ count = (ikye-ikys+1) * (ikxe-ikxs+1) * (ize-izs+1);
CALL MPI_COMM_RANK(comm_p,world_rank,ierr)
CALL MPI_COMM_SIZE(comm_p,world_size,ierr)
@@ -124,11 +124,11 @@ SUBROUTINE manual_3D_bcast(field_)
!! Broadcast phi to the other processes on the same k range (communicator along p)
IF (world_rank .EQ. root) THEN
! Fill the buffer
- DO ikx = ikxs,ikxe
- DO iky = ikys,ikye
- DO iz = izs,ize
- buffer(ikx,iky,iz) = field_(ikx,iky,iz)
- ENDDO
+ DO iz = izs,ize
+ DO ikx = ikxs,ikxe
+ DO iky = ikys,ikye
+ buffer(iky,ikx,iz) = field_(iky,ikx,iz)
+ ENDDO
ENDDO
ENDDO
! Send it to all the other processes
@@ -140,10 +140,10 @@ SUBROUTINE manual_3D_bcast(field_)
! Recieve buffer from root
CALL MPI_RECV(buffer, count, MPI_DOUBLE_COMPLEX, root, 0, comm_p, MPI_STATUS_IGNORE, ierr)
! Write it in phi
- DO ikx = ikxs,ikxe
- DO iky = ikys,ikye
- DO iz = izs,ize
- field_(ikx,iky,iz) = buffer(ikx,iky,iz)
+ DO iz = izs,ize
+ DO ikx = ikxs,ikxe
+ DO iky = ikys,ikye
+ field_(iky,ikx,iz) = buffer(iky,ikx,iz)
ENDDO
ENDDO
ENDDO
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 2bcaece23872762f602ba06b3cd5c572c2b7a085..eba8d519a25149f53b33c9a05065c2deb10b5869 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -10,7 +10,7 @@ system(['mkdir -p ',MISCDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
system(CMD);
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 10;
+JOBNUMMIN = 00; JOBNUMMAX = 10;
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
@@ -50,7 +50,7 @@ options.PLAN = 'xy';
% options.PLAN = 'sx';
options.COMP = 'avg';
% options.TIME = dat.Ts5D;
-options.TIME = 0:1:80;
+options.TIME = 100:1:200;
data.EPS = 0.1;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
@@ -62,17 +62,17 @@ if 0
options.INTERP = 0;
options.POLARPLOT = 0;
options.AXISEQUAL = 1;
-% options.NAME = '\phi';
-options.NAME = 'n_i';
+options.NAME = '\phi';
+% options.NAME = 'n_i';
% options.NAME = 'N_i^{00}';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'xy';
+options.PLAN = 'kxky';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
-options.COMP = 8;
-options.TIME = [500 700 900];
+options.COMP = 1;
+options.TIME = [1];
data.a = data.EPS * 1000;
fig = photomaton(data,options);
save_figure(data,fig)
@@ -155,7 +155,7 @@ options.K2PLOT = 1;
options.TIME = 5:1:15;
options.NMA = 1;
options.NMODES = 5;
-options.iz = 8;
+options.iz = 1;
fig = mode_growth_meter(data,options);
save_figure(data,fig)
end
@@ -189,4 +189,4 @@ options.kzkx = 0;
options.kzky = 1;
[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
save_figure(data,fig)
-end
\ No newline at end of file
+end
diff --git a/wk/analysis_header.m b/wk/analysis_header.m
index f2149b60de00f70270205d3889ae3459c74c3c1f..3970fcc89874200ad09092a74155a4988fb67046 100644
--- a/wk/analysis_header.m
+++ b/wk/analysis_header.m
@@ -2,9 +2,11 @@
helazdir = '/home/ahoffman/HeLaZ/';
% Directory of the simulation (from results)
% if 1% Local results
-outfile ='2D_Zpinch/LDDK';
+outfile ='';
+outfile ='quick_run/64x64_5x3_L_120_kN_2.0_kT_0.5_nu_1e-01_SGGK';
+% outfile ='pedestal/64x64x16x2x1_L_300_LnT_20_nu_0.1';
% outfile ='quick_run/32x32x16_5x3_L_300_q0_2.5_e_0.18_kN_20_kT_20_nu_1e-01_DGGK';
% outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_1.0/';
% outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8_to_1.0/';
% outfile ='pedestal/128x128x16x4x2_L_120_LnT_40_nuDG_0.1';
-run analysis_3D
\ No newline at end of file
+run analysis_3D
diff --git a/wk/quick_run.m b/wk/quick_run.m
index 64400569c20e5f2a3d8288f3e907527402e34540..e7525829beff8ae20ce42b8b41387e759dc4ddea 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -15,8 +15,8 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_N = 1.8;%1.9;%2.22; % Density gradient drive
-K_T = 0;%0.25*K_N; % Temperature '''
+K_N = 1.9;%2.22; % Density gradient drive
+K_T = 0.25*K_N; % Temperature '''
K_E = 0.0; % Electrostat '''
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
@@ -49,10 +49,10 @@ SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
SIMID = 'quick_run'; % Name of the simulation
-LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+LINEARITY = 'nonlinear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO = 'DG';
+CO = 'SG';
GKCO = 1; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
@@ -73,11 +73,11 @@ HD_CO = 0.0; % Hyper diffusivity cutoff ratio
kmax = NX*pi/LX;% Highest fourier mode
MU = 0.0; % Hyperdiffusivity coefficient
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
-MU_X = MU; %
-MU_Y = MU; %
+MU_X = MU; %
+MU_Y = MU; %
MU_Z = 0.0; %
MU_P = 0.0; %
-MU_J = 0.0; %
+MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
@@ -89,7 +89,7 @@ setup
system(['rm fort*.90']);
% Run linear simulation
if RUN
- system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',HELAZDIR,'bin/',EXECNAME,' 2 1 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 3 2 0; cd ../../../wk'])
@@ -100,7 +100,7 @@ end
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [HELAZDIR,'results/',filename,'/'];
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 00;
+JOBNUMMIN = 00; JOBNUMMAX = 00;
data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
%% Short analysis