diff --git a/Python_analysis/Postprocessing/load_analyze_plot.py b/Python_analysis/Postprocessing/load_analyze_plot.py
index 6b7bc4cc..ce58d664 100644
--- a/Python_analysis/Postprocessing/load_analyze_plot.py
+++ b/Python_analysis/Postprocessing/load_analyze_plot.py
@@ -16,7 +16,7 @@
 
 #%% open Eric's matlab file for DIII-D magnetic field
 
-on_nersc=1
+on_nersc=0
 if on_nersc==0:
     filename='/home/21b/ExperimentalInputs/DIIID/200236/Matts_Bfield_data.mat'
 
@@ -42,7 +42,8 @@
 #run_directory=['/home/21b/KORC_RUNS/FROM_PERLMUTTER/TEST21']
 #run_directory=['/home/21b/KORC_RUNS/FROM_PERLMUTTER/TEST20b_rr']
 #run_directory=['/pscratch/sd/m/mbeidler/KORC_GPU_RUNS/DIIID_177031_GPU_TEST21a']
-run_directory=['/home/21b/KORC/test/elong_trans/rank_1']
+#run_directory=['/home/21b/KORC/test/elong_trans/rank_1']
+run_directory=['/home/21b/KORC/test/mars_EM/rank_1']
 
 for kk in range(0,dir_num):
 
@@ -341,16 +342,13 @@
 
 #%% open KORC field file
 
-useextfield=0
+useextfield=1
 
 if useextfield==1:
 
     #filename=r"/home/21b/KORC_RUNS/LOCAL/TEST2/AORSA_D3D_171089_200MHz_EFITgrid.h5"
-    
-    if on_nersc==0:
-        filename=r"../"+run_directory[0]+"/DIII_200236_MARS_CaseA.h5"
-    else:
-        filename=r"../../"+run_directory[0]+"/DIII_200236_MARS_CaseA.h5"
+    #filename=r"../../"+run_directory[0]+"/DIII_200236_MARS_CaseA.h5"
+    filename=r"/home/21b/KORC/test/mars_EM/ARC_42533M20_MARS.h5"
     
     with h5py.File(filename,'r') as f:
         NR=int(f['NR'][0])           
@@ -388,7 +386,7 @@
     interp_br = RegularGridInterpolator((Rg,Zg), np.transpose(BR_g[:,:]),method='cubic')
     interp_bz = RegularGridInterpolator((Rg,Zg), np.transpose(BZ_g[:,:]),method='cubic')
     
-    Rq=1.7
+    Rq=4.
     Zq=0.
     
     psip_q=interp_psip([Rq,Zq])
@@ -475,19 +473,23 @@
     vperp=vmag*np.sin(np.radians(eta))
     vpll=vmag*np.cos(np.radians(eta))
 
-rm=np.sqrt((R-R0)**2+((zz-Z0)/kappa)**2)
+try:
+    rm=np.sqrt((R-R0)**2+((zz-Z0)/kappa)**2)
+except: 
+    rm=np.sqrt((R-R0)**2+(zz-Z0)**2)
 
-flagDecon=np.zeros(np.shape(flagCon))
-flagDecon[(flagCon<1) & (flagRE>0)]=1
-confined=np.sum(flagCon,axis=1)
-deconfined=np.sum(flagDecon,axis=1)
-
-flagTherm=np.zeros(np.shape(flagCol))
-flagTherm[flagCol<1 & (flagRE>0)]=1
-Thermal=np.sum(flagTherm,axis=1)
-Energetic=np.sum(flagCol,axis=1)
 
 if 'flagRE' in outputs_list:
+  flagDecon=np.zeros(np.shape(flagCon))
+  flagDecon[(flagCon<1) & (flagRE>0)]=1
+  confined=np.sum(flagCon,axis=1)
+  deconfined=np.sum(flagDecon,axis=1)
+
+  flagTherm=np.zeros(np.shape(flagCol))
+  flagTherm[flagCol<1 & (flagRE>0)]=1
+  Thermal=np.sum(flagTherm,axis=1)
+  Energetic=np.sum(flagCol,axis=1) 
+    
   Total=np.sum(flagRE,axis=1)
     
   flagActive=flagRE*flagCon*flagCol
@@ -503,16 +505,15 @@
   flagSecondary[:,flagRE[-1,:]<1]=0
   Secondary=np.sum(flagSecondary,axis=1)
 
-Ipart=qe*vpll*bPHI/(2*np.pi*R*bmag)
-Ipart[flagActive==0]=0
-Itot=np.sum(Ipart,axis=1)
-Ipri=np.sum(Ipart*flagPrimary,axis=1)
-Isec=np.sum(Ipart*flagSecondary,axis=1)
-
-IPOLpart=qe*vpll*bpol/(2*np.pi*R*bmag)
-IPOLpart[flagActive==0]=0
-IPOL=np.sum(IPOLpart,axis=1)
-
+  Ipart=qe*vpll*bPHI/(2*np.pi*R*bmag)
+  Ipart[flagActive==0]=0
+  Itot=np.sum(Ipart,axis=1)
+  Ipri=np.sum(Ipart*flagPrimary,axis=1)
+  Isec=np.sum(Ipart*flagSecondary,axis=1)
+    
+  IPOLpart=qe*vpll*bpol/(2*np.pi*R*bmag)
+  IPOLpart[flagActive==0]=0
+  IPOL=np.sum(IPOLpart,axis=1)
 
 #KE=me*c**2/qe*np.sum((g-1)*flagActive)
 #KEtot=me*c**2/qe*np.sum((1-1)*flagRE)
