From 3e3e448161113b349573826ef3272cc2d0301ab7 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Fri, 30 Oct 2020 15:48:39 +0100
Subject: [PATCH] Utility changes

---
 src/advance_field.F90        |  2 +-
 src/advance_field_adapt.F90  |  2 +-
 src/auxval.F90               | 10 ++++----
 src/basic_mod.F90            | 30 +++++++++++++++++++++-
 src/control.F90              | 16 ++++++++++--
 src/diagnose.F90             | 18 +++++---------
 src/diagnostics_par_mod.F90  |  8 +++---
 src/fourier_mod.F90          | 24 +++++++++---------
 src/grid_mod.F90             |  4 +--
 src/inital.F90               | 48 ++++++++++++------------------------
 src/moments_eq_rhs.F90       |  6 ++---
 src/poisson.F90              |  2 +-
 src/prec_const_mod.F90       |  6 ++---
 src/readinputs.F90           |  5 ++++
 src/stepon.F90               |  2 --
 src/time_integration_mod.F90 |  2 +-
 wk/fort.90                   | 24 +++++++++---------
 wk/parameters_ZP.m           | 41 +++++++++++++++---------------
 18 files changed, 135 insertions(+), 115 deletions(-)

diff --git a/src/advance_field.F90 b/src/advance_field.F90
index c615d87b..0c27b593 100644
--- a/src/advance_field.F90
+++ b/src/advance_field.F90
@@ -12,7 +12,7 @@ CONTAINS
     use prec_const
     IMPLICIT NONE
 
-    call set_updatetlevel(mod(updatetlevel,ntimelevel)+1)
+    CALL set_updatetlevel(mod(updatetlevel,ntimelevel)+1)
 
   END SUBROUTINE advance_time_level
 
diff --git a/src/advance_field_adapt.F90 b/src/advance_field_adapt.F90
index d2d18f26..23809044 100644
--- a/src/advance_field_adapt.F90
+++ b/src/advance_field_adapt.F90
@@ -12,7 +12,7 @@ CONTAINS
     use prec_const
     IMPLICIT NONE
 
-    call set_updatetlevel(mod(updatetlevel,ntimelevel)+1)
+    CALL set_updatetlevel(mod(updatetlevel,ntimelevel)+1)
 
   END SUBROUTINE advance_time_level
 
diff --git a/src/auxval.F90 b/src/auxval.F90
index f0c451c9..c0cb1b6d 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -11,11 +11,11 @@ subroutine auxval
   INTEGER :: irows,irowe, irow, icol
   WRITE(*,*) '=== Set auxiliary values ==='
 
-  call set_rgrid
-  call set_zgrid
-  call set_krgrid
-  call set_kzgrid
-  call set_pj
+  CALL set_rgrid
+  CALL set_zgrid
+  CALL set_krgrid
+  CALL set_kzgrid
+  CALL set_pj
 
   CALL memory ! Allocate memory for global arrays
 END SUBROUTINE auxval
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index 735a3b45..eb0541ed 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -26,6 +26,7 @@ MODULE basic
   INTEGER :: lu_job  = 91              ! myjob file
 
   real :: start, finish ! To measure computation time
+  real :: start_est, finish_est, time_est
 
   INTERFACE allocate_array
     MODULE PROCEDURE allocate_array_dp1,allocate_array_dp2,allocate_array_dp3,allocate_array_dp4
@@ -66,9 +67,36 @@ CONTAINS
     !
   END SUBROUTINE daytim
   !================================================================================
