From 9df7bc48d8fc3aa5f7d1313f7ec39f10423def72 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 10 May 2022 09:40:35 +0200
Subject: [PATCH] added the degree closure on moments

---
 src/auxval.F90             |  12 +++
 src/closure_mod.F90        |  38 +++----
 src/moments_eq_rhs_mod.F90 | 216 ++++++++++++++++++-------------------
 src/processing_mod.F90     |   4 +-
 wk/analysis_3D.m           |  15 +--
 wk/analysis_header.m       |   5 +-
 6 files changed, 151 insertions(+), 139 deletions(-)

diff --git a/src/auxval.F90 b/src/auxval.F90
index 7d6fd969..db695e67 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -78,4 +78,16 @@ subroutine auxval
   ENDDO
   CALL mpi_barrier(MPI_COMM_WORLD, ierr)
 
+  IF((my_id.EQ.0) .AND. (CLOS .EQ. 1)) THEN
+  IF(KIN_E) &
+  write(*,*) 'Closure = 1 -> Maximal Nepj degree is min(Pmaxe,2*Jmaxe+1): De = ', dmaxi
+  write(*,*) 'Closure = 1 -> Maximal Nipj degree is min(Pmaxi,2*Jmaxi+1): Di = ', dmaxi
+  ENDIF
+  DO ip = ips_i,ipe_i
+    DO ij = ijs_i,ije_i
+      IF((parray_i(ip)+2*jarray_i(ij) .LE. dmaxi) .AND. (rank_ky + rank_z .EQ. 0))&
+      print*, '(',parray_i(ip),',',jarray_i(ij),')'
+    ENDDO
+  ENDDO
+
 END SUBROUTINE auxval
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index b316e442..5d48206c 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -23,29 +23,25 @@ SUBROUTINE apply_closure_model
 
 
   ELSEIF (CLOS .EQ. 1) THEN
-    ! Truncation at highest fully represented kinetic moment
+    ! zero truncation, An+1=0 for n+1>nmax only
+    CALL ghosts_upper_truncation
+    ! Additional truncation at highest fully represented kinetic moment
     ! e.g. Dmax = 3 means
-    ! all Napj s.t. p+2j <= 3
+    ! only Napj s.t. p+2j <= 3 are evolved
     ! -> (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
-          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,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,iky,ikx,iz,updatetlevel) = 0._dp
-            ENDDO
-          ENDDO
-        ENDDO
+    ! =>> Dmax = min(Pmax,Jmax+1)
+    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,ikys:ikye,ikxs:ikxe,izgs:izge,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,ikys:ikye,ikxs:ikxe,izgs:izge,updatetlevel) = 0._dp
       ENDDO
     ENDDO
     ! + ghosts truncation
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 2fe9694a..d2e316fe 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -60,61 +60,61 @@ SUBROUTINE moments_eq_rhs_e
             kperp2= kparray(iky,ikx,iz,eo)**2
 
           IF((CLOS .EQ. 1) .AND. (p_int+2*j_int .LE. dmaxe)) THEN
-          !! 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,iky,ikx,iz)
-          ! term propto n_e^{p+2,j}
-          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,iky,ikx,iz)
-          ! term propto n_e^{p,j+1}
-          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,iky,ikx,iz)
-          ! Parallel dynamic
-          ! ddz derivative for Landau damping term
-          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  ,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  ,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)
+            !! 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,iky,ikx,iz)
+            ! term propto n_e^{p+2,j}
+            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,iky,ikx,iz)
+            ! term propto n_e^{p,j+1}
+            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,iky,ikx,iz)
+            ! Parallel dynamic
+            ! ddz derivative for Landau damping term
+            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  ,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  ,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  ,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
+            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  ,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,iky,ikx,iz,updatetlevel) = &
-              ! Perpendicular magnetic gradient/curvature effects
-              - 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)
-              - gradzB(iz,eo)* Tmir  *gradz_coeff(iz,eo) &
-              ! 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,iky,ikx,iz,updatetlevel) &
-              ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)"
-              + mu_z * diff_dz_coeff * ddz2_Nepj(ip,ij,iky,ikx,iz) &
-              ! Collision term
-              + TColl_e(ip,ij,iky,ikx,iz) &
-              ! Nonlinear term
-              - Sepj(ip,ij,iky,ikx,iz)
+            !! Sum of all RHS terms
+            moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = &
+                ! Perpendicular magnetic gradient/curvature effects
+                - 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)
+                - gradzB(iz,eo)* Tmir  *gradz_coeff(iz,eo) &
+                ! 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,iky,ikx,iz,updatetlevel) &
+                ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)"
+                + mu_z * diff_dz_coeff * ddz2_Nepj(ip,ij,iky,ikx,iz) &
+                ! Collision term
+                + TColl_e(ip,ij,iky,ikx,iz) &
+                ! Nonlinear term
+                - Sepj(ip,ij,iky,ikx,iz)
           ELSE
             moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp
           ENDIF
@@ -180,64 +180,64 @@ SUBROUTINE moments_eq_rhs_i
             eo    = MODULO(p_int,2) ! Indicates if we are on odd or even z grid
             kperp2= kparray(iky,ikx,iz,eo)**2
             IF((CLOS .EQ. 1) .AND. (p_int+2*j_int .LE. dmaxi)) THEN