@@ -1055,7 +1056,7 @@ def velocity_magnetic_to_collapsed_cartesian(v, eta, chi, Bx, By, Bz, phi_raw, p
 plt.rc('legend', fontsize=SMALL_SIZE)
 plt.rc('figure', titlesize=SMALL_SIZE)
 
-plot_histrm=1
+plot_histrm=0
 plot_LAC_ParamScaling=0
 plot_LAC_Escaling=0
 plot_GPUscaling=0
@@ -1074,7 +1075,7 @@ def velocity_magnetic_to_collapsed_cartesian(v, eta, chi, Bx, By, Bz, phi_raw, p
 plot_bz=0
 plot_bphi=0
 plot_psip=0
-plot_histRZ_ext=0
+plot_histRZ_ext=1
 plot_histKeta=0
 plot_histK=0
 plot_histeta=0
@@ -1097,7 +1098,7 @@ def velocity_magnetic_to_collapsed_cartesian(v, eta, chi, Bx, By, Bz, phi_raw, p
 plot_evoI = 0
 plot_evoRE = 0
 
-timeind_p=100
+timeind_p=0
 timeind_g=0
 
 #tloss=time[12]
diff --git a/build_scripts/build_nvhpc_gpu.sh b/build_scripts/build_nvhpc_gpu.sh
index 01939c17..1c8ca4e4 100755
--- a/build_scripts/build_nvhpc_gpu.sh
+++ b/build_scripts/build_nvhpc_gpu.sh
@@ -18,7 +18,7 @@ cmake \
   -DCMAKE_Fortran_FLAGS="-acc=gpu -gpu=deepcopy,cc89 -Mfree -fPIC -c++libs -Mpreprocess -DDOUBLE_PRECISION" \
   -DCMAKE_C_FLAGS="-mp -DDOUBLE_PRECISION" \
   -DCMAKE_CXX_FLAGS="-std=c++11 -mp -DDOUBLE_PRECISION" \
-  -DCMAKE_Fortran_FLAGS_DEBUG='-g -gpu=debug -O1 -Minfo=accel' \
+  -DCMAKE_Fortran_FLAGS_DEBUG='-g -gpu=debug -O1' \
   ..
 
 make -j VERBOSE=1
diff --git a/src/korc_fields.f90 b/src/korc_fields.f90
index d626d0ed..f518fa4e 100755
--- a/src/korc_fields.f90
+++ b/src/korc_fields.f90
@@ -384,9 +384,9 @@ subroutine analytical_fields_p_ACC(T_R,T_T,T_Z, &
     Br = Br + sqrt(dBr_norm_squared)
   end if
 
-  B_X = Bzeta*cZ - Bp*sT*sZ + Br*cT*sZ
-  B_Y = -Bzeta*sZ - Bp*sT*cZ + Br*cT*cZ
-  B_Z = Bp*cT + Br*sT
+  B_X = Bzeta*cZ - Bp*sT*sZ/kappa + Br*cT*sZ
+  B_Y = -Bzeta*sZ - Bp*sT*cZ/kappa + Br*cT*cZ
+  B_Z = Bp*cT + Br*sT/kappa
 
   Ezeta = -E0/( 1.0_rp + eta*cT)
 
@@ -2210,6 +2210,9 @@ subroutine load_field_data_from_hdf5(params,F)
         if (F%dims(2).gt.0) then
           dset = '/GR'
           call load_array_from_hdf5(h5file_id,dset,F%GR)
+
+          dset = '/FR'
+          call load_array_from_hdf5(h5file_id,dset,F%FR)
         endif
 
       end if
diff --git a/src/korc_interp.f90 b/src/korc_interp.f90
index ea70ef56..ea4d4f99 100755
--- a/src/korc_interp.f90
+++ b/src/korc_interp.f90
@@ -486,7 +486,7 @@ subroutine initialize_fields_interpolant(params,F)
 
     if(F%B1field) then
 
-      if (params%field_model(10:13).eq.'MARS') then
+      if (params%field_model(10:16).eq.'MARS   ') then
 
         write(output_unit_write,*) '2D n=1 MARS magnetic fields' 
 
@@ -1707,7 +1707,7 @@ subroutine initialize_fields_interpolant(params,F)
 
        if(F%E1field) then
 
-          if (params%field_model(10:13).eq.'MARS_EM') then
+          if (params%field_model(10:16).eq.'MARS_EM') then
 
             write(output_unit_write,*) '2D MARS NL EM electric fields'
             flush(output_unit_write)
@@ -3530,12 +3530,20 @@ subroutine provide_ezspline_marsEM_ACC(bfield_2d_local, &
   b1Imfield_2d_local_2,b1Refield_2d_local_3,b1Imfield_2d_local_3, &
   b1Refield_2d_local_4,b1Imfield_2d_local_4,b1Refield_2d_local_5, &
   b1Imfield_2d_local_5,b1Refield_2d_local_6,b1Imfield_2d_local_6, &
+  e1Refield_2d_local_1,e1Imfield_2d_local_1,e1Refield_2d_local_2, &
+  e1Imfield_2d_local_2,e1Refield_2d_local_3,e1Imfield_2d_local_3, &
+  e1Refield_2d_local_4,e1Imfield_2d_local_4,e1Refield_2d_local_5, &
+  e1Imfield_2d_local_5,e1Refield_2d_local_6,e1Imfield_2d_local_6, &
   fields_domain_local)
   TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(OUT)      :: bfield_2d_local
   TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(OUT)      :: b1Refield_2d_local_1,b1Refield_2d_local_2,b1Refield_2d_local_3
   TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(OUT)      :: b1Imfield_2d_local_1,b1Imfield_2d_local_2,b1Imfield_2d_local_3
   TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(OUT)      :: b1Refield_2d_local_4,b1Refield_2d_local_5,b1Refield_2d_local_6
   TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(OUT)      :: b1Imfield_2d_local_4,b1Imfield_2d_local_5,b1Imfield_2d_local_6
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(OUT)      :: e1Refield_2d_local_1,e1Refield_2d_local_2,e1Refield_2d_local_3
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(OUT)      :: e1Imfield_2d_local_1,e1Imfield_2d_local_2,e1Imfield_2d_local_3
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(OUT)      :: e1Refield_2d_local_4,e1Refield_2d_local_5,e1Refield_2d_local_6
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(OUT)      :: e1Imfield_2d_local_4,e1Imfield_2d_local_5,e1Imfield_2d_local_6
   TYPE(KORC_INTERPOLANT_DOMAIN),INTENT(OUT)        :: fields_domain_local
 
   bfield_2d_local=bfield_2d  
@@ -3551,6 +3559,18 @@ subroutine provide_ezspline_marsEM_ACC(bfield_2d_local, &
   b1Imfield_2d_local_5=b1Imfield_2d_5
   b1Refield_2d_local_6=b1Refield_2d_6
   b1Imfield_2d_local_6=b1Imfield_2d_6
+  e1Refield_2d_local_1=e1Refield_2d_1
+  e1Imfield_2d_local_1=e1Imfield_2d_1
+  e1Refield_2d_local_2=e1Refield_2d_2
+  e1Imfield_2d_local_2=e1Imfield_2d_2
+  e1Refield_2d_local_3=e1Refield_2d_3
+  e1Imfield_2d_local_3=e1Imfield_2d_3
+  e1Refield_2d_local_4=e1Refield_2d_4
+  e1Imfield_2d_local_4=e1Imfield_2d_4
+  e1Refield_2d_local_5=e1Refield_2d_5
+  e1Imfield_2d_local_5=e1Imfield_2d_5
+  e1Refield_2d_local_6=e1Refield_2d_6
+  e1Imfield_2d_local_6=e1Imfield_2d_6
   fields_domain_local=fields_domain
 
 end subroutine provide_ezspline_marsEM_ACC
@@ -3830,17 +3850,17 @@ subroutine interp_FOfields_marsEM_p(pchunk,F,Y_R,Y_PHI,Y_Z, &
     cPshift=cos(F%X%PHI(1)*Y_PHI(cc)-phase-fr(1)*time)
     sPshift=sin(F%X%PHI(1)*Y_PHI(cc)-phase-fr(1)*time)
 
-    B1_R = gr(1)*MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(B1Re_R*cPshift-B1Im_R*sPshift)
-    B1_PHI = gr(1)*MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
-    B1_Z = gr(1)*MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(B1Re_Z*cPshift-B1Im_Z*sPshift)
+    B1_R = gr(1)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(B1Re_R*cPshift-B1Im_R*sPshift)
+    B1_PHI = gr(1)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
+    B1_Z = gr(1)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(B1Re_Z*cPshift-B1Im_Z*sPshift)
 
     B_R = B0_R+B1_R
     B_PHI = B0_PHI+B1_PHI
     B_Z(cc) = B0_Z+B1_Z
 
-    E1_R = gr(1)*MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(E1Re_R*cPshift-E1Im_R*sPshift)
-    E1_PHI = gr(1)*MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
-    E1_Z = gr(1)*MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(E1Re_Z*cPshift-E1Im_Z*sPshift)
+    E1_R = gr(1)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(E1Re_R*cPshift-E1Im_R*sPshift)
+    E1_PHI = gr(1)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
+    E1_Z = gr(1)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(E1Re_Z*cPshift-E1Im_Z*sPshift)
 
     E_R = E1_R
     E_PHI = E1_PHI
@@ -3859,17 +3879,17 @@ subroutine interp_FOfields_marsEM_p(pchunk,F,Y_R,Y_PHI,Y_Z, &
     cPshift=cos(F%X%PHI(2)*Y_PHI(cc)-phase-fr(2)*time)
     sPshift=sin(F%X%PHI(2)*Y_PHI(cc)-phase-fr(2)*time)
 
-    B1_R = MARS_max/(1+MARS_max*exp(-time*gr(2)))*amp(2)*(B1Re_R*cPshift-B1Im_R*sPshift)
-    B1_PHI = MARS_max/(1+MARS_max*exp(-time*gr(2)))*amp(2)*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
-    B1_Z = MARS_max/(1+MARS_max*exp(-time*gr(2)))*amp(2)*(B1Re_Z*cPshift-B1Im_Z*sPshift)
+    B1_R = gr(2)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(2)))*(B1Re_R*cPshift-B1Im_R*sPshift)
+    B1_PHI = gr(2)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(2)))*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
+    B1_Z = gr(2)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(2)))*(B1Re_Z*cPshift-B1Im_Z*sPshift)
 
     B_R = B_R+B1_R
     B_PHI = B_PHI+B1_PHI
     B_Z(cc) = B_Z(cc)+B1_Z
 
-    E1_R = MARS_max/(1+MARS_max*exp(-time*gr(2)))*amp(2)*(E1Re_R*cPshift-E1Im_R*sPshift)
-    E1_PHI = MARS_max/(1+MARS_max*exp(-time*gr(2)))*amp(2)*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
-    E1_Z = MARS_max/(1+MARS_max*exp(-time*gr(2)))*amp(2)*(E1Re_Z*cPshift-E1Im_Z*sPshift)
+    E1_R = gr(2)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(2)))*(E1Re_R*cPshift-E1Im_R*sPshift)
+    E1_PHI = gr(2)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(2)))*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
+    E1_Z = gr(2)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(2)))*(E1Re_Z*cPshift-E1Im_Z*sPshift)
 
     E_R = E_R+E1_R
     E_PHI = E_PHI+E1_PHI
@@ -3888,17 +3908,17 @@ subroutine interp_FOfields_marsEM_p(pchunk,F,Y_R,Y_PHI,Y_Z, &
     cPshift=cos(F%X%PHI(3)*Y_PHI(cc)-phase-fr(3)*time)
     sPshift=sin(F%X%PHI(3)*Y_PHI(cc)-phase-fr(3)*time)
 
-    B1_R = MARS_max/(1+MARS_max*exp(-time*gr(3)))*amp(3)*(B1Re_R*cPshift-B1Im_R*sPshift)
-    B1_PHI = MARS_max/(1+MARS_max*exp(-time*gr(3)))*amp(3)*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
-    B1_Z = MARS_max/(1+MARS_max*exp(-time*gr(3)))*amp(3)*(B1Re_Z*cPshift-B1Im_Z*sPshift)
+    B1_R = gr(3)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(3)))*(B1Re_R*cPshift-B1Im_R*sPshift)
+    B1_PHI = gr(3)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(3)))*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
+    B1_Z = gr(3)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(3)))*(B1Re_Z*cPshift-B1Im_Z*sPshift)
 
     B_R = B_R+B1_R
     B_PHI = B_PHI+B1_PHI
     B_Z(cc) = B_Z(cc)+B1_Z
 
-    E1_R = MARS_max/(1+MARS_max*exp(-time*gr(3)))*amp(3)*(E1Re_R*cPshift-E1Im_R*sPshift)
-    E1_PHI = MARS_max/(1+MARS_max*exp(-time*gr(3)))*amp(3)*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
-    E1_Z = MARS_max/(1+MARS_max*exp(-time*gr(3)))*amp(3)*(E1Re_Z*cPshift-E1Im_Z*sPshift)
+    E1_R = gr(3)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(3)))*(E1Re_R*cPshift-E1Im_R*sPshift)
+    E1_PHI = gr(3)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(3)))*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
+    E1_Z = gr(3)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(3)))*(E1Re_Z*cPshift-E1Im_Z*sPshift)
 
     E_R = E_R+E1_R
     E_PHI = E_PHI+E1_PHI
