adiabatic ion model (not tested)

src/model_mod.F90

@@ -23,7 +23,8 @@ MODULE model
   REAL(xp), PUBLIC, PROTECTED :: lambdaD =  0._xp     ! Debye length
   REAL(xp), PUBLIC, PROTECTED ::    beta =  0._xp     ! electron plasma Beta (8piNT_e/B0^2)
   LOGICAL,  PUBLIC            :: ADIAB_E =  .false.   ! adiabatic electron model
-  REAL(xp), PUBLIC, PROTECTED ::   tau_e =  1.0        ! electron temperature ratio for adiabatic electrons
+  LOGICAL,  PUBLIC            :: ADIAB_I =  .false.   ! adiabatic ion model
+  REAL(xp), PUBLIC, PROTECTED ::   tau_i =  1.0       ! electron-ion temperature ratio for ion adiabatic model
   ! Auxiliary variable
   LOGICAL,  PUBLIC, PROTECTED ::      EM =  .false.   ! Electromagnetic effects flag
   LOGICAL,  PUBLIC, PROTECTED ::  MHD_PD =  .false.   ! MHD pressure drift
@@ -43,18 +44,23 @@ CONTAINS
 
     NAMELIST /MODEL_PAR/ KERN, LINEARITY, RM_LD_T_EQ, &
                          mu_x, mu_y, N_HD, HDz_h, mu_z, mu_p, mu_j, HYP_V, Na,&
-                         nu, k_gB, k_cB, lambdaD, MHD_PD, beta, ADIAB_E, tau_e
+                         nu, k_gB, k_cB, lambdaD, MHD_PD, beta, ADIAB_E, ADIAB_I, tau_i
 
     READ(lu_in,model_par)
 
-    IF((HYP_V .EQ. 'dvpar4') .AND. (num_procs_p .GT. 1)) THEN
+    IF (ADIAB_E .AND. ADIAB_I) &
+      ERROR STOP '>> ERROR << cannot have both adiab e and adiab i models'
+
+    IF((HYP_V .EQ. 'dvpar4') .AND. (num_procs_p .GT. 1)) &
       ERROR STOP '>> ERROR << dvpar4 velocity dissipation is not compatible with current p parallelization'
-    ENDIF
 
     IF(Na .EQ. 1) THEN
-      IF(.NOT. ADIAB_E) ERROR STOP "With one species, ADIAB_E must be set to .true. STOP"
-      CALL speak('Adiabatic electron model -> beta = 0')
-      beta = 0._xp
+      IF((.NOT. ADIAB_E) .AND. (.NOT. ADIAB_I)) ERROR STOP "With one species, ADIAB_E or ADIAB_I must be set to .true. STOP"
+      IF(ADIAB_E) THEN
+        CALL speak('Adiabatic electron model -> beta = 0')
+        beta = 0._xp
+      ENDIF
+      IF(ADIAB_I) CALL speak('Adiabatic ion model -> beta = 0')
     ENDIF
 
     IF(beta .GT. 0) THEN
@@ -91,7 +97,8 @@ CONTAINS
     CALL attach(fid, TRIM(str),    "MHD_PD",  MHD_PD)
     CALL attach(fid, TRIM(str),      "beta",    beta)
     CALL attach(fid, TRIM(str),   "ADIAB_E", ADIAB_E)
-    CALL attach(fid, TRIM(str),     "tau_e",   tau_e)
+    CALL attach(fid, TRIM(str),   "ADIAB_I", ADIAB_I)
+    CALL attach(fid, TRIM(str),     "tau_i",   tau_i)
   END SUBROUTINE model_outputinputs
 
 END MODULE model

src/numerics_mod.F90

@@ -134,10 +134,10 @@ SUBROUTINE evaluate_poisson_op
   USE grid,    ONLY : local_na, local_nkx, local_nky, local_nz,&
                       kxarray, kyarray, local_nj, ngj, ngz, ieven
   USE species, ONLY : q2_tau