+  SUBROUTINE display_h_min_s(time)
+    real :: time
+    integer :: days, hours, mins, secs
+    days = FLOOR(time/24./3600.);
+    hours= FLOOR(time/3600.);
+    mins = FLOOR(time/60.);
+    secs = FLOOR(time);
+
+    IF ( days .GT. 0 ) THEN !display day h min s
+      hours = (time/3600./24. - days) * 24
+      mins  = (time/3600. - days*24. - hours) * 60
+      secs  = (time/60. - days*24.*60 - hours*60 - mins) * 60
+      WRITE(*,*) 'CPU Time = ', days, '[day]', hours, '[h]', mins, '[min]', secs, '[s]'
+
+    ELSEIF ( hours .GT. 0 ) THEN !display h min s
+      mins  = (time/3600. - hours) * 60
+      secs  = (time/60. - hours*60 - mins) * 60
+      WRITE(*,*) 'CPU Time = ', hours, '[h]', mins, '[min]', secs, '[s]'
+
+    ELSEIF ( mins .GT. 0 ) THEN !display min s
+      secs  = (time/60. - mins) * 60
+      WRITE(*,*) 'CPU Time = ', mins, '[min]', secs, '[s]'
+
+    ELSE ! display s
+      WRITE(*,*) 'CPU Time = ', FLOOR(time), '[s]'
+    ENDIF
+  END SUBROUTINE display_h_min_s
+  !================================================================================
 
 ! To allocate arrays of doubles, integers, etc. at run time