@@ -3917,17 +3937,17 @@ subroutine interp_FOfields_marsEM_p(pchunk,F,Y_R,Y_PHI,Y_Z, &
     cPshift=cos(F%X%PHI(4)*Y_PHI(cc)-phase-fr(4)*time)
     sPshift=sin(F%X%PHI(4)*Y_PHI(cc)-phase-fr(4)*time)
 
-    B1_R = MARS_max/(1+MARS_max*exp(-time*gr(4)))*amp(4)*(B1Re_R*cPshift-B1Im_R*sPshift)
-    B1_PHI = MARS_max/(1+MARS_max*exp(-time*gr(4)))*amp(4)*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
-    B1_Z = MARS_max/(1+MARS_max*exp(-time*gr(4)))*amp(4)*(B1Re_Z*cPshift-B1Im_Z*sPshift)
+    B1_R = gr(4)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(4)))*(B1Re_R*cPshift-B1Im_R*sPshift)
+    B1_PHI = gr(4)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(4)))*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
+    B1_Z = gr(4)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(4)))*(B1Re_Z*cPshift-B1Im_Z*sPshift)
 
     B_R = B_R+B1_R
     B_PHI = B_PHI+B1_PHI
     B_Z(cc) = B_Z(cc)+B1_Z
 
-    E1_R = MARS_max/(1+MARS_max*exp(-time*gr(4)))*amp(4)*(E1Re_R*cPshift-E1Im_R*sPshift)
-    E1_PHI = MARS_max/(1+MARS_max*exp(-time*gr(4)))*amp(4)*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
-    E1_Z = MARS_max/(1+MARS_max*exp(-time*gr(4)))*amp(4)*(E1Re_Z*cPshift-E1Im_Z*sPshift)
+    E1_R = gr(4)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(4)))*(E1Re_R*cPshift-E1Im_R*sPshift)
+    E1_PHI = gr(4)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(4)))*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
+    E1_Z = gr(4)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(4)))*(E1Re_Z*cPshift-E1Im_Z*sPshift)
 
     E_R = E_R+E1_R
     E_PHI = E_PHI+E1_PHI
@@ -3946,17 +3966,17 @@ subroutine interp_FOfields_marsEM_p(pchunk,F,Y_R,Y_PHI,Y_Z, &
     cPshift=cos(F%X%PHI(5)*Y_PHI(cc)-phase-fr(5)*time)
     sPshift=sin(F%X%PHI(5)*Y_PHI(cc)-phase-fr(5)*time)
 
-    B1_R = MARS_max/(1+MARS_max*exp(-time*gr(5)))*amp(5)*(B1Re_R*cPshift-B1Im_R*sPshift)
-    B1_PHI = MARS_max/(1+MARS_max*exp(-time*gr(5)))*amp(5)*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
-    B1_Z = MARS_max/(1+MARS_max*exp(-time*gr(5)))*amp(5)*(B1Re_Z*cPshift-B1Im_Z*sPshift)
+    B1_R = gr(5)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(5)))*(B1Re_R*cPshift-B1Im_R*sPshift)
+    B1_PHI = gr(5)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(5)))*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
+    B1_Z = gr(5)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(5)))*(B1Re_Z*cPshift-B1Im_Z*sPshift)
 
     B_R = B_R+B1_R
     B_PHI = B_PHI+B1_PHI
     B_Z(cc) = B_Z(cc)+B1_Z
 
-    E1_R = MARS_max/(1+MARS_max*exp(-time*gr(5)))*amp(5)*(E1Re_R*cPshift-E1Im_R*sPshift)
-    E1_PHI = MARS_max/(1+MARS_max*exp(-time*gr(5)))*amp(5)*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
-    E1_Z = MARS_max/(1+MARS_max*exp(-time*gr(5)))*amp(5)*(E1Re_Z*cPshift-E1Im_Z*sPshift)
+    E1_R = gr(5)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(5)))*(E1Re_R*cPshift-E1Im_R*sPshift)
+    E1_PHI = gr(5)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(5)))*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
+    E1_Z = gr(5)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(5)))*(E1Re_Z*cPshift-E1Im_Z*sPshift)
 
     E_R = E_R+E1_R
     E_PHI = E_PHI+E1_PHI
@@ -3975,17 +3995,17 @@ subroutine interp_FOfields_marsEM_p(pchunk,F,Y_R,Y_PHI,Y_Z, &
     cPshift=cos(F%X%PHI(6)*Y_PHI(cc)-phase-fr(6)*time)
     sPshift=sin(F%X%PHI(6)*Y_PHI(cc)-phase-fr(6)*time)
 
-    B1_R = MARS_max/(1+MARS_max*exp(-time*gr(6)))*amp(6)*(B1Re_R*cPshift-B1Im_R*sPshift)
-    B1_PHI = MARS_max/(1+MARS_max*exp(-time*gr(6)))*amp(6)*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
-    B1_Z = MARS_max/(1+MARS_max*exp(-time*gr(6)))*amp(6)*(B1Re_Z*cPshift-B1Im_Z*sPshift)
+    B1_R = gr(6)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(6)))*(B1Re_R*cPshift-B1Im_R*sPshift)
+    B1_PHI = gr(6)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(6)))*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
+    B1_Z = gr(6)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(6)))*(B1Re_Z*cPshift-B1Im_Z*sPshift)
 
     B_R = B_R+B1_R
     B_PHI = B_PHI+B1_PHI
     B_Z(cc) = B_Z(cc)+B1_Z
 
-    E1_R = MARS_max/(1+MARS_max*exp(-time*gr(6)))*amp(6)*(E1Re_R*cPshift-E1Im_R*sPshift)
-    E1_PHI = MARS_max/(1+MARS_max*exp(-time*gr(6)))*amp(6)*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
-    E1_Z = MARS_max/(1+MARS_max*exp(-time*gr(6)))*amp(6)*(E1Re_Z*cPshift-E1Im_Z*sPshift)
+    E1_R = gr(6)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(6)))*(E1Re_R*cPshift-E1Im_R*sPshift)
+    E1_PHI = gr(6)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(6)))*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
+    E1_Z = gr(6)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(6)))*(E1Re_Z*cPshift-E1Im_Z*sPshift)
 
     E_R = E_R+E1_R
     E_PHI = E_PHI+E1_PHI
@@ -4163,9 +4183,21 @@ end subroutine interp_FOfields_marsNL_p_ACC
 
 subroutine interp_FOfields_marsEM_p_ACC(time,bfield_2d_local, &
   b1Refield_2d_local_1,b1Imfield_2d_local_1, &
-  psip_conv,amp,gr,nmode,phase,MARS_max,Bo,Ro,Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z,PSIp)
+  b1Refield_2d_local_2,b1Imfield_2d_local_2, &
+  b1Refield_2d_local_3,b1Imfield_2d_local_3, &
+  b1Refield_2d_local_4,b1Imfield_2d_local_4, &
+  b1Refield_2d_local_5,b1Imfield_2d_local_5, &
+  b1Refield_2d_local_6,b1Imfield_2d_local_6, &
+  e1Refield_2d_local_1,e1Imfield_2d_local_1, &
+  e1Refield_2d_local_2,e1Imfield_2d_local_2, &
+  e1Refield_2d_local_3,e1Imfield_2d_local_3, &
+  e1Refield_2d_local_4,e1Imfield_2d_local_4, &
+  e1Refield_2d_local_5,e1Imfield_2d_local_5, &
+  e1Refield_2d_local_6,e1Imfield_2d_local_6, &
+  psip_conv,amp,gr,fr,nmode,phase,MARS_max,MARS_quas_fac,Bo,Ro,Y_R,Y_PHI,Y_Z, &
+  B_X,B_Y,B_Z,E_X,E_Y,E_Z,PSIp)
   !$acc routine seq
-  REAL(rp),INTENT(IN)  :: time,MARS_max
+  REAL(rp),INTENT(IN)  :: time,MARS_max,MARS_quas_fac
   REAL(rp),INTENT(IN)   :: Y_R,Y_PHI,Y_Z
   REAL(rp),INTENT(OUT)   :: B_X,B_Y,B_Z
   REAL(rp)   :: B_R,B_PHI