-            !! 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  ,iky,ikx,iz)
-            ! term propto n_i^{p+2,j}
-            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  ,iky,ikx,iz)
-            ! term propto n_i^{p,j+1}
-            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,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,iky,ikx,iz) &
-                      + xnipm1j(ip) * ddz_nipj(ip-1,ij,iky,ikx,iz)
-            ! Mirror terms
-            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  ,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)
+              !! 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  ,iky,ikx,iz)
+              ! term propto n_i^{p+2,j}
+              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  ,iky,ikx,iz)
+              ! term propto n_i^{p,j+1}
+              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,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,iky,ikx,iz) &
+                        + xnipm1j(ip) * ddz_nipj(ip-1,ij,iky,ikx,iz)
+              ! Mirror terms
+              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  ,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
+              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  ,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
+              !! Electrical potential term
+              IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
+                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,iky,ikx,iz,updatetlevel) = &
-                ! Perpendicular magnetic gradient/curvature effects
-                - imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp &
-                ! Parallel coupling (Landau damping)
-                - gradz_coeff(iz,eo) * Tpar &
-                ! Mirror term (parallel magnetic gradient)
-                - gradzB(iz,eo) * gradz_coeff(iz,eo) * Tmir &
-                ! 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,iky,ikx,iz,updatetlevel) &
-                ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)"
-                + mu_z * diff_dz_coeff * ddz2_Nipj(ip,ij,iky,ikx,iz) &
-                ! Collision term
-                + TColl_i(ip,ij,iky,ikx,iz)&
-                ! Nonlinear term
-                - Sipj(ip,ij,iky,ikx,iz)
+              !! Sum of all RHS terms
+              moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = &
+                  ! Perpendicular magnetic gradient/curvature effects
+                  - imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp &
+                  ! Parallel coupling (Landau damping)
+                  - gradz_coeff(iz,eo) * Tpar &
+                  ! Mirror term (parallel magnetic gradient)
+                  - gradzB(iz,eo) * gradz_coeff(iz,eo) * Tmir &
+                  ! 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,iky,ikx,iz,updatetlevel) &
+                  ! Numerical parallel hyperdiffusion "+ (mu_z*kz**4)"
+                  + mu_z * diff_dz_coeff * ddz2_Nipj(ip,ij,iky,ikx,iz) &
+                  ! Collision term
+                  + TColl_i(ip,ij,iky,ikx,iz)&
+                  ! Nonlinear term
+                  - Sipj(ip,ij,iky,ikx,iz)
           ELSE
             moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp
           ENDIF
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index d4997a59..d98334cb 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -170,10 +170,10 @@ SUBROUTINE compute_nadiab_moments_z_gradients_and_interp
                                ddz_nepj, ddz2_Nepj, interp_nepj,&
                                ddz_nipj, ddz2_Nipj, interp_nipj
   USE time_integration, ONLY : updatetlevel
-  USE model,            ONLY : qe_taue, qi_taui, KIN_E
+  USE model,            ONLY : qe_taue, qi_taui, KIN_E, CLOS
   USE calculus,         ONLY : grad_z, grad_z2, interp_z
   IMPLICIT NONE
-  INTEGER :: eo, p_int
+  INTEGER :: eo, p_int, j_int
   CALL cpu_time(t0_process)
 
   ! Electron non adiab moments
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 2add3a76..0f2cdef4 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -24,7 +24,8 @@ if 1
 options.TAVG_0   = 0.8*data.Ts3D(end); 
 options.TAVG_1   = data.Ts3D(end); % Averaging times duration
 options.NMVA     = 1;              % Moving average for time traces
-options.ST_FIELD = '\Gamma_x';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+% options.ST_FIELD = '\Gamma_x';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+options.ST_FIELD = '\phi';          % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
 options.INTERP   = 1;
 fig = plot_radial_transport_and_spacetime(data,options);
 save_figure(data,fig)
@@ -41,18 +42,18 @@ if 0
 % Options
 options.INTERP    = 1;
 options.POLARPLOT = 0;
-% options.NAME      = '\phi';
+options.NAME      = '\phi';
 % options.NAME      = 'N_i^{00}';
 % options.NAME      = 'v_y';
 % options.NAME      = 'n_i^{NZ}';
-options.NAME      = '\Gamma_x';
+% options.NAME      = '\Gamma_x';
 % options.NAME      = 'n_i';
-options.PLAN      = 'kxky';
+options.PLAN      = 'xy';
 % options.NAME      = 'f_e';
 % options.PLAN      = 'sx';
 options.COMP      = 'avg';
 % options.TIME      = dat.Ts5D;
-options.TIME      = 00:1:200;
+options.TIME      = 00:1:250;
 data.EPS          = 0.1;
 data.a = data.EPS * 2000;
 create_film(data,options,'.gif')
@@ -99,7 +100,7 @@ options.XPERP     = linspace( 0,6,64);
 % options.SPAR      = vp';
 % options.XPERP     = mu';
 options.Z         = 1;
-options.T         = 5500;
+options.T         = 200;
 options.CTR       = 1;
 options.ONED      = 0;
 fig = plot_fa(data,options);
@@ -109,7 +110,7 @@ end
 if 0
 %% Hermite-Laguerre spectrum
 % options.TIME = 'avg';
-options.P2J  = 0;
+options.P2J  = 1;
 options.ST   = 0;
 options.NORMALIZED = 0;
 fig = show_moments_spectrum(data,options);
diff --git a/wk/analysis_header.m b/wk/analysis_header.m
index b3cf14ff..5188876d 100644
--- a/wk/analysis_header.m
+++ b/wk/analysis_header.m
@@ -3,10 +3,13 @@ helazdir = '/home/ahoffman/HeLaZ/';
 % Directory of the simulation (from results)
 % if 1% Local results
 outfile ='';
+outfile ='';
+outfile ='';
+outfile ='shearless_cyclone/128x128x16xdmax_L_120_CBC_1.0';
 % outfile ='quick_run/CLOS_1_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/128x128x16x8x4_L_120_CTC_1.0/';
+% outfile ='shearless_cyclone/128x128x16x8x4_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
-- 
GitLab