-
   SUBROUTINE allocate_array_dp1(a,is1,ie1)
     IMPLICIT NONE
     real(dp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: a
diff --git a/src/control.F90 b/src/control.F90
index 59abdcf2..92638b01 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -5,7 +5,7 @@ SUBROUTINE control
   use prec_const
   IMPLICIT NONE
 
-  call cpu_time(start)
+  CALL cpu_time(start)
   !________________________________________________________________________________
   !              1.   Prologue
   !                   1.1     Initialize the parallel environment
@@ -14,7 +14,7 @@ SUBROUTINE control
   WRITE(*,'(a/)') '...MPI initialized'
 
 
-  call daytim('Start at ')
+  CALL daytim('Start at ')
 
   !                   1.2     Define data specific to run
   WRITE(*,*) 'Load basic data...'
@@ -37,6 +37,18 @@ SUBROUTINE control
   CALL inital
   WRITE(*,'(a/)') '...initial state created'
 
+  ! !                   1.5.1   Computation time estimation
+  ! WRITE(*,*) 'Estimation of computation time...'
+  ! CALL cpu_time(start_est)
+  ! CALL stepon
+  ! CALL cpu_time(finish_est)
+  ! time_est = 4.*tmax/dt*(finish_est-start_est)
+  ! CALL display_h_min_s(time_est)
+  ! WRITE(*,*) '... reinitializing'
+  ! CALL inital
+  ! WRITE(*,'(a/)') '... done'
+
+
   !                   1.6     Initial diagnostics
   WRITE(*,*) 'Initial diagnostics...'
   CALL diagnose(0)
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index e975693f..8c91e560 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -40,7 +40,7 @@ SUBROUTINE diagnose(kstep)
 
 
      WRITE(*,'(3x,a,a)') TRIM(resfile), ' created'
-     call flush(6)
+     CALL flush(6)
 
 
      !  Data group
@@ -222,16 +222,10 @@ SUBROUTINE diagnose(kstep)
      CALL attach(fidres, "/data/input","cpu_time",finish-start)
 
      ! Display computational time cost
-     IF     ( FLOOR((finish-start)/3600.) .GT. 0 ) THEN !display h min s
-       WRITE(*,*) 'CPU Time = ', FLOOR((finish-start)/3600.), '[h]', &
-       FLOOR((finish-start)/3600. - FLOOR((finish-start)/3600.))*60., '[min]',&
-       FLOOR((finish-start)/60.   - FLOOR((finish-start)/3600. - FLOOR((finish-start)/3600.))*60.)*60, '[s]'
-     ELSEIF ( FLOOR((finish-start)/60.)   .GT. 0 ) THEN !display min s
-       WRITE(*,*) 'CPU Time = ', FLOOR((finish-start)/60.), '[min]', &
-       FLOOR((finish-start)/60. - FLOOR(finish-start)/60.)*60., '[s]'
-     ELSE ! display s
-       WRITE(*,*) 'CPU Time = ', FLOOR((finish-start)), '[s]'
-     ENDIF
+     CALL display_h_min_s(finish-start)
+     ! WRITE(*,*) 'Time estimated :'
+     ! CALL display_h_min_s(time_est)
+
      !   Close all diagnostic files
      CALL closef(fidres)
   END IF
@@ -251,7 +245,7 @@ SUBROUTINE diagnose_0d
           '*** Timestep (this run/total) =', step, '/', cstep, 'Time =', time, 'dt =', dt
   WRITE(*,*)
 
-  ! flush stdout of all ranks. Usually ONLY rank 0 should write, but error messages might be written from other ranks as well
+  ! flush stdout of all ranks. Usually ONLY rank 0 should WRITE, but error messages might be written from other ranks as well
   CALL FLUSH(stdout)
 
 END SUBROUTINE diagnose_0d
diff --git a/src/diagnostics_par_mod.F90 b/src/diagnostics_par_mod.F90
index 6f065ecf..30a6424b 100644
--- a/src/diagnostics_par_mod.F90
+++ b/src/diagnostics_par_mod.F90
@@ -15,12 +15,12 @@ MODULE diagnostics_par
   INTEGER, PUBLIC :: nsave_cp = 1e4
 
   !  HDF5 file
-  CHARACTER(len=128), PUBLIC :: resfile0 = "results"   ! Head of main result file name
-  CHARACTER(len=128), PUBLIC :: resfile                ! Main result file
+  CHARACTER(len=256), PUBLIC :: resfile0 = "results"   ! Head of main result file name
+  CHARACTER(len=256), PUBLIC :: resfile                ! Main result file
   INTEGER, PUBLIC            :: job2load               ! jobnum of the checkpoint to load
   INTEGER, PUBLIC            :: fidres                 ! FID for resfile
-  CHARACTER(len=128), PUBLIC :: rstfile0 = "restart"   ! Head of restart file name
-  CHARACTER(len=128), PUBLIC :: rstfile                ! Full restart file
+  CHARACTER(len=256), PUBLIC :: rstfile0 = "restart"   ! Head of restart file name
+  CHARACTER(len=256), PUBLIC :: rstfile                ! Full restart file
   INTEGER, PUBLIC            :: fidrst                 ! FID for restart file
 
   PUBLIC :: output_par_readinputs, output_par_outputinputs
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index 48eaac9a..76179631 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -23,9 +23,9 @@ MODULE fourier
     COMPLEX(dp), DIMENSION(Nkr), INTENT(OUT):: Fk_out
     integer*8 plan
 
-    call dfftw_plan_dft_r2c_1d(plan,Nr,fx_in,Fk_out,FFTW_FORWARD,FFTW_ESTIMATE)
-    call dfftw_execute_dft_r2c(plan, fx_in, Fk_out)
-    call dfftw_destroy_plan(plan)
+    CALL dfftw_plan_dft_r2c_1d(plan,Nr,fx_in,Fk_out,FFTW_FORWARD,FFTW_ESTIMATE)
+    CALL dfftw_execute_dft_r2c(plan, fx_in, Fk_out)
+    CALL dfftw_destroy_plan(plan)
 
   END SUBROUTINE fft_r2cc
 
@@ -38,9 +38,9 @@ MODULE fourier
     REAL(dp),    DIMENSION(Nr),   INTENT(OUT):: fx_out
     integer*8 plan
 
-    call dfftw_plan_dft_c2r_1d(plan,Nr,Fk_in,fx_out,FFTW_BACKWARD,FFTW_ESTIMATE)
-    call dfftw_execute_dft_c2r(plan, Fk_in, fx_out)
-    call dfftw_destroy_plan(plan)
+    CALL dfftw_plan_dft_c2r_1d(plan,Nr,Fk_in,fx_out,FFTW_BACKWARD,FFTW_ESTIMATE)
+    CALL dfftw_execute_dft_c2r(plan, Fk_in, fx_out)
+    CALL dfftw_destroy_plan(plan)
     fx_out = fx_out/Nr
 
   END SUBROUTINE ifft_cc2r
@@ -55,9 +55,9 @@ MODULE fourier
       integer*8 plan
 
       !!! 2D Forward FFT ________________________!
-      call dfftw_plan_dft_r2c_2d(plan,Nr,Nz,ffx_in,FFk_out,FFTW_FORWARD,FFTW_ESTIMATE)
-      call dfftw_execute_dft_r2c(plan,ffx_in,FFk_out)
-      call dfftw_destroy_plan(plan)
+      CALL dfftw_plan_dft_r2c_2d(plan,Nr,Nz,ffx_in,FFk_out,FFTW_FORWARD,FFTW_ESTIMATE)
+      CALL dfftw_execute_dft_r2c(plan,ffx_in,FFk_out)
+      CALL dfftw_destroy_plan(plan)
 
   END SUBROUTINE fft2_r2cc
 
@@ -73,9 +73,9 @@ MODULE fourier
 
       tmp_c = FFk_in
       !!! 2D Backward FFT ________________________!
-      call dfftw_plan_dft_c2r_2d(plan,Nr,Nz,tmp_c,ffx_out,FFTW_BACKWARD,FFTW_ESTIMATE)
-      call dfftw_execute_dft_c2r(plan,tmp_c,ffx_out)
-      call dfftw_destroy_plan(plan)
+      CALL dfftw_plan_dft_c2r_2d(plan,Nr,Nz,tmp_c,ffx_out,FFTW_BACKWARD,FFTW_ESTIMATE)
+      CALL dfftw_execute_dft_c2r(plan,tmp_c,ffx_out)
+      CALL dfftw_destroy_plan(plan)
       ffx_out = ffx_out/Nr/Nz
 
   END SUBROUTINE ifft2_cc2r
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 125fa4a9..9c98a1a7 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -196,8 +196,8 @@ CONTAINS
     ! Put kz0RT to the nearest grid point on kz
     ikz0KH = NINT(kr0KH/deltakr)+1
     kr0KH  = kzarray(ikz0KH)
-    write(*,*) 'ikz0KH = ', ikz0KH
-    write(*,*) 'kr0KH = ', kr0KH
+    WRITE(*,*) 'ikz0KH = ', ikz0KH
+    WRITE(*,*) 'kr0KH = ', kr0KH
 
   END SUBROUTINE set_kzgrid
 
diff --git a/src/inital.F90 b/src/inital.F90
index b3d9fc9b..f1fa2efe 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -12,30 +12,34 @@ SUBROUTINE inital
 
   implicit none
 
+  real :: start_init, end_init, time_estimation
+
   CALL set_updatetlevel(1)
   !!!!!! Set the moments arrays Nepj, Nipj !!!!!!
-  write(*,*) 'Init moments'
+  ! WRITE(*,*) 'Init moments'
   IF ( RESTART ) THEN
     CALL load_cp
   ELSE
     CALL init_moments
   ENDIF
   !!!!!! Set phi !!!!!!
-  write(*,*) 'Init phi'
+  ! WRITE(*,*) 'Init phi'
   CALL poisson
 
   !!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!!
   IF ( NON_LIN .OR. (A0KH .NE. 0)) THEN;
-    write(*,*) 'Init Sapj'
+    ! WRITE(*,*) 'Init Sapj'
     CALL compute_Sapj
-    WRITE(*,*) 'Building Dnjs table'
+
+    ! WRITE(*,*) 'Building Dnjs table'
     CALL build_dnjs_table
   ENDIF
   !!!!!! Load the full coulomb collision operator coefficients !!!!!!
   IF (CO .EQ. -1) THEN
-    WRITE(*,*) '=== Load Full Coulomb matrix ==='
+    ! WRITE(*,*) '=== Load Full Coulomb matrix ==='
     CALL load_FC_mat
   ENDIF
+
 END SUBROUTINE inital
 !******************************************************************************!
 
@@ -69,10 +73,10 @@ SUBROUTINE init_moments
     END DO
 
   ELSE
-    sigma    = 5._dp     ! Gaussian sigma
-    gain     = 0.5_dp    ! Gaussian mean
-    kz_shift = 0.5_dp    ! Gaussian z shift
     !**** Gaussian initialization (Hakim 2017) *********************************
+    ! sigma    = 5._dp     ! Gaussian sigma
+    ! gain     = 0.5_dp    ! Gaussian mean
+    ! kz_shift = 0.5_dp    ! Gaussian z shift
     ! moments_i = 0; moments_e = 0;
     !   DO ikr=ikrs,ikre
     !     kr = krarray(ikr)
@@ -89,8 +93,8 @@ SUBROUTINE init_moments
         DO ikr=ikrs,ikre
           DO ikz=ikzs,ikze
             CALL RANDOM_NUMBER(noise)
-            moments_e( ip,ij,     ikr,    ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) &
-                                                    * AA_r(ikr) * AA_z(ikz)
+            moments_e( ip,ij,     ikr,    ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) !&
+                                                    ! * AA_r(ikr) * AA_z(ikz)
           END DO
         END DO
         DO ikz=2,Nkz/2 !symmetry at kr = 0
@@ -104,8 +108,8 @@ SUBROUTINE init_moments
         DO ikr=ikrs,ikre
           DO ikz=ikzs,ikze
             CALL RANDOM_NUMBER(noise)
-            moments_i( ip,ij,ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) &
-                                                    * AA_r(ikr) * AA_z(ikz)
+            moments_i( ip,ij,ikr,ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp) !&
+                                                    ! * AA_r(ikr) * AA_z(ikz)
           END DO
         END DO
         DO ikz=2,Nkz/2 !symmetry at kr = 0
@@ -115,26 +119,6 @@ SUBROUTINE init_moments
       END DO
     END DO
 
-    !**** Sinusoidal phi initialization for Kelvin-Helmholtz *******************
-  !   phi(1,FLOOR(0.5_dp/deltakr)) = 1._dp ! Trigger only mode kr = 1
-  !   ! moments_e( :,:, 1,FLOOR(0.5_dp/deltakr), :) = 1._dp
-  !   ! moments_i( :,:, 1,FLOOR(0.5_dp/deltakr), :) = 1._dp
-  !
-  !   DO ikr=ikrs,ikre
-  !     DO ikz=ikzs,ikze
-  !       CALL RANDOM_NUMBER(noise)
-  !       moments_e( :,:, ikr,ikz, :) = moments_e( :,:, ikr,ikz, :) &
-  !       + initnoise_moments*(noise-0.5_dp) ! adding noise
-  !
-  !       CALL RANDOM_NUMBER(noise)
-  !       moments_i( :,:, ikr,ikz, :) = moments_i( :,:, ikr,ikz, :) &
-  !       + initnoise_moments*(noise-0.5_dp) ! adding noise
-  !
-  !       CALL RANDOM_NUMBER(noise)
-  !       phi(ikr,ikz) = phi(ikr,ikz) + initnoise_moments*(noise-0.5_dp) ! adding noise
-  !     ENDDO
-  !   ENDDO
-  !   CALL moments_eq_rhs
   ENDIF
 END SUBROUTINE init_moments
 !******************************************************************************!
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index bdecea05..35cc155a 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -234,7 +234,7 @@ SUBROUTINE moments_eq_rhs
           ! Adding non linearity
           IF ( NON_LIN .OR. (A0KH .NE. 0) ) THEN
             moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = &
-              moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) - Sepj(ip,ij,ikr,ikz)*Nkr*Nkz
+              moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) - Sepj(ip,ij,ikr,ikz)
           ENDIF
 
         END DO kzloope