@@ -4173,6 +4205,11 @@ subroutine interp_FOfields_marsEM_p_ACC(time,bfield_2d_local, &
   REAL(rp)   :: B1_R,B1_PHI,B1_Z
   REAL(rp)   :: B1Re_R,B1Re_PHI,B1Re_Z
   REAL(rp)   :: B1Im_R,B1Im_PHI,B1Im_Z
+  REAL(rp),INTENT(OUT)   :: E_X,E_Y,E_Z
+  REAL(rp)   :: E_R,E_PHI
+  REAL(rp)   :: E1_R,E1_PHI,E1_Z
+  REAL(rp)   :: E1Re_R,E1Re_PHI,E1Re_Z
+  REAL(rp)   :: E1Im_R,E1Im_PHI,E1Im_Z
   REAL(rp),INTENT(OUT)   :: PSIp
   REAL(rp)   :: cP,sP,cPshift,sPshift
   REAL(rp), DIMENSION(3)  :: A
@@ -4180,87 +4217,221 @@ subroutine interp_FOfields_marsEM_p_ACC(time,bfield_2d_local, &
   INTEGER                                      :: cc
   !! Particle chunk iterator.
   REAL(rp), INTENT(IN) :: psip_conv,phase,Bo,Ro
-  REAL(rp),DIMENSION(3),INTENT(IN) :: amp,gr,nmode
+  REAL(rp),DIMENSION(6),INTENT(IN) :: amp,gr,fr,nmode
   TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: bfield_2d_local
-  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: b1Refield_2d_local_1!,b1Refield_2d_local_2,b1Refield_2d_local_3
-  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: b1Imfield_2d_local_1!,b1Imfield_2d_local_2,b1Imfield_2d_local_3
-
-  !$acc routine (EZspline_interp2_FOmars) seq
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: b1Refield_2d_local_1,b1Refield_2d_local_2
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: b1Refield_2d_local_3,b1Refield_2d_local_4
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: b1Refield_2d_local_5,b1Refield_2d_local_6
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: b1Imfield_2d_local_1,b1Imfield_2d_local_2
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: b1Imfield_2d_local_3,b1Imfield_2d_local_4
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: b1Imfield_2d_local_5,b1Imfield_2d_local_6
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: e1Refield_2d_local_1,e1Refield_2d_local_2
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: e1Refield_2d_local_3,e1Refield_2d_local_4
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: e1Refield_2d_local_5,e1Refield_2d_local_6
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: e1Imfield_2d_local_1,e1Imfield_2d_local_2
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: e1Imfield_2d_local_3,e1Imfield_2d_local_4
+  TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN)      :: e1Imfield_2d_local_5,e1Imfield_2d_local_6
+
+  !$acc routine (EZspline_interp2_FOmarsEM) seq
   !$acc routine (EZspline_error) seq
      
   !write(6,*) Y_R
   !write(6,*) Y_Z
 
-  call EZspline_interp2_FOmars(bfield_2d_local%A,b1Refield_2d_local_1%R,b1Refield_2d_local_1%PHI, &
-    b1Refield_2d_local_1%Z,b1Imfield_2d_local_1%R,b1Imfield_2d_local_1%PHI,b1Imfield_2d_local_1%Z, &
-    Y_R,Y_Z,A,B1Re_R,B1Re_PHI,B1Re_Z,B1Im_R,B1Im_PHI,B1Im_Z,ezerr_local)
+  call EZspline_interp2_FOmarsEM(bfield_2d_local%A, &
+    b1Refield_2d_local_1%R,b1Refield_2d_local_1%PHI,b1Refield_2d_local_1%Z, &
+    b1Imfield_2d_local_1%R,b1Imfield_2d_local_1%PHI,b1Imfield_2d_local_1%Z, &
+    e1Refield_2d_local_1%R,e1Refield_2d_local_1%PHI,e1Refield_2d_local_1%Z, &
+    e1Imfield_2d_local_1%R,e1Imfield_2d_local_1%PHI,e1Imfield_2d_local_1%Z, &
+    Y_R,Y_Z,A, &
+    B1Re_R,B1Re_PHI,B1Re_Z,B1Im_R,B1Im_PHI,B1Im_Z, &
+    E1Re_R,E1Re_PHI,E1Re_Z,E1Im_R,E1Im_PHI,E1Im_Z, &
+    ezerr_local)
   call EZspline_error(ezerr_local)
 
-  !write(6,*) B1Re_R
-  !write(6,*) B1Im_R
-  !write(6,*) B1Re_PHI
-  !write(6,*) B1Im_PHI
-  !write(6,*) B1Re_Z
-  !write(6,*) B1Im_Z
-
   PSIp=A(1)
 
   B0_R = psip_conv*A(3)/Y_R
   B0_PHI = -Bo*Ro/Y_R
   B0_Z = -psip_conv*A(2)/Y_R
 
-  !write(6,*) A(1)
-  !write(6,*) A(2)
-  !write(6,*) A(3)
-
   cP=cos(Y_PHI)
   sP=sin(Y_PHI)
-  cPshift=cos(nmode(1)*Y_PHI-phase)
-  sPshift=sin(nmode(1)*Y_PHI-phase)
+  cPshift=cos(nmode(1)*Y_PHI-phase-fr(1)*time)
+  sPshift=sin(nmode(1)*Y_PHI-phase-fr(1)*time)
 
-  B1_R = MARS_max/(1+MARS_max*exp(-time*gr(1)))*amp(1)*(B1Re_R*cPshift-B1Im_R*sPshift)
-  B1_PHI = MARS_max/(1+MARS_max*exp(-time*gr(1)))*amp(1)*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
-  B1_Z = MARS_max/(1+MARS_max*exp(-time*gr(1)))*amp(1)*(B1Re_Z*cPshift-B1Im_Z*sPshift)
+  B1_R = gr(1)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(B1Re_R*cPshift-B1Im_R*sPshift)
+  B1_PHI = gr(1)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
+  B1_Z = gr(1)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(B1Re_Z*cPshift-B1Im_Z*sPshift)
 
   B_R = B0_R+B1_R
   B_PHI = B0_PHI+B1_PHI
   B_Z = B0_Z+B1_Z
 
-  !call EZspline_interp2_FOmars(bfield_2d_local%A,b1Refield_2d_local_2%R,b1Refield_2d_local_2%PHI, &
-  !b1Refield_2d_local_2%Z,b1Imfield_2d_local_2%R,b1Imfield_2d_local_2%PHI,b1Imfield_2d_local_2%Z, &
-  !Y_R,Y_Z,A,B1Re_R,B1Re_PHI,B1Re_Z,B1Im_R,B1Im_PHI,B1Im_Z,ezerr_local)
-  !call EZspline_error(ezerr_local)
+  E1_R = gr(1)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(E1Re_R*cPshift-E1Im_R*sPshift)
+  E1_PHI = gr(1)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
+  E1_Z = gr(1)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(1)))*(E1Re_Z*cPshift-E1Im_Z*sPshift)
+
+  E_R = E1_R
+  E_PHI = E1_PHI
+  E_Z = E1_Z
+
+  call EZspline_interp2_FOmarsEM(bfield_2d_local%A, &
+    b1Refield_2d_local_2%R,b1Refield_2d_local_2%PHI,b1Refield_2d_local_2%Z, &
+    b1Imfield_2d_local_2%R,b1Imfield_2d_local_2%PHI,b1Imfield_2d_local_2%Z, &
+    e1Refield_2d_local_2%R,e1Refield_2d_local_2%PHI,e1Refield_2d_local_2%Z, &
+    e1Imfield_2d_local_2%R,e1Imfield_2d_local_2%PHI,e1Imfield_2d_local_2%Z, &
+    Y_R,Y_Z,A, &
+    B1Re_R,B1Re_PHI,B1Re_Z,B1Im_R,B1Im_PHI,B1Im_Z, &
+    E1Re_R,E1Re_PHI,E1Re_Z,E1Im_R,E1Im_PHI,E1Im_Z, &
+    ezerr_local)
+  call EZspline_error(ezerr_local)
+
+  cPshift=cos(nmode(2)*Y_PHI-phase-fr(2)*time)
+  sPshift=sin(nmode(2)*Y_PHI-phase-fr(2)*time)
+
+  B1_R = gr(2)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(2)))*(B1Re_R*cPshift-B1Im_R*sPshift)
+  B1_PHI = gr(2)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(2)))*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
+  B1_Z = gr(2)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(2)))*(B1Re_Z*cPshift-B1Im_Z*sPshift)
+
+  B_R = B_R+B1_R
+  B_PHI = B_PHI+B1_PHI
+  B_Z = B_Z+B1_Z
+
+  E1_R = gr(2)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(2)))*(E1Re_R*cPshift-E1Im_R*sPshift)
+  E1_PHI = gr(2)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(2)))*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
+  E1_Z = gr(2)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(2)))*(E1Re_Z*cPshift-E1Im_Z*sPshift)
+
+  E_R = E_R+E1_R
+  E_PHI = E_PHI+E1_PHI
+  E_Z = E_Z+E1_Z
+
+  call EZspline_interp2_FOmarsEM(bfield_2d_local%A, &
+    b1Refield_2d_local_3%R,b1Refield_2d_local_3%PHI,b1Refield_2d_local_3%Z, &
+    b1Imfield_2d_local_3%R,b1Imfield_2d_local_3%PHI,b1Imfield_2d_local_3%Z, &
+    e1Refield_2d_local_3%R,e1Refield_2d_local_3%PHI,e1Refield_2d_local_3%Z, &
+    e1Imfield_2d_local_3%R,e1Imfield_2d_local_3%PHI,e1Imfield_2d_local_3%Z, &
+    Y_R,Y_Z,A, &
+    B1Re_R,B1Re_PHI,B1Re_Z,B1Im_R,B1Im_PHI,B1Im_Z, &
+    E1Re_R,E1Re_PHI,E1Re_Z,E1Im_R,E1Im_PHI,E1Im_Z, &
+    ezerr_local)
+  call EZspline_error(ezerr_local)
 
-  !cPshift=cos(nmode(2)*Y_PHI-phase)
-  !sPshift=sin(nmode(2)*Y_PHI-phase)
+  cPshift=cos(nmode(3)*Y_PHI-phase-fr(3)*time)
+  sPshift=sin(nmode(3)*Y_PHI-phase-fr(3)*time)
 
-  !B1_R = MARS_max/(1+MARS_max*exp(-time*gr(2)))*amp(2)*(B1Re_R*cPshift-B1Im_R*sPshift)
-  !B1_PHI = MARS_max/(1+MARS_max*exp(-time*gr(2)))*amp(2)*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
-  !B1_Z = MARS_max/(1+MARS_max*exp(-time*gr(2)))*amp(2)*(B1Re_Z*cPshift-B1Im_Z*sPshift)
+  B1_R = gr(3)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(3)))*(B1Re_R*cPshift-B1Im_R*sPshift)
+  B1_PHI = gr(3)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(3)))*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
+  B1_Z = gr(3)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(3)))*(B1Re_Z*cPshift-B1Im_Z*sPshift)
 
-  !B_R = B_R+B1_R
-  !B_PHI = B_PHI+B1_PHI
-  !B_Z = B_Z+B1_Z
+  B_R = B_R+B1_R
+  B_PHI = B_PHI+B1_PHI
+  B_Z = B_Z+B1_Z
 