-  USE model,   ONLY : ADIAB_E
+  USE model,   ONLY : ADIAB_E, ADIAB_I, tau_i
   USE prec_const, ONLY: xp
   IMPLICIT NONE
-  REAL(xp)    :: pol_tot, operator, operator_ion     ! (Z^2/tau (1-sum_n kernel_na^2))
+  REAL(xp)    :: pol_tot, operator_ion     ! (Z^2/tau (1-sum_n kernel_na^2))
   INTEGER     :: in,ikx,iky,iz,ia
   REAL(xp)    :: sumker     ! (Z_a^2/tau_a (1-sum_n kernel_na^2))
 
@@ -149,7 +149,7 @@ SUBROUTINE evaluate_poisson_op
   IF( (kxarray(ikx).EQ.0._xp) .AND. (kyarray(iky).EQ.0._xp) ) THEN
       inv_poisson_op(iky, ikx, iz) =  0._xp
       inv_pol_ion   (iky, ikx, iz) =  0._xp
-ELSE
+  ELSE
     ! loop over n only up to the max polynomial degree
     pol_tot = 0._xp  ! total polarisation term
     a:DO ia = 1,local_na ! sum over species
@@ -161,12 +161,15 @@ ELSE
       pol_tot = pol_tot + q2_tau(ia) - sumker
     ENDDO a
     operator_ion = pol_tot
-    IF(ADIAB_E) THEN ! Adiabatic electron model
+
+    IF(ADIAB_E) & ! Adiabatic electron model
       pol_tot = pol_tot + 1._xp
-    ENDIF
-    operator = pol_tot
+      
     inv_poisson_op(iky, ikx, iz) =  1._xp/pol_tot
     inv_pol_ion   (iky, ikx, iz) =  1._xp/operator_ion
+
+    IF(ADIAB_I) & ! adiabatic ions
+      inv_poisson_op(iky, ikx, iz)  = tau_i
   ENDIF
   END DO zloop
   END DO kyloop

src/solve_EM_fields.F90

@@ -22,7 +22,7 @@ CONTAINS
                                 ip0, total_nj, ngj
     USE calculus,         ONLY: simpson_rule_z
     USE parallel,         ONLY: manual_3D_bcast
-    USE model,            ONLY: lambdaD, ADIAB_E
+    USE model,            ONLY: lambdaD, ADIAB_E, ADIAB_I, tau_i
     use species,          ONLY: q
     USE processing,       ONLY: compute_density
     USE geometry,         ONLY: iInt_Jacobian, Jacobian
@@ -67,9 +67,6 @@ CONTAINS
               ENDIF
               rho = rho + fsa_phi
             ENDIF
-            !!!!!!!!!!!!!!! adiabatic ions ?
-            ! IF (ADIAB_I) THEN
-            ! ENDIF
             !!!!!!!!!!!!!!! Inverting the poisson equation
             DO iz = 1,local_nz
               phi(iky,ikx,iz+ngz/2) = inv_poisson_op(iky,ikx,iz)*rho(iz)

src/species_mod.F90

@@ -37,7 +37,7 @@ CONTAINS
   SUBROUTINE species_readinputs
     !    Read the input parameters
     USE basic, ONLY : lu_in
-    USE model, ONLY : Na, nu, ADIAB_E, MHD_PD
+    USE model, ONLY : Na, nu, ADIAB_E, ADIAB_I, MHD_PD
     USE prec_const
     IMPLICIT NONE
     INTEGER :: ia,ib
@@ -80,6 +80,8 @@ CONTAINS
       SELECT CASE (name_)
       CASE ('electrons','e','electron')
         ADIAB_E = .FALSE.
+      CASE ('ions','i','ion')
+        ADIAB_I = .FALSE.
       END SELECT
     ENDDO
     IF (.NOT. MHD_PD) Ptot = 0._dp