@@ -428,11 +428,11 @@ SUBROUTINE moments_eq_rhs
               -imagu * kpar*(TNapp1j + TNapm1j + xphijpar*kernelj*phi(ikr,ikz)) &
               - mu*kperp2**2 * moments_i(ip,ij,ikr,ikz,updatetlevel) &
                + TColl
-               
+
           ! Adding non linearity
           IF ( NON_LIN .OR. (A0KH .NE. 0) ) THEN
            moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = &
-             moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) - Sipj(ip,ij,ikr,ikz)*Nkr*Nkz
+             moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) - Sipj(ip,ij,ikr,ikz)
           ENDIF
 
         END DO kzloopi
diff --git a/src/poisson.F90 b/src/poisson.F90
index 760a9fe1..c0fdbcfb 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -80,7 +80,7 @@ SUBROUTINE poisson
       gammaD_phi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i
 
       IF ( (gammaD .EQ. 0) .AND. (abs(kr)+abs(kz) .NE. 0._dp) ) THEN
-        write(*,*) 'Warning gammaD = 0', sum_kernel2_i
+        WRITE(*,*) 'Warning gammaD = 0', sum_kernel2_i
       ENDIF
 
       phi(ikr, ikz) =  gammaD_phi/gammaD
diff --git a/src/prec_const_mod.F90 b/src/prec_const_mod.F90
index a7257cc0..7282015b 100644
--- a/src/prec_const_mod.F90
+++ b/src/prec_const_mod.F90
@@ -58,11 +58,11 @@ MODULE prec_const
       dp_r = range(b)
       dp_p = precision(b)
 