-  !call EZspline_interp2_FOmars(bfield_2d_local%A,b1Refield_2d_local_3%R,b1Refield_2d_local_3%PHI, &
-  !b1Refield_2d_local_3%Z,b1Imfield_2d_local_3%R,b1Imfield_2d_local_3%PHI,b1Imfield_2d_local_3%Z, &
-  !Y_R,Y_Z,A,B1Re_R,B1Re_PHI,B1Re_Z,B1Im_R,B1Im_PHI,B1Im_Z,ezerr_local)
-  !call EZspline_error(ezerr_local)
+  E1_R = gr(3)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(3)))*(E1Re_R*cPshift-E1Im_R*sPshift)
+  E1_PHI = gr(3)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(3)))*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
+  E1_Z = gr(3)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(3)))*(E1Re_Z*cPshift-E1Im_Z*sPshift)
+
+  E_R = E_R+E1_R
+  E_PHI = E_PHI+E1_PHI
+  E_Z = E_Z+E1_Z
+
+  call EZspline_interp2_FOmarsEM(bfield_2d_local%A, &
+    b1Refield_2d_local_4%R,b1Refield_2d_local_4%PHI,b1Refield_2d_local_4%Z, &
+    b1Imfield_2d_local_4%R,b1Imfield_2d_local_4%PHI,b1Imfield_2d_local_4%Z, &
+    e1Refield_2d_local_4%R,e1Refield_2d_local_4%PHI,e1Refield_2d_local_4%Z, &
+    e1Imfield_2d_local_4%R,e1Imfield_2d_local_4%PHI,e1Imfield_2d_local_4%Z, &
+    Y_R,Y_Z,A, &
+    B1Re_R,B1Re_PHI,B1Re_Z,B1Im_R,B1Im_PHI,B1Im_Z, &
+    E1Re_R,E1Re_PHI,E1Re_Z,E1Im_R,E1Im_PHI,E1Im_Z, &
+    ezerr_local)
+  call EZspline_error(ezerr_local)
 
-  !cPshift=cos(nmode(3)*Y_PHI-phase)
-  !sPshift=sin(nmode(3)*Y_PHI-phase)
+  cPshift=cos(nmode(4)*Y_PHI-phase-fr(4)*time)
+  sPshift=sin(nmode(4)*Y_PHI-phase-fr(4)*time)
 
-  !B1_R = MARS_max/(1+MARS_max*exp(-time*gr(3)))*amp(3)*(B1Re_R*cPshift-B1Im_R*sPshift)
-  !B1_PHI = MARS_max/(1+MARS_max*exp(-time*gr(3)))*amp(3)*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
-  !B1_Z = MARS_max/(1+MARS_max*exp(-time*gr(3)))*amp(3)*(B1Re_Z*cPshift-B1Im_Z*sPshift)
+  B1_R = gr(4)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(4)))*(B1Re_R*cPshift-B1Im_R*sPshift)
+  B1_PHI = gr(4)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(4)))*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
+  B1_Z = gr(4)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(4)))*(B1Re_Z*cPshift-B1Im_Z*sPshift)
 
-  !B_R = B_R+B1_R
-  !B_PHI = B_PHI+B1_PHI
-  !B_Z = B_Z+B1_Z
+  B_R = B_R+B1_R
+  B_PHI = B_PHI+B1_PHI
+  B_Z = B_Z+B1_Z
+
+  E1_R = gr(4)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(4)))*(E1Re_R*cPshift-E1Im_R*sPshift)
+  E1_PHI = gr(4)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(4)))*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
+  E1_Z = gr(4)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(4)))*(E1Re_Z*cPshift-E1Im_Z*sPshift)
+
+  E_R = E_R+E1_R
+  E_PHI = E_PHI+E1_PHI
+  E_Z = E_Z+E1_Z
+
+  call EZspline_interp2_FOmarsEM(bfield_2d_local%A, &
+    b1Refield_2d_local_5%R,b1Refield_2d_local_5%PHI,b1Refield_2d_local_5%Z, &
+    b1Imfield_2d_local_5%R,b1Imfield_2d_local_5%PHI,b1Imfield_2d_local_5%Z, &
+    e1Refield_2d_local_5%R,e1Refield_2d_local_5%PHI,e1Refield_2d_local_5%Z, &
+    e1Imfield_2d_local_5%R,e1Imfield_2d_local_5%PHI,e1Imfield_2d_local_5%Z, &
+    Y_R,Y_Z,A, &
+    B1Re_R,B1Re_PHI,B1Re_Z,B1Im_R,B1Im_PHI,B1Im_Z, &
+    E1Re_R,E1Re_PHI,E1Re_Z,E1Im_R,E1Im_PHI,E1Im_Z, &
+    ezerr_local)
+  call EZspline_error(ezerr_local)
+
+  cPshift=cos(nmode(5)*Y_PHI-phase-fr(5)*time)
+  sPshift=sin(nmode(5)*Y_PHI-phase-fr(5)*time)
+
+  B1_R = gr(5)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(5)))*(B1Re_R*cPshift-B1Im_R*sPshift)
+  B1_PHI = gr(5)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(5)))*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
+  B1_Z = gr(5)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(5)))*(B1Re_Z*cPshift-B1Im_Z*sPshift)
+
+  B_R = B_R+B1_R
+  B_PHI = B_PHI+B1_PHI
+  B_Z = B_Z+B1_Z
+
+  E1_R = gr(5)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(5)))*(E1Re_R*cPshift-E1Im_R*sPshift)
+  E1_PHI = gr(5)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(5)))*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
+  E1_Z = gr(5)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(5)))*(E1Re_Z*cPshift-E1Im_Z*sPshift)
+
+  E_R = E_R+E1_R
+  E_PHI = E_PHI+E1_PHI
+  E_Z = E_Z+E1_Z
+
+  call EZspline_interp2_FOmarsEM(bfield_2d_local%A, &
+    b1Refield_2d_local_6%R,b1Refield_2d_local_6%PHI,b1Refield_2d_local_6%Z, &
+    b1Imfield_2d_local_6%R,b1Imfield_2d_local_6%PHI,b1Imfield_2d_local_6%Z, &
+    e1Refield_2d_local_6%R,e1Refield_2d_local_6%PHI,e1Refield_2d_local_6%Z, &
+    e1Imfield_2d_local_6%R,e1Imfield_2d_local_6%PHI,e1Imfield_2d_local_6%Z, &
+    Y_R,Y_Z,A, &
+    B1Re_R,B1Re_PHI,B1Re_Z,B1Im_R,B1Im_PHI,B1Im_Z, &
+    E1Re_R,E1Re_PHI,E1Re_Z,E1Im_R,E1Im_PHI,E1Im_Z, &
+    ezerr_local)
+  call EZspline_error(ezerr_local)
+
+  cPshift=cos(nmode(6)*Y_PHI-phase-fr(6)*time)
+  sPshift=sin(nmode(6)*Y_PHI-phase-fr(6)*time)
+
+  B1_R = gr(6)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(6)))*(B1Re_R*cPshift-B1Im_R*sPshift)
+  B1_PHI = gr(6)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(6)))*(B1Re_PHI*cPshift-B1Im_PHI*sPshift)
+  B1_Z = gr(6)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(6)))*(B1Re_Z*cPshift-B1Im_Z*sPshift)
+
+  B_R = B_R+B1_R
+  B_PHI = B_PHI+B1_PHI
+  B_Z = B_Z+B1_Z
+
+  E1_R = gr(6)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(6)))*(E1Re_R*cPshift-E1Im_R*sPshift)
+  E1_PHI = gr(6)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(6)))*(E1Re_PHI*cPshift-E1Im_PHI*sPshift)
+  E1_Z = gr(6)/MARS_quas_fac/(1+(MARS_max-1)*exp(-time*gr(6)))*(E1Re_Z*cPshift-E1Im_Z*sPshift)
+
+  E_R = E_R+E1_R
+  E_PHI = E_PHI+E1_PHI
+  E_Z = E_Z+E1_Z
 
   B_X = B_R*cP - B_PHI*sP
   B_Y = B_R*sP + B_PHI*cP
 
