Trying to solve

Still I don't understand this basic thing...

src/modules/solver/electromagnetic/moduleEM.f90

@@ -170,8 +170,18 @@ MODULE moduleEM
       REAL(8), INTENT(in):: phi(1:n)
       REAL(8):: n_e(1:n)
       REAL(8):: n_e0 = 1.0D16, phi_0 = -520.0D0, T_e = 11604.0
+      INTEGER:: i
 
-      n_e = n_e0 / n_ref * EXP(-qe * (phi*Volt_ref - phi_0) / (kb * T_e))
+      DO i =1, n 
+        IF (phi(i)*Volt_ref >= phi_0) THEN
+          n_e(i) = n_e0 / n_ref * EXP(-qe * (phi(i)*Volt_ref - phi_0) / (kb * T_e))
+
+        ELSE
+          n_e(i) = 0.D0
+
+        END IF
+
+      END DO
 
       RETURN
 
@@ -185,7 +195,7 @@ MODULE moduleEM
       INTEGER, SAVE:: INFO
       INTEGER:: n
       REAL(8), ALLOCATABLE, SAVE:: tempF(:)
-      REAL(8), ALLOCATABLE:: n_e(:), phi_old(:)
+      REAL(8), ALLOCATABLE, SAVE:: n_e(:), phi_old(:)
       INTEGER:: k
       EXTERNAL:: dgetrs
 
@@ -193,21 +203,19 @@ MODULE moduleEM
       ALLOCATE(tempF(1:mesh%numNodes))
       ALLOCATE(n_e(1:mesh%numNodes))
       ALLOCATE(phi_old(1:mesh%numNodes))
-      DO n = 1, mesh%numNodes
-        tempF(n) = mesh%nodes(n)%obj%emData%phi
-
-      END DO
-      n_e = BoltzmannElectron(tempF, mesh%numNodes)
-      !$OMP END SINGLE
-
+      n_e = 0.D0
       CALL assembleSourceVector(tempF, n_e)
 
+      !$OMP END SINGLE
+
       !$OMP SINGLE
-      DO k = 1, 5
+      DO k = 1, 10
         phi_old = tempF
         CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, &
                                 mesh%IPIV, tempF,  mesh%numNodes, info)
 
+        PRINT *, MAXVAL(n_e),   MINVAL(n_e)
+        PRINT *, MAXVAL(tempF), MINVAL(tempF)
         PRINT*, k, "diff = ", MAXVAL(ABS(tempF - phi_old))
         n_e = BoltzmannElectron(tempF, mesh%numNodes)
         CALL assembleSourceVector(tempF, n_e)
@@ -232,7 +240,7 @@ MODULE moduleEM
       END IF
 
       !$OMP SINGLE
-      DEALLOCATE(tempF)
+      DEALLOCATE(tempF, n_e, phi_old)
       !$OMP END SINGLE
 
     END SUBROUTINE solveElecFieldBoltzmann