-      call mpi_comm_rank(MPI_COMM_WORLD,me,ierr)
+      CALL mpi_comm_rank(MPI_COMM_WORLD,me,ierr)
 
       !Create MPI datatypes that support the specific size
-      call MPI_Type_create_f90_real(sp_p,sp_r,MPI_sp,ierr)
-      call MPI_Type_create_f90_real(dp_p,dp_r,MPI_dp,ierr)
+      CALL MPI_Type_create_f90_real(sp_p,sp_r,MPI_sp,ierr)
+      CALL MPI_Type_create_f90_real(dp_p,dp_r,MPI_dp,ierr)
 
     END SUBROUTINE INIT_PREC_CONST
 
diff --git a/src/readinputs.F90 b/src/readinputs.F90
index 44bfcf8e..776f9009 100644
--- a/src/readinputs.F90
+++ b/src/readinputs.F90
@@ -16,17 +16,22 @@ SUBROUTINE readinputs
 
   ! Load grid data from input file
   CALL grid_readinputs
+  WRITE(*,*) 'check'
 
   ! Load diagnostic options from input file
   CALL output_par_readinputs
+  WRITE(*,*) 'check'
 
   ! Load model parameters from input file
   CALL model_readinputs
+  WRITE(*,*) 'check'
 
   ! Load initial condition parameters from input file
   CALL initial_readinputs