+  E_X = E_R*cP - E_PHI*sP
+  E_Y = E_R*sP + E_PHI*cP
+
 end subroutine interp_FOfields_marsEM_p_ACC
 
 subroutine provide_ezspline_aorsa_ACC(bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local, &
diff --git a/src/korc_ppusher.f90 b/src/korc_ppusher.f90
index 67778247..92ae3b29 100755
--- a/src/korc_ppusher.f90
+++ b/src/korc_ppusher.f90
@@ -31,6 +31,18 @@ module korc_ppusher
   TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: b1Imfield_2d_local_5
   TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: b1Refield_2d_local_6
   TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: b1Imfield_2d_local_6
+  TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: e1Refield_2d_local_1
+  TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: e1Imfield_2d_local_1
+  TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: e1Refield_2d_local_2
+  TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: e1Imfield_2d_local_2
+  TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: e1Refield_2d_local_3
+  TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: e1Imfield_2d_local_3
+  TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: e1Refield_2d_local_4
+  TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: e1Imfield_2d_local_4
+  TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: e1Refield_2d_local_5
+  TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: e1Imfield_2d_local_5
+  TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: e1Refield_2d_local_6
+  TYPE(KORC_2D_FIELDS_INTERPOLANT)      :: e1Imfield_2d_local_6
   TYPE(KORC_INTERPOLANT_DOMAIN)        :: fields_domain_local
   TYPE(KORC_2DX_FIELDS_INTERPOLANT)      :: b1Refield_2dx_local
   TYPE(KORC_2DX_FIELDS_INTERPOLANT)      :: b1Imfield_2dx_local
@@ -1340,12 +1352,12 @@ subroutine FO_init_marsEM_ACC(params,F,spp,output,step)
   REAL(rp) :: E_X,E_Y,E_Z
   REAL(rp) :: PSIp
   REAL(rp) :: m_cache,q_cache,psip_conv,phase,Ro,Bo,circumradius,ntiles,mmode,omega
-  REAL(rp),DIMENSION(3) :: amp,gr,nmode
+  REAL(rp),DIMENSION(6) :: amp,gr,fr,nmode
   INTEGER(is) :: flagCon,flagCol
   LOGICAL :: Analytic_D3D_IWL,useDiMES,Dim2x1t
   REAL(rp),DIMENSION(2) :: DiMESdims
   REAL(rp),DIMENSION(3) :: DiMESloc_cyl
-  REAL(rp) :: time,MARS_max
+  REAL(rp) :: time,MARS_max,MARS_quas_fac
 
   !$acc routine (cart_to_cyl_p_ACC) seq
   !$acc routine (interp_FOfields_marsEM_p_ACC) seq
@@ -1359,9 +1371,11 @@ subroutine FO_init_marsEM_ACC(params,F,spp,output,step)
     psip_conv=F%psip_conv
     amp=F%AMP
     gr=F%GR
+    fr=F%FR
     nmode=F%X%PHI
     phase=F%MARS_phase
     MARS_max=F%MARS_max
+    MARS_quas_fac=F%MARS_quas_fac
     Ro=F%Ro
     Bo=F%Bo
 
@@ -1378,14 +1392,22 @@ subroutine FO_init_marsEM_ACC(params,F,spp,output,step)
     if(output) then
 
 
-      call provide_ezspline_marsNL_ACC(bfield_2d_local,b1Refield_2d_local_1,b1Imfield_2d_local_1, &
+      call provide_ezspline_marsEM_ACC(bfield_2d_local,b1Refield_2d_local_1,b1Imfield_2d_local_1, &
         b1Refield_2d_local_2,b1Imfield_2d_local_2,b1Refield_2d_local_3,b1Imfield_2d_local_3, &
-        fields_domain_local)
+        b1Refield_2d_local_4,b1Imfield_2d_local_4,b1Refield_2d_local_5,b1Imfield_2d_local_5, &
+        b1Refield_2d_local_6,b1Imfield_2d_local_6,e1Refield_2d_local_1,e1Imfield_2d_local_1, &
+        e1Refield_2d_local_2,e1Imfield_2d_local_2,e1Refield_2d_local_3,e1Imfield_2d_local_3, &
+        e1Refield_2d_local_4,e1Imfield_2d_local_4,e1Refield_2d_local_5,e1Imfield_2d_local_5, &
+        e1Refield_2d_local_6,e1Imfield_2d_local_6,fields_domain_local)
 
 
       !$acc  enter data copyin(bfield_2d_local,b1Refield_2d_local_1,b1Imfield_2d_local_1, &
       !$acc& b1Refield_2d_local_2,b1Imfield_2d_local_2,b1Refield_2d_local_3,b1Imfield_2d_local_3, &
-      !$acc& fields_domain_local)
+      !$acc& b1Refield_2d_local_4,b1Imfield_2d_local_4,b1Refield_2d_local_5,b1Imfield_2d_local_5, &
+      !$acc& b1Refield_2d_local_6,b1Imfield_2d_local_6,e1Refield_2d_local_1,e1Imfield_2d_local_1, &
+      !$acc& e1Refield_2d_local_2,e1Imfield_2d_local_2,e1Refield_2d_local_3,e1Imfield_2d_local_3, &
+      !$acc& e1Refield_2d_local_4,e1Imfield_2d_local_4,e1Refield_2d_local_5,e1Imfield_2d_local_5, &
+      !$acc& e1Refield_2d_local_6,e1Imfield_2d_local_6,fields_domain_local)
 
       !$acc  parallel loop
       do pp=1_idef,spp(ii)%ppp
@@ -1423,10 +1445,21 @@ subroutine FO_init_marsEM_ACC(params,F,spp,output,step)
           Dim2x1t,Analytic_D3D_IWL,circumradius, &
           ntiles,useDiMES,DiMESloc_cyl,DiMESdims,Y_R,Y_PHI,Y_Z,flagCon)
 
-
-        call interp_FOfields_marsNL_p_ACC(time,bfield_2d_local,b1Refield_2d_local_1,b1Imfield_2d_local_1, &
-          b1Refield_2d_local_2,b1Imfield_2d_local_2,b1Refield_2d_local_3,b1Imfield_2d_local_3, &
-          psip_conv,amp,gr,nmode,phase,MARS_max,Bo,Ro,Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z,PSIp)
+        call interp_FOfields_marsEM_p_ACC(time,bfield_2d_local, &
+          b1Refield_2d_local_1,b1Imfield_2d_local_1, &
+          b1Refield_2d_local_2,b1Imfield_2d_local_2, &
+          b1Refield_2d_local_3,b1Imfield_2d_local_3, &
+          b1Refield_2d_local_4,b1Imfield_2d_local_4, &
+          b1Refield_2d_local_5,b1Imfield_2d_local_5, &
+          b1Refield_2d_local_6,b1Imfield_2d_local_6, &
+          e1Refield_2d_local_1,e1Imfield_2d_local_1, &
+          e1Refield_2d_local_2,e1Imfield_2d_local_2, &
+          e1Refield_2d_local_3,e1Imfield_2d_local_3, &
+          e1Refield_2d_local_4,e1Imfield_2d_local_4, &
+          e1Refield_2d_local_5,e1Imfield_2d_local_5, &
+          e1Refield_2d_local_6,e1Imfield_2d_local_6, &
+          psip_conv,amp,gr,fr,nmode,phase,MARS_max,MARS_quas_fac,Bo,Ro,Y_R,Y_PHI,Y_Z, &
+          B_X,B_Y,B_Z,E_X,E_Y,E_Z,PSIp)
 
 
 #endif PSPLINE
@@ -1511,7 +1544,11 @@ subroutine FO_init_marsEM_ACC(params,F,spp,output,step)
 
       !$acc  exit data delete(bfield_2d_local,b1Refield_2d_local_1,b1Imfield_2d_local_1, &
       !$acc& b1Refield_2d_local_2,b1Imfield_2d_local_2,b1Refield_2d_local_3,b1Imfield_2d_local_3, &
-      !$acc& fields_domain_local)
+      !$acc& b1Refield_2d_local_4,b1Imfield_2d_local_4,b1Refield_2d_local_5,b1Imfield_2d_local_5, &
+      !$acc& b1Refield_2d_local_6,b1Imfield_2d_local_6,e1Refield_2d_local_1,e1Imfield_2d_local_1, &
+      !$acc& e1Refield_2d_local_2,e1Imfield_2d_local_2,e1Refield_2d_local_3,e1Imfield_2d_local_3, &
+      !$acc& e1Refield_2d_local_4,e1Imfield_2d_local_4,e1Refield_2d_local_5,e1Imfield_2d_local_5, &
+      !$acc& e1Refield_2d_local_6,e1Imfield_2d_local_6,fields_domain_local)
 
     endif !(if output)
 
@@ -3830,8 +3867,8 @@ subroutine adv_FOinterp_marsEM_top_ACC(params,F,P,spp)
   INTEGER(is) :: flagCon,flagCol
   INTEGER(ip) :: tskip,it
   REAL(rp) :: a,m_cache,q_cache,psip_conv,phase
-  REAL(rp),DIMENSION(3) :: amp,gr,nmode
-  REAL(rp) :: Ro,Bo,circumradius,ntiles,dt,init_time,time,MARS_max,t_norm
+  REAL(rp),DIMENSION(6) :: amp,gr,fr,nmode
+  REAL(rp) :: Ro,Bo,circumradius,ntiles,dt,init_time,time,MARS_max,MARS_quas_fac,t_norm
   INTEGER  :: ii,pp,ss,tt,ppp
   LOGICAL :: Analytic_D3D_IWL,useDiMES,Dim2x1t
   REAL(rp),DIMENSION(2) :: DiMESdims
@@ -3859,11 +3896,13 @@ subroutine adv_FOinterp_marsEM_top_ACC(params,F,P,spp)
     psip_conv=F%psip_conv
     amp=F%AMP
     gr=F%GR
+    fr=F%FR
     nmode=F%X%PHI
     phase=F%MARS_phase
     Ro=F%Ro
     Bo=F%Bo
     MARS_max=F%MARS_max
+    MARS_quas_fac=F%MARS_quas_fac
 
     Dim2x1t=F%Dim2x1t
     Analytic_D3D_IWL=F%Analytic_D3D_IWL
@@ -3880,6 +3919,12 @@ subroutine adv_FOinterp_marsEM_top_ACC(params,F,P,spp)
     !$acc& b1Refield_2d_local_4,b1Imfield_2d_local_4, &
     !$acc& b1Refield_2d_local_5,b1Imfield_2d_local_5, &
     !$acc& b1Refield_2d_local_6,b1Imfield_2d_local_6, &
+    !$acc& e1Refield_2d_local_1,e1Imfield_2d_local_1, &
+    !$acc& e1Refield_2d_local_2,e1Imfield_2d_local_2, &
+    !$acc& e1Refield_2d_local_3,e1Imfield_2d_local_3, &
+    !$acc& e1Refield_2d_local_4,e1Imfield_2d_local_4, &
+    !$acc& e1Refield_2d_local_5,e1Imfield_2d_local_5, &
+    !$acc& e1Refield_2d_local_6,e1Imfield_2d_local_6, &
     !$acc& fields_domain_local)
 
     !$acc parallel loop 