+  WRITE(*,*) 'check'
 
   ! Load parameters for time integration from input file
   CALL time_integration_readinputs
+  WRITE(*,*) 'check'
 
 END SUBROUTINE readinputs
diff --git a/src/stepon.F90 b/src/stepon.F90
index a900af8f..f23c5ff4 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -14,9 +14,7 @@ SUBROUTINE stepon
 
   INTEGER :: num_step
 
-
    DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
-
       ! Compute right hand side of moments hierarchy equation
       CALL moments_eq_rhs
       ! Advance from updatetlevel to updatetlevel+1 (according to num. scheme)
diff --git a/src/time_integration_mod.F90 b/src/time_integration_mod.F90
index 3f473707..03e7098b 100644
--- a/src/time_integration_mod.F90
+++ b/src/time_integration_mod.F90
@@ -71,7 +71,7 @@ CONTAINS
        WRITE(*,*) 'Cannot initialize time integration scheme. Name invalid.'
     END SELECT
 
-    write(*,*) " Time integration with ", numerical_scheme
+    WRITE(*,*) " Time integration with ", numerical_scheme
 
   END SUBROUTINE set_numerical_scheme
 
diff --git a/wk/fort.90 b/wk/fort.90
index ca350602..c849db79 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -1,7 +1,7 @@
 &BASIC
   nrun   = 100000000
   dt     = 0.01
-  tmax   = 150
+  tmax   = 200
   RESTART = .false.
 /
 &GRID
@@ -10,32 +10,32 @@
   pmaxi = 2
   jmaxi = 1
   Nr   = 256
-  Lr = 40
+  Lr = 66
   Nz   = 256
-  Lz = 40
+  Lz = 66
   kpar = 0
 /
 &OUTPUT_PAR
   nsave_0d = 0
   nsave_1d = 0
-  nsave_2d = 20
-  nsave_5d = 200
+  nsave_2d = 50
+  nsave_5d = 1000
   write_Na00    = .true.
   write_moments = .true.
   write_phi     = .true.
   write_non_lin     = .true.
   write_doubleprecision = .true.