@@ -3923,8 +3968,21 @@ subroutine adv_FOinterp_marsEM_top_ACC(params,F,P,spp)
 
         time=(init_time+(it-1+tt)*dt)*t_norm
 
-        call interp_FOfields_marsEM_p_ACC(time,bfield_2d_local,b1Refield_2d_local_1,b1Imfield_2d_local_1, &
-          psip_conv,amp,gr,nmode,phase,MARS_max,Bo,Ro,Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z,PSIp)
+        call interp_FOfields_marsEM_p_ACC(time,bfield_2d_local, &
+          b1Refield_2d_local_1,b1Imfield_2d_local_1, &
+          b1Refield_2d_local_2,b1Imfield_2d_local_2, &
+          b1Refield_2d_local_3,b1Imfield_2d_local_3, &
+          b1Refield_2d_local_4,b1Imfield_2d_local_4, &
+          b1Refield_2d_local_5,b1Imfield_2d_local_5, &
+          b1Refield_2d_local_6,b1Imfield_2d_local_6, &
+          e1Refield_2d_local_1,e1Imfield_2d_local_1, &
+          e1Refield_2d_local_2,e1Imfield_2d_local_2, &
+          e1Refield_2d_local_3,e1Imfield_2d_local_3, &
+          e1Refield_2d_local_4,e1Imfield_2d_local_4, &
+          e1Refield_2d_local_5,e1Imfield_2d_local_5, &
+          e1Refield_2d_local_6,e1Imfield_2d_local_6, &
+          psip_conv,amp,gr,fr,nmode,phase,MARS_max,MARS_quas_fac,Bo,Ro,Y_R,Y_PHI,Y_Z, &
+          B_X,B_Y,B_Z,E_X,E_Y,E_Z,PSIp)
 
         call advance_FO_vars_ACC(dt,tt,a,q_cache,m_cache, &
             X_X,X_Y,X_Z,V_X,V_Y,V_Z,B_X,B_Y,B_Z,E_X,E_Y,E_Z, &
@@ -4024,6 +4082,12 @@ subroutine adv_FOinterp_marsEM_top_ACC(params,F,P,spp)
     !$acc& b1Refield_2d_local_4,b1Imfield_2d_local_4, &
     !$acc& b1Refield_2d_local_5,b1Imfield_2d_local_5, &
     !$acc& b1Refield_2d_local_6,b1Imfield_2d_local_6, &
+    !$acc& e1Refield_2d_local_1,e1Imfield_2d_local_1, &
+    !$acc& e1Refield_2d_local_2,e1Imfield_2d_local_2, &
+    !$acc& e1Refield_2d_local_3,e1Imfield_2d_local_3, &
+    !$acc& e1Refield_2d_local_4,e1Imfield_2d_local_4, &
+    !$acc& e1Refield_2d_local_5,e1Imfield_2d_local_5, &
+    !$acc& e1Refield_2d_local_6,e1Imfield_2d_local_6, &
     !$acc& fields_domain_local)
 
   end do !species iterator
diff --git a/src/korc_units.f90 b/src/korc_units.f90
index 8263fc49..a8444326 100755
--- a/src/korc_units.f90
+++ b/src/korc_units.f90
@@ -394,7 +394,7 @@ subroutine normalize_variables(params,spp,F,P)
 
      if (F%B1field) then
 
-        if (params%field_model(10:13).eq.'MARS') then
+        if (params%field_model(10:14).eq.'MARS ') then
 
            if (ALLOCATED(F%B1Re_2D%R)) F%B1Re_2D%R = F%B1Re_2D%R/ &
                 params%cpp%Bo
@@ -459,6 +459,22 @@ subroutine normalize_variables(params,spp,F,P)
         if (ALLOCATED(F%E1Im_2DX%Z)) F%E1Im_2DX%Z = F%E1Im_2DX%Z/ &
              params%cpp%Eo
 
+        if (params%field_model(10:16).eq.'MARS_EM') then
+
+           if (ALLOCATED(F%E1Re_3D%R)) F%E1Re_3D%R = F%E1Re_3D%R/ &
+                params%cpp%Eo
+           if (ALLOCATED(F%E1Re_3D%PHI)) F%E1Re_3D%PHI = F%E1Re_3D%PHI/ &
+                params%cpp%Eo
+           if (ALLOCATED(F%E1Re_3D%Z)) F%E1Re_3D%Z = F%E1Re_3D%Z/ &
+                params%cpp%Eo
+           if (ALLOCATED(F%E1Im_3D%R)) F%E1Im_3D%R = F%E1Im_3D%R/ &
+                params%cpp%Eo
+           if (ALLOCATED(F%E1Im_3D%PHI)) F%E1Im_3D%PHI = F%E1Im_3D%PHI/ &
+                params%cpp%Eo
+           if (ALLOCATED(F%E1Im_3D%Z)) F%E1Im_3D%Z = F%E1Im_3D%Z/ &
+                params%cpp%Eo     
+        endif 
+
      end if
 
      F%X%R = F%X%R/params%cpp%length
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 7c5ec8b2..ae041fa9 100755
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -23,14 +23,14 @@ set_property(TARGET xtest PROPERTY LINKER_LANGUAGE Fortran)
 configure_file(${CMAKE_SOURCE_DIR}/test/egyro/korc_egyro.sh.in ${CMAKE_BINARY_DIR}/egyro_test/korc_egyro.sh)
 configure_file(${CMAKE_SOURCE_DIR}/test/elong_trans/korc_elong_trans.sh.in ${CMAKE_BINARY_DIR}/elong_trans_test/korc_elong_trans.sh)
 configure_file(${CMAKE_SOURCE_DIR}/test/mars/korc_mars.sh.in ${CMAKE_BINARY_DIR}/mars_test/korc_mars.sh)
-#configure_file(${CMAKE_SOURCE_DIR}/test/mars_EM/korc_mars_EM.sh.in ${CMAKE_BINARY_DIR}/mars_EM_test/korc_mars_EM.sh)
+configure_file(${CMAKE_SOURCE_DIR}/test/mars_EM/korc_mars_EM.sh.in ${CMAKE_BINARY_DIR}/mars_EM_test/korc_mars_EM.sh)
 configure_file(${CMAKE_SOURCE_DIR}/test/aorsa/korc_aorsa.sh.in ${CMAKE_BINARY_DIR}/aorsa_test/korc_aorsa.sh)
 configure_file(${CMAKE_SOURCE_DIR}/test/fio_m3dc1/korc_fio_m3dc1.sh.in ${CMAKE_BINARY_DIR}/fio_m3dc1_test/korc_fio_m3dc1.sh)
 
 file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bin/egyro_test)
 file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bin/elong_trans_test)
 file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bin/mars_test)
-#file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bin/mars_EM_test)
+file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bin/mars_EM_test)
 file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bin/aorsa_test)
 file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bin/fio_m3dc1_test)
 
@@ -52,10 +52,10 @@ foreach(RANK IN ITEMS 1 2 4 8)
                 WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
       set_tests_properties (mars_test_${RANK} PROPERTIES PROCESSORS ${RANK} ENVIRONMENT OMP_NUM_THREADS=1)
 
-      #add_test (NAME mars_EM_test_${RANK}
-      #    COMMAND ${CMAKE_BINARY_DIR}/mars_EM_test/korc_mars_EM.sh ${RANK}
-      #          WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
-      #set_tests_properties (mars_test_EM_${RANK} PROPERTIES PROCESSORS ${RANK} ENVIRONMENT OMP_NUM_THREADS=1)
+      add_test (NAME mars_EM_test_${RANK}
+          COMMAND ${CMAKE_BINARY_DIR}/mars_EM_test/korc_mars_EM.sh ${RANK}
+                WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
+      set_tests_properties (mars_EM_test_${RANK} PROPERTIES PROCESSORS ${RANK} ENVIRONMENT OMP_NUM_THREADS=1)
 
       add_test (NAME aorsa_test_${RANK}
       COMMAND ${CMAKE_BINARY_DIR}/aorsa_test/korc_aorsa.sh ${RANK}
diff --git a/test/mars_EM/input_file_MARS_EM.korc b/test/mars_EM/input_file_MARS_EM.korc
index ea45f1a6..f2b17398 100644
--- a/test/mars_EM/input_file_MARS_EM.korc
+++ b/test/mars_EM/input_file_MARS_EM.korc
@@ -5,17 +5,17 @@
     ! Append simulation results after previous simulation_time reached
   reinit = F
     ! Begin a new simulation, reinitializing from restart file state
-  simulation_time = 0.E-6
+  simulation_time = 1.E-9
     ! Total aimed simulation time in seconds   
     ! Run 10 mu s If transients exist put 5 extra mu s.
-  snapshot_frequency = 1.E-7
+  snapshot_frequency = 1.E-9
     ! Time between snapshots in seconds
   dt = 1.E-2
     ! Time step as fraction of relativistic gyro-period
   minimum_particle_energy = 1.0E6 
     ! Minimum energy of simulated particles in eV
   radiation = F
-  GC_rad_model='SDE'
+  GC_rad_SDE= F
   collisions = F
   bound_electron_model = 'HESSLOW' 
     ! Options are: 'NO_BOUND', 'HESSLOW', and 'ROSENBLUTH'
@@ -41,17 +41,19 @@
     ! Set for plasma_model='ANALYTICAL'. Can be 'interp' or 'eqn',
     ! where 'eqn' evaluates particle fields at particle positions and
     ! 'interp' interpolates precomputed fields.
+  SameRandSeed = T
 /
 
 &plasma_species
   runaway = T
-  ppp = 10000
+  ppp = 100
+  pinit = 100
     ! Number of particles per process (mpi)
   q = -1.0
     ! Electric charge
   m = 1.0
     ! In units of electron mass
-  spatial_distribution = 'MH_psi'
+  spatial_distribution = 'TRACER'
     ! Options are: 'UNIFORM', 'DISK', 'TORUS', 'EXPONENTIAL-TORUS',
     ! 'GAUSSIAN-TORUS', 'ELLIPTIC-TORUS', 'EXPONENTIAL-ELLIPTIC-TORUS',
     ! 'GAUSSIAN-ELLIPTICAL-TORUS', '2D-GAUSSIAN-ELLIPTIC-TORUS-MH',
@@ -66,7 +68,7 @@
   sigmaR = 1.e6
   sigmaZ = 0.2
   theta_gauss = 0.0
-  psi_max=1.
+  psi_max=0.8
     ! goes as R^2 for HOLLMANN-3D, is psiN_max for HOLLMANN-3D-PSI
   energy_distribution = 'MONOENERGETIC'
     ! Options are: 'MONOENERGETIC', 'THERMAL', 'AVALANCHE', 
@@ -77,11 +79,7 @@
     ! Initial energy in eV
   etao = 170.
     ! Initial pitch angle
-  Eo_lims = 1.0E6,50.0E6
-    ! Lower and upper limit of simulated energy range, in eV.
-  etao_lims = 0.0,20.0
-    ! Lower and upper limit of simulated pitch-angle range, in degrees.
-  Xtrace = 1.6,0.0,-0.75
+  Xtrace = 5.,0.0,-0.0
     ! Initial position of tracer particle for debugging with
     ! spatial_distribution='TRACER'
   dth = 3.
@@ -116,27 +114,5 @@
   MARS_AMP_Scale=1.
   MARS_phase=0.
   MARS_max=100. !ratio of init to max
-/
-
-&plasmaProfiles
-  axisymmetric = T
-  filename = 'DIII_200236_MARS_CaseA.h5' 
-/
-
-&HollmannPDF
-  filename_Hollmann = 'Hollmann_PDF_DIII_200236.h5'
-  rho_ind = 4
-  max_pitch_angle_Hollmann = 90. 
-    ! In degrees
-  min_pitch_angle_Hollmann = 0. 
-    ! In degrees
-  min_energy_Hollmann = 2.E3 
-    ! In eV
-  max_energy_Hollmann = 3.E7
-    ! For Hollmann_PDF_HR.h5, needs to be less than 80 MeV
-    ! In eV
-  current_direction_Hollmann = 'ANTIPARALLEL'
-  Zeff_Hollmann = 1.
-  E_Hollmann = 1.
-
+  MARS_quas_fac = 1.E8
 /
\ No newline at end of file
diff --git a/test/mars_EM/rank_1/file_0.h5 b/test/mars_EM/rank_1/file_0.h5
new file mode 100644
index 00000000..c6efad8b
Binary files /dev/null and b/test/mars_EM/rank_1/file_0.h5 differ
diff --git a/test/mars_EM/rank_1/output.korc b/test/mars_EM/rank_1/output.korc
new file mode 100644
index 00000000..3537e719
Binary files /dev/null and b/test/mars_EM/rank_1/output.korc differ
diff --git a/test/mars_EM/rank_1/restart_file.h5 b/test/mars_EM/rank_1/restart_file.h5
new file mode 100644
index 00000000..aacb07a0
Binary files /dev/null and b/test/mars_EM/rank_1/restart_file.h5 differ
diff --git a/test/mars_EM/rank_1/simulation_parameters.h5 b/test/mars_EM/rank_1/simulation_parameters.h5
new file mode 100644
index 00000000..10108da5
Binary files /dev/null and b/test/mars_EM/rank_1/simulation_parameters.h5 differ
diff --git a/test/mars_EM/rank_2/file_0.h5 b/test/mars_EM/rank_2/file_0.h5
new file mode 100644
index 00000000..43e0e2aa
Binary files /dev/null and b/test/mars_EM/rank_2/file_0.h5 differ
diff --git a/test/mars_EM/rank_2/file_1.h5 b/test/mars_EM/rank_2/file_1.h5
new file mode 100644
index 00000000..43e0e2aa
Binary files /dev/null and b/test/mars_EM/rank_2/file_1.h5 differ
diff --git a/test/mars_EM/rank_2/output.korc b/test/mars_EM/rank_2/output.korc
new file mode 100644
index 00000000..26956c3c
Binary files /dev/null and b/test/mars_EM/rank_2/output.korc differ
diff --git a/test/mars_EM/rank_2/restart_file.h5 b/test/mars_EM/rank_2/restart_file.h5
new file mode 100644
index 00000000..8f189255
Binary files /dev/null and b/test/mars_EM/rank_2/restart_file.h5 differ
diff --git a/test/mars_EM/rank_2/simulation_parameters.h5 b/test/mars_EM/rank_2/simulation_parameters.h5
new file mode 100644
index 00000000..47a9d583
Binary files /dev/null and b/test/mars_EM/rank_2/simulation_parameters.h5 differ
diff --git a/test/mars_EM/rank_4/file_0.h5 b/test/mars_EM/rank_4/file_0.h5
new file mode 100644
index 00000000..1bd75268
Binary files /dev/null and b/test/mars_EM/rank_4/file_0.h5 differ
diff --git a/test/mars_EM/rank_4/file_1.h5 b/test/mars_EM/rank_4/file_1.h5
new file mode 100644
index 00000000..1bd75268
Binary files /dev/null and b/test/mars_EM/rank_4/file_1.h5 differ
diff --git a/test/mars_EM/rank_4/file_2.h5 b/test/mars_EM/rank_4/file_2.h5
new file mode 100644
index 00000000..1bd75268
Binary files /dev/null and b/test/mars_EM/rank_4/file_2.h5 differ
diff --git a/test/mars_EM/rank_4/file_3.h5 b/test/mars_EM/rank_4/file_3.h5
new file mode 100644
index 00000000..1bd75268
Binary files /dev/null and b/test/mars_EM/rank_4/file_3.h5 differ
diff --git a/test/mars_EM/rank_4/output.korc b/test/mars_EM/rank_4/output.korc
new file mode 100644
index 00000000..1634272a
Binary files /dev/null and b/test/mars_EM/rank_4/output.korc differ
diff --git a/test/mars_EM/rank_4/restart_file.h5 b/test/mars_EM/rank_4/restart_file.h5
new file mode 100644
index 00000000..82ff5be9
Binary files /dev/null and b/test/mars_EM/rank_4/restart_file.h5 differ
diff --git a/test/mars_EM/rank_4/simulation_parameters.h5 b/test/mars_EM/rank_4/simulation_parameters.h5
new file mode 100644
index 00000000..504918a6
Binary files /dev/null and b/test/mars_EM/rank_4/simulation_parameters.h5 differ
diff --git a/test/mars_EM/rank_8/file_0.h5 b/test/mars_EM/rank_8/file_0.h5
new file mode 100644
index 00000000..18303391
Binary files /dev/null and b/test/mars_EM/rank_8/file_0.h5 differ
diff --git a/test/mars_EM/rank_8/file_1.h5 b/test/mars_EM/rank_8/file_1.h5
new file mode 100644
index 00000000..18303391
Binary files /dev/null and b/test/mars_EM/rank_8/file_1.h5 differ
diff --git a/test/mars_EM/rank_8/file_2.h5 b/test/mars_EM/rank_8/file_2.h5
new file mode 100644
index 00000000..18303391
Binary files /dev/null and b/test/mars_EM/rank_8/file_2.h5 differ
diff --git a/test/mars_EM/rank_8/file_3.h5 b/test/mars_EM/rank_8/file_3.h5
new file mode 100644
index 00000000..18303391
Binary files /dev/null and b/test/mars_EM/rank_8/file_3.h5 differ
diff --git a/test/mars_EM/rank_8/file_4.h5 b/test/mars_EM/rank_8/file_4.h5
new file mode 100644
index 00000000..18303391
Binary files /dev/null and b/test/mars_EM/rank_8/file_4.h5 differ
diff --git a/test/mars_EM/rank_8/file_5.h5 b/test/mars_EM/rank_8/file_5.h5
new file mode 100644
index 00000000..18303391
Binary files /dev/null and b/test/mars_EM/rank_8/file_5.h5 differ
diff --git a/test/mars_EM/rank_8/file_6.h5 b/test/mars_EM/rank_8/file_6.h5
new file mode 100644
index 00000000..06e8e918
Binary files /dev/null and b/test/mars_EM/rank_8/file_6.h5 differ
diff --git a/test/mars_EM/rank_8/file_7.h5 b/test/mars_EM/rank_8/file_7.h5
new file mode 100644
index 00000000..06e8e918
Binary files /dev/null and b/test/mars_EM/rank_8/file_7.h5 differ
diff --git a/test/mars_EM/rank_8/output.korc b/test/mars_EM/rank_8/output.korc
new file mode 100644
index 00000000..5c0d12dd
Binary files /dev/null and b/test/mars_EM/rank_8/output.korc differ
diff --git a/test/mars_EM/rank_8/restart_file.h5 b/test/mars_EM/rank_8/restart_file.h5
new file mode 100644
index 00000000..2ed0e61c
Binary files /dev/null and b/test/mars_EM/rank_8/restart_file.h5 differ
diff --git a/test/mars_EM/rank_8/simulation_parameters.h5 b/test/mars_EM/rank_8/simulation_parameters.h5
new file mode 100644
index 00000000..9bfd801a
Binary files /dev/null and b/test/mars_EM/rank_8/simulation_parameters.h5 differ