-  resfile0      = 'ZP_forced_sym_256x128_L_40_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-03__mu_1e-04_'
-  rstfile0      = '../checkpoint/cp_ZP_forced_sym_256x128_L_40_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-03__mu_1e-04_'
-  job2load      = 0
+  resfile0      = '../results/ZP/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-04/outputs'
+  rstfile0      = '../results/ZP/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-04/checkpoint'
+  job2load      = 1
 /
 &MODEL_PAR
   ! Collisionality
-  CO      = 0
+  CO      = -2
   DK      = .false.
   NON_LIN = .true.
-  mu      = 0.0001
-  nu      = 0.001
+  mu      = 0.0005
+  nu      = 0.01
   tau_e   = 1
   tau_i   = 1
   sigma_e = 0.023338
@@ -52,7 +52,7 @@
 &INITIAL_CON
   only_Na00         =.false.
   initback_moments  =0
-  initnoise_moments =5e-05
+  initnoise_moments =1e-05
   iseed             =42
   selfmat_file      ='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10.h5'
   eimat_file        ='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_2_Jmaxe_1_Pmaxi_2_Jmaxi_1_pamaxx_10_tau_1.0000_mu_0.0233.h5'
diff --git a/wk/parameters_ZP.m b/wk/parameters_ZP.m
index 45a749ac..4033fbd2 100644
--- a/wk/parameters_ZP.m
+++ b/wk/parameters_ZP.m
@@ -5,40 +5,39 @@ default_plots_options
 %% Set Up parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 1e-3;   % Collision frequency
+NU      = 1e-2;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
 ETAB    = 0.5;    % Magnetic gradient
 ETAN    = 1.0;    % Density gradient
 ETAT    = 0.0;    % Temperature gradient
-MU      = 1e-4;   % Hyper diffusivity coefficient
-LAMBDAD = 0.0; 
-NOISE0  = 5.0e-5;
+MU      = 5e-4;   % Hyper diffusivity coefficient
+NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
 N       = 256;     % Frequency gridpoints (Nkr = N/2)
-L       = 40;     % Size of the squared frequency domain
-KREQ0   = 0;      % put kr = 0
-PMAXE   = 02;     % Highest electron Hermite polynomial degree
-JMAXE   = 01;     % Highest ''       Laguerre ''
-PMAXI   = 02;     % Highest ion      Hermite polynomial degree
-JMAXI   = 01;     % Highest ''       Laguerre ''
-KPAR    = 0.0;    % Parellel wave vector component
+L       = 66;     % Size of the squared frequency domain
+PMAXE   = 2;     % Highest electron Hermite polynomial degree
+JMAXE   = 1;     % Highest ''       Laguerre ''
+PMAXI   = 2;     % Highest ion      Hermite polynomial degree
+JMAXI   = 1;     % Highest ''       Laguerre ''
 %% TIME PARAMETERS 
-TMAX    = 150;  % Maximal time unit
-DT    = 1e-2;   % Time step
-SPS2D   = 5;      % Sampling per time unit for 2D arrays
-SPS5D   = 0.5;    % Sampling per time unit for 5D arrays
+TMAX    = 200;  % Maximal time unit
+DT      = 1e-2;   % Time step
+SPS2D   = 2;      % Sampling per time unit for 2D arrays
+SPS5D   = 0.1;    % Sampling per time unit for 5D arrays
 RESTART = 0;      % To restart from last checkpoint
-JOB2LOAD= 00;
+JOB2LOAD= 01;
 %% OPTIONS
-SIMID   = 'ZP_forced_sym';  % Name of the simulation
-NON_LIN = 1 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
-CO      = 0;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+SIMID   = 'ZP';  % Name of the simulation
+CO      = -2;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% unused
+%% unused
 KR0KH   = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
 NO_E    = 0;  % Remove electrons dynamic
 % DK    = 0;  % Drift kinetic model (put every kernel_n to 0 except n=0 to 1)
-
+KREQ0   = 0;      % put kr = 0
+KPAR    = 0.0;    % Parellel wave vector component
+LAMBDAD = 0.0; 
+NON_LIN = 1 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
 
 setup
-- 
GitLab