diff --git a/maxima/g0/gk_collisionless/dg_gk-surf.mac b/maxima/g0/gk_collisionless/dg_gk-surf.mac
index 77901135..d8d97ea6 100644
--- a/maxima/g0/gk_collisionless/dg_gk-surf.mac
+++ b/maxima/g0/gk_collisionless/dg_gk-surf.mac
@@ -26,7 +26,7 @@ calcGKSurfUpdateInDir(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, vars
    BstarXdBmag_e,BstarYdBmag_e,BstarZdBmag_e,BstardBmag_e,
    hamil_e,alphaSurfL_e,alphaSurfR_e,
    fl_e,fc_e,fr_e,fUpL_e,fUpR_e,GhatL_c,GhatR_c,GhatL_e,GhatR_e,incrL_c,incrR_c,pOrderCFL,
-   fnodal_l_e, fnodal_r_e, fmodproj_e],
+   fnodal_l_e,fnodal_r_e,fmodproj_e,numC,vmap_prime_fac_l,vmap_prime_fac_c,vmap_prime_fac_r,vmap_prime_e,basisNodal,i],
 
   kill(varsC,varsP,bC,bP),
   pDim : cdim+vdim,
@@ -143,7 +143,8 @@ calcGKBoundarySurfUpdateInDir(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrd
    rdx2vec,rdv2vec,rdSurfVar2,bmagBasis,ignoreVars,inFlds_e,cmag_e,b_x_e,b_y_e,b_z_e,jacobTotInv_e,
    BstarXdBmag_e,BstarYdBmag_e,BstarZdBmag_e,BstardBmag_e,
    hamil_e,alphaUpL_e,alphaSurfL_e,alphaUpSurfL_e,alphaUpR_e,alphaSurfR_e,alphaUpSurfR_e,
-   fEdge_e,fSkin_e,fUpL_e,fUpR_e,GhatL_c,GhatR_c,GhatL_e,GhatR_e,incrL_c,incrR_c,pOrderCFL],
+   fEdge_e,fSkin_e,fUpL_e,fUpR_e,GhatL_c,GhatR_c,GhatL_e,GhatR_e,incrL_c,incrR_c,pOrderCFL,
+   numC,vmap_prime_fac_edge,vmap_prime_fac_skin,vmap_prime_e,basisNodal,i,fnodal_l_e,fnodal_r_e,fmodproj_e],
 
   kill(varsC,varsP,bC,bP),
   pDim : cdim+vdim,
diff --git a/maxima/g0/gk_collisionless/dg_gk-vol.mac b/maxima/g0/gk_collisionless/dg_gk-vol.mac
index 01905ad3..71e05911 100644
--- a/maxima/g0/gk_collisionless/dg_gk-vol.mac
+++ b/maxima/g0/gk_collisionless/dg_gk-vol.mac
@@ -11,18 +11,23 @@ load("utilities")$
 fpprec : 24$
 
 buildGKVolKernel(fh, funcNm, cdim, vdim, basisFun, polyOrder, varsInB, no_by) := block(
-  [pDim,varsC,bC,varsP,bP,varsV,vSub,numC,numP,varLabel,d,rdx2vec,rdv2vec,allVarLabelsC,
-   bmagBasis,ignoreVars,inFlds_e,cmag_e,b_x_e,b_y_e,b_z_e,jacobTotInv_e,vmap_e,BstardBmag_e,
-   hamil_e,pbAuxFlds,alphaSum_e,vd,dir,dirLabel,wDir,rdDirVar2,vmap_prime_fac,dirVar,
-   dirVar_phys,alpha_e,alpha_c,alphaLabel,alphaNoZero_c,alphaDotGradBasis_e,f_e,volTerm_c, dH_dz_e, alphaJf_e, Jf_e, replaceListHamil, replaceListVpar,hamil2_c,isqlist,mvpar_e,mvparsq_e],
+  [pdim,varsC,bC,varsP,bP,vSub,numC,numP,varLabel,d,rdx2vec,rdv2vec,allVarLabelsC,
+   bmagBasis,vmap_e,hamil_e,dir,dirLabel,alpha_e,alpha_c,alphaLabel,alphaNoZero_c,volTerm_c, 
+   alphaJf_e, Jf_e, replaceListHamil, replaceListVpar,isqlist,mvpar_e,mvparsq_e,
+   replaceList,dvparSimp,phi_e,bmag_e,rtg33inv_e,
+   dualcurlbhatoverB_1_e,dualcurlbhatoverB_2_e,dualcurlbhatoverB_3_e,dualcurlbhatoverB_vec,
+   bioverJB_1_e,bioverJB_2_e,bioverJB_3_e,bioverJB_vec,
+   vmapSq_e,vmap_prime_e,hamilCvar,hamilNoZero_c,hamil_c,
+   xidx,yidx,zidx,vidx,vpardim,dH_dx,dH_dy,dH_dz,dH_dvpar,gradH_vec, gradpsi, dpsidvpar,
+   boverBxgradH_x,boverBxgradH_y,boverBxgradH_z,boverBxgradH_vec,i,k,clst],
 
   kill(varsC,varsP,bC,bP),
-  pDim : cdim+vdim,
+  pdim : cdim+vdim,
 
   [varsC,bC,varsP,bP,vSub] : loadGkBasis(basisFun, cdim, vdim, polyOrder),
   numC : length(bC),  numP : length(bP), 
 
-  varLabel : makelist(string(varsP[d]),d,1,pDim),
+  varLabel : makelist(string(varsP[d]),d,1,pdim),
 
   print("Working on ", funcNm),
   printf(fh, "GKYL_CU_DH double ~a(const double *w, const double *dxv, const double *vmap, const double *vmapSq,
@@ -41,32 +46,41 @@ buildGKVolKernel(fh, funcNm, cdim, vdim, basisFun, polyOrder, varsInB, no_by) :=
   printf(fh, "~%"),
 
   /* Declare cell-center variables and variables multiplying gradients. */
-  for d : 1 thru pDim do (
+  for d : 1 thru pdim do (
     printf(fh, "  double rd~a2 = 2.0/dxv[~a];~%", varLabel[d], d-1)
   ),
   printf(fh, "~%"),
   rdx2vec : makelist(eval_string(sconcat("rd",varLabel[i],"2")),i,1,cdim),
-  rdv2vec : makelist(eval_string(sconcat("rd",varLabel[i],"2")),i,cdim+1,pDim),
+  rdv2vec : makelist(eval_string(sconcat("rd",varLabel[i],"2")),i,cdim+1,pdim),
 
   /* Declare variables with squared of cell centers and rdx2 variables (only need vpar^2). */
   printf(fh, "  double rdvpar2Sq = rdvpar2*rdvpar2;~%"),
   printf(fh, "  double dvparSq = dxv[~a]*dxv[~a];~%", cdim, cdim),
   printf(fh, "~%"),
   replaceList : [rdvpar2^2=rdvpar2Sq,dxv[cdim]^2=dvparSq,rdvpar2Sq=4/dvparSq],
-  dvparSimp : append(makelist(dxv[i-1]=2/eval_string(sconcat("rd",varLabel[i],"2")),i,1,pDim),
+  dvparSimp : append(makelist(dxv[i-1]=2/eval_string(sconcat("rd",varLabel[i],"2")),i,1,pdim),
                      [dvparSq=4/rdvpar2Sq]),
 
+  /* Store indices of the coordinate directions. */
+  if cdim = 3 then (
+    xidx : 1, yidx : 2, zidx : 3, vidx : 4
+  ) else if cdim = 2 then (
+    xidx : 1, zidx : 2, vidx : 3
+  ) else if cdim = 1 then (
+    zidx : 1, vidx : 2
+  ),
+
   /* Create pointers to the components of b_i. */
   allVarLabelsC : ["x","y","z"],
   for d : 1 thru 3 do (
-    printf(fh, "  const double *bioverJB_~a = &bioverJB[~a]; ~%", allVarLabelsC[d], numC*(d-1))
+    printf(fh, "  const double *bioverJB_~a = &bioverJB[~a]; ~%", d, numC*(d-1))
   ),
   printf(fh, "~%"),
 
   /* Create pointers to the components of dualcurlbhatoverB. */
   allVarLabelsC : ["x","y","z"],
   for d : 1 thru 3 do (
-    printf(fh, "  const double *dualcurlbhatoverB_~a = &dualcurlbhatoverB[~a]; ~%", allVarLabelsC[d], numC*(d-1))
+    printf(fh, "  const double *dualcurlbhatoverB_~a = &dualcurlbhatoverB[~a]; ~%", d, numC*(d-1))
   ),
   printf(fh, "~%"),
 
@@ -75,16 +89,24 @@ buildGKVolKernel(fh, funcNm, cdim, vdim, basisFun, polyOrder, varsInB, no_by) :=
   /* Expand input fields for Hamiltonian calculation */
   phi_e     : doExpand1(phi,bC),
   bmag_e : doExpand1(bmag, bmagBasis),
-  dualcurlbhatoverB_x_e : doExpand1(dualcurlbhatoverB_x, bmagBasis),
-  dualcurlbhatoverB_y_e : doExpand1(dualcurlbhatoverB_y, bmagBasis),
-  dualcurlbhatoverB_z_e : doExpand1(dualcurlbhatoverB_z, bmagBasis),
+  dualcurlbhatoverB_1_e : doExpand1(dualcurlbhatoverB_1, bmagBasis),
+  dualcurlbhatoverB_2_e : doExpand1(dualcurlbhatoverB_2, bmagBasis),
+  dualcurlbhatoverB_3_e : doExpand1(dualcurlbhatoverB_3, bmagBasis),
   rtg33inv_e : doExpand1(rtg33inv, bmagBasis),
-  bioverJB_x_e : doExpand1(bioverJB_x, bmagBasis),
-  bioverJB_y_e : doExpand1(bioverJB_y, bmagBasis),
-  bioverJB_z_e : doExpand1(bioverJB_z, bmagBasis),
-
-  dualcurlbhatoverB_list : [dualcurlbhatoverB_x_e, dualcurlbhatoverB_y_e, dualcurlbhatoverB_z_e],
-  bioverJB_list : [bioverJB_x_e, bioverJB_y_e, bioverJB_z_e],
+  bioverJB_1_e : doExpand1(bioverJB_1, bmagBasis),
+  bioverJB_2_e : doExpand1(bioverJB_2, bmagBasis),
+  bioverJB_3_e : doExpand1(bioverJB_3, bmagBasis),
+
+  if cdim = 3 then (
+    dualcurlbhatoverB_vec : [dualcurlbhatoverB_1_e, dualcurlbhatoverB_2_e, dualcurlbhatoverB_3_e],
+    bioverJB_vec : [bioverJB_1_e, bioverJB_2_e, bioverJB_3_e]
+  ) else if cdim = 2 then (
+    dualcurlbhatoverB_vec : [dualcurlbhatoverB_1_e, dualcurlbhatoverB_3_e],
+    bioverJB_vec : [bioverJB_1_e, bioverJB_3_e]
+  ) else if cdim = 1 then (
+    dualcurlbhatoverB_vec : [dualcurlbhatoverB_3_e],
+    bioverJB_vec : [bioverJB_3_e]
+  ),
 
   /* Velocity mapping fields. */
   [vmap_e,vmapSq_e,vmap_prime_e] : expandVmapFields(varsP),
@@ -110,16 +132,38 @@ buildGKVolKernel(fh, funcNm, cdim, vdim, basisFun, polyOrder, varsInB, no_by) :=
   Jf_e : doExpand1(fin,bP),
 
   /* Calculate expressions for dericatives of the hamiltonian*/
-  vpardim : pDim-1,
-  if vdim = 1 then ( vpardim : pDim ),
-  dH_dz_e : makelist(0, i, 1, pDim),  
-  for i : 1 thru vpardim do (
-    if i = vpardim then (
-      dH_dz_e[i] : diff(hamil_e,varsP[i])
-    )
-    else (
-      dH_dz_e[i] : diff(hamil_e*rdx2vec[i],varsP[i])
-    )
+  vpardim : pdim-1,
+  if vdim = 1 then ( vpardim : pdim ),
+  dH_dvpar : diff(hamil_e,varsP[vpardim])/vmap_prime_e[1], /* We do not multply by rdv2vec because it is inside vmap_prime */
+  if cdim = 3 then (
+    dH_dx : diff(hamil_e * rdx2vec[xidx],varsP[xidx]),
+    dH_dy : diff(hamil_e * rdx2vec[yidx],varsP[yidx]),
+    dH_dz : diff(hamil_e * rdx2vec[zidx],varsP[zidx]),
+    gradH_vec : [dH_dx, dH_dy, dH_dz]
+  ) else if cdim = 2 then (
+    dH_dx : diff(hamil_e * rdx2vec[xidx],varsP[xidx]),
+    dH_dy : 0,
+    dH_dz : diff(hamil_e * rdx2vec[zidx],varsP[zidx]),
+    gradH_vec : [dH_dx, dH_dz]
+  ) else if cdim = 1 then (
+    dH_dx : 0,
+    dH_dy : 0,
+    dH_dz : diff(hamil_e * rdx2vec[zidx],varsP[zidx]),
+    gradH_vec : [dH_dz]
+  ),
+
+  if cdim = 3 then (
+    boverBxgradH_x : bioverJB_2_e*dH_dz - bioverJB_3_e*dH_dy,
+    boverBxgradH_y : bioverJB_3_e*dH_dx - bioverJB_1_e*dH_dz,
+    boverBxgradH_z : bioverJB_1_e*dH_dy - bioverJB_2_e*dH_dx,
+    boverBxgradH_vec : [boverBxgradH_x, boverBxgradH_y, boverBxgradH_z]
+  ) else if cdim = 2 then (
+    boverBxgradH_x : bioverJB_2_e*dH_dz - bioverJB_3_e*dH_dy,
+    boverBxgradH_z : bioverJB_1_e*dH_dy - bioverJB_2_e*dH_dx,
+    boverBxgradH_vec : [boverBxgradH_x, boverBxgradH_z]
+  ) else if cdim = 1 then (
+    boverBxgradH_z : bioverJB_1_e*dH_dy - bioverJB_2_e*dH_dx,
+    boverBxgradH_vec : [boverBxgradH_z]
   ),
 
   /*Make sure to avoid having hamil[i]^2 or vmap[i]^2 in expressions*/
@@ -127,7 +171,8 @@ buildGKVolKernel(fh, funcNm, cdim, vdim, basisFun, polyOrder, varsInB, no_by) :=
   printf(fh, "  double vmap2 = vmap[1]*vmap[1]; ~%"),
   printf(fh, "~%"),
 
-  mvpar_e : dH_dz_e[vpardim]/vmap_prime_e[1],
+  mvpar_e : dH_dvpar,
+
   mvparsq_e : mvpar_e*mvpar_e/m_,
   isqlist : [],
   for i : 1 thru numP do (
@@ -144,86 +189,312 @@ buildGKVolKernel(fh, funcNm, cdim, vdim, basisFun, polyOrder, varsInB, no_by) :=
   ),
   printf(fh, "~%"),
 
+  /* Auxiliary variables to improve reading */
+  gradpsi : rdx2vec,
+  dpsidvpar : 1 / vmap_prime_e[1],
+
   /* Note: no contribution from mu. */
   for dir : 1 thru cdim+1 do (
 
-    dirLabel : varLabel[dir],
-
-    wDir      : eval_string(sconcat("w",dirLabel)),
-    rdDirVar2 : eval_string(sconcat("rd",dirLabel,"2")),
-  
-    dirVar    : varsP[dir],         /* Variable in current direction. */
+    dirLabel : varLabel[dir],  
 
     if dir = cdim then (
-      alpha_e : rtg33inv_e*dH_dz_e[vpardim]/vmap_prime_e[1]/m_
+      alpha_e : rtg33inv_e * dH_dvpar / m_ * gradpsi[dir] /* Contribution from B_0 . dH/dvpar . ∇ψ in z*/
     )
     else if dir = vpardim then (
-      alpha_e : -rtg33inv_e * dH_dz_e[cdim]/m_
+      alpha_e : -rtg33inv_e * gradH_vec[zidx]/m_ * dpsidvpar  /* Contribution from B_0 . ∇H . dψ/dvpar in z */
     )
     else (
-      alpha_e : 0
+      alpha_e : 0 /* all other B_0 contributions are 0 since B_0_x,y = 0 */
     ),
 
-
     if no_by = false then (
-      if cdim = 3 then (
-        curvdriftdir : dir
-      ),
-      if cdim = 2 then (
-        if dir = 1 then (
-          curvdriftdir : dir
-        ),
-        if dir = 2 then (
-          curvdriftdir : 3
-        )
-      ),
-      if cdim = 1 then (
-        curvdriftdir : 3
-      ),
 
+      /* Add curvature drift terms (m vpar/q * ∇ x b)/mB  = (m vpar * (∇ x b)/B)/qm */
       if dir < vpardim then (
-        alpha_e : alpha_e + dualcurlbhatoverB_list[curvdriftdir]*dH_dz_e[vpardim]/vmap_prime_e[1]*dH_dz_e[vpardim]/vmap_prime_e[1]/m_/q_
-      ),
-      if cdim = 3 then (
-        if dir = 1 then (
-          alpha_e : alpha_e + 1/q_ * (bioverJB_list[2]*dH_dz_e[3] - bioverJB_list[3]*dH_dz_e[2])
-        ),
-        if dir = 2 then (
-          alpha_e : alpha_e + 1/q_ * (bioverJB_list[3]*dH_dz_e[1] - bioverJB_list[1]*dH_dz_e[3])
-        ),
-        if dir = 3 then (
-          alpha_e : alpha_e + 1/q_ * (bioverJB_list[1]*dH_dz_e[2] - bioverJB_list[2]*dH_dz_e[1])
-        )
-      ),
-      if cdim = 2 then (
-        if dir = 1 then (
-          alpha_e : alpha_e + 1/q_ * (bioverJB_list[2]*dH_dz_e[2])
-        ),
-        if dir = 2 then (
-          alpha_e : alpha_e - 1/q_ * (bioverJB_list[2]*dH_dz_e[1])
-        )
-      ),
-      if dir = vpardim then (
-        if cdim = 3 then (
-          for k : 1 thru cdim do (
-            alpha_e : alpha_e - dualcurlbhatoverB_list[k]*dH_dz_e[k]*dH_dz_e[vpardim]/vmap_prime_e[1]/q_/m_
-          )
-        ),
-        if cdim = 2 then (
-            alpha_e : alpha_e - dualcurlbhatoverB_list[1]*dH_dz_e[1]*dH_dz_e[vpardim]/vmap_prime_e[1]/q_/m_ - dualcurlbhatoverB_list[3]*dH_dz_e[2]*dH_dz_e[vpardim]/vmap_prime_e[1]/q_/m_
-        ),
-        if cdim = 1 then (
-            alpha_e : alpha_e - dualcurlbhatoverB_list[3]*dH_dz_e[1]*dH_dz_e[vpardim]/vmap_prime_e[1]/q_/m_
-        )
+        /* 1/m vpar/q [(∇ x b)/B] dH/dvpar . ∇ψ */
+        alpha_e : alpha_e + 1/m_ * mvpar_e/q_ * dualcurlbhatoverB_vec[dir] * dH_dvpar * gradpsi[dir],
+        /* 1/q (b/B x ∇H) . ∇ψ */
+        alpha_e : alpha_e + 1/q_ * boverBxgradH_vec[dir] * gradpsi[dir]
+
+      ) else (
+        /* 1/m (m vpar/q [(∇ x b)/B]) . ∇H dψ/dvpar */
+        for k : 1 thru cdim do (
+          alpha_e : alpha_e - 1/m_ * mvpar_e/q_ * dualcurlbhatoverB_vec[k] * gradH_vec[k] * dpsidvpar
+        ) 
       )
-
     ),
 
-    if dir < vpardim then (
-      alpha_e : alpha_e*rdx2vec[dir]
+    /* Project alpha on basis and write to array. */
+    printf(fh, "  double alpha~a[~a] = {0.}; ~%", dirLabel, numP),
+    alpha_c : fullratsimp(calcInnerProdList(varsP, 1, bP, alpha_e)),
+    alpha_c : subst(replaceList, alpha_c),
+    alpha_c : subst(replaceListHamil, alpha_c),
+    alpha_c : subst(replaceListVpar, alpha_c),
+    alpha_c : subst(dvparSimp, alpha_c),
+    alphaLabel : eval_string(sconcat(alpha, dirLabel)),
+    clst : [rdx2vec, rdv2vec, m_, q_, wvpar, rdvpar2Sq,
+            makelist(dxv[i-1],i,1,pdim), makelist(vmap[i-1],i,1,2*length(vmap_e[1]))],
+    writeCExprsCollect1(alphaLabel, alpha_c, clst),
+    printf(fh, "~%"),
+    flush_output(fh),
+    alphaNoZero_c : makelistNoZeros1(alpha_c, alphaLabel),
+    alpha_e       : doExpand(alphaNoZero_c, bP),
+
+    alphaJf_e : alpha_e*Jf_e,
+
+    printf(fh, "~%"),
+    volTerm_c : fullratsimp(calcInnerProdList(varsP, 1, diff(bP,varsP[dir]), alphaJf_e)),
+    volTerm_c : subst(replaceList, volTerm_c),
+    writeCIncrExprsNoExpand(gcfac(float(expand(volTerm_c)))),
+    flush_output(fh),
+    printf(fh, "~%")
+
+  ),
+
+  printf(fh, "  return 0.; ~%"),
+  printf(fh, "} ~%")
+
+)$
+
+addAparGKEMVolKernel(fh, funcNm, cdim, vdim, basisFun, polyOrder, varsInB, no_by) := block(
+  [pdim,varsC,bC,varsP,bP,vSub,numC,numP,varLabel,d,rdx2vec,rdv2vec,allVarLabelsC,
+   bmagBasis,b_1_e,b_2_e,b_3_e,jacobTotInv_e,vmap_e,hamil_e,dir,dirLabel,gradpsi, dpsidvpar,
+   alpha_e,alpha_c,alphaLabel,alphaNoZero_c,volTerm_c, alphaJf_e, Jf_e, 
+   replaceListHamil, replaceListVpar,isqlist,mvpar_e,mvparsq_e, xidx, yidx, zidx, vidx,
+   replaceList,dvparSimp,phi_e,bmag_e,rtg33inv_e,bioverJB_1_e,bioverJB_2_e,bioverJB_3_e,bioverJB_vec,
+   dualcurlbhatoverB_1_e,dualcurlbhatoverB_2_e,dualcurlbhatoverB_3_e,dualcurlbhatoverB_vec, dH_dx, dH_dy, dH_dz, dH_dvpar, gradH_vec,
+   vmapSq_e,vmap_prime_e,hamil_c,hamilCvar,hamilNoZero_c,vpardim,i,k,Apar_e, dA_dx, dA_dy, dA_dz, gradA_vec,
+   rotAbovermB_x,rotAbovermB_y,rotAbovermB_z,rotAbovermB_vec, gradAxbhatoverB_x,gradAxbhatoverB_y,gradAxbhatoverB_z,gradAxboverB_vec,
+   clst],
+
+  kill(varsC,varsP,bC,bP),
+  pdim : cdim+vdim,
+
+  [varsC,bC,varsP,bP,vSub] : loadGkBasis(basisFun, cdim, vdim, polyOrder),
+  numC : length(bC),  numP : length(bP), 
+
+  varLabel : makelist(string(varsP[d]),d,1,pdim),
+
+  print("Working on ", funcNm),
+  printf(fh, "GKYL_CU_DH double ~a(const double *w, const double *dxv, const double *vmap, const double *vmapSq,
+    const double q_, const double m_, const double *bmag, const double *jacobtot_inv, const double *dualcurlbhatoverB, const double *bioverJB,
+    const double *b_i, const double *phi, const double *apar, const double *fin, double* GKYL_RESTRICT out) ~%{ ~%", funcNm),
+  printf(fh, "  // w[NDIM]: cell-center.~%"),
+  printf(fh, "  // dxv[NDIM]: cell length.~%"),
+  printf(fh, "  // vmap: velocity space mapping.~%"),
+  printf(fh, "  // vmapSq: velocity space mapping squared.~%"),
+  printf(fh, "  // q_,m_: species charge and mass.~%"),
+  printf(fh, "  // bmag: magnetic field amplitude.~%"),
+  printf(fh, "  // jacobtot_inv: reciprocal of the conf-space jacobian time the guiding center coordinate Jacobian.~%"),
+  printf(fh, "  // dualcurlbhatoverB: dual curl of bhat over B.~%"),
+  printf(fh, "  // bioverJB: b_i over J times B.~%"),
+  printf(fh, "  // b_i: covariant components of the field aligned unit vector.~%"),
+  printf(fh, "  // apar: parallel component of magnetic vector potential.~%"),
+  printf(fh, "  // phi: electrostatic potential .~%"),
+  printf(fh, "  // fin: Distribution function.~%"),
+  printf(fh, "  // out: output increment.~%"),
+  printf(fh, "~%"),
+
+  /* Declare cell-center variables and variables multiplying gradients. */
+  for d : 1 thru pdim do (
+    printf(fh, "  double rd~a2 = 2.0/dxv[~a];~%", varLabel[d], d-1)
+  ),
+  printf(fh, "~%"),
+  rdx2vec : makelist(eval_string(sconcat("rd",varLabel[i],"2")),i,1,cdim),
+  rdv2vec : makelist(eval_string(sconcat("rd",varLabel[i],"2")),i,cdim+1,pdim),
+
+  /* Declare variables with squared of cell centers and rdx2 variables (only need vpar^2). */
+  /* printf(fh, "  double rdvpar2Sq = rdvpar2*rdvpar2;~%"), */
+  /* printf(fh, "  double dvparSq = dxv[~a]*dxv[~a];~%", cdim, cdim), */
+  printf(fh, "~%"),
+  replaceList : [rdx^2=dx2Sq, rdy^2=dy2Sq, rdz^2=dx2Sq, rdvpar2^2=rdvpar2Sq,dxv[cdim]^2=dvparSq,rdvpar2Sq=4/dvparSq],
+  dvparSimp : append(makelist(dxv[i-1]=2/eval_string(sconcat("rd",varLabel[i],"2")),i,1,pdim),
+                     [dvparSq=4/rdvpar2Sq]),
+
+  /* Store indices of the coordinate directions. */
+  if cdim = 3 then (
+    xidx : 1, yidx : 2, zidx : 3, vidx : 4
+  ) else if cdim = 2 then (
+    xidx : 1, zidx : 2, vidx : 3
+  ) else if cdim = 1 then (
+    zidx : 1, vidx : 2
+  ),
+
+  /* Create pointers to the components of b_i. */
+  allVarLabelsC : ["x","y","z"],
+  for d : 1 thru 3 do (
+    printf(fh, "  const double *b_~a = &b_i[~a];~%", d, numC*(d-1))
+  ),
+  printf(fh, "~%"),
+
+  for d : 1 thru 3 do (
+    printf(fh, "  const double *bioverJB_~a = &bioverJB[~a]; ~%", d, numC*(d-1))
+  ),
+  printf(fh, "~%"),
+
+  /* Create pointers to the components of dualcurlbhatoverB. */
+  for d : 1 thru 3 do (
+    printf(fh, "  const double *dualcurlbhatoverB_~a = &dualcurlbhatoverB[~a]; ~%", d, numC*(d-1))
+  ),
+  printf(fh, "~%"),
+
+  /* Axisymmetric basis (independent of y). */
+  bmagBasis : getAxisymmetricConfBasis(bC),
+  /* Expand input fields for Hamiltonian calculation */
+  phi_e     : doExpand1(phi,bC),
+  bmag_e : doExpand1(bmag, bmagBasis),
+  rtg33inv_e : doExpand1(rtg33inv, bmagBasis),
+  dualcurlbhatoverB_1_e : doExpand1(dualcurlbhatoverB_1, bmagBasis),
+  dualcurlbhatoverB_2_e : doExpand1(dualcurlbhatoverB_2, bmagBasis),
+  dualcurlbhatoverB_3_e : doExpand1(dualcurlbhatoverB_3, bmagBasis),
+  bioverJB_1_e : doExpand1(bioverJB_1, bmagBasis),
+  bioverJB_2_e : doExpand1(bioverJB_2, bmagBasis),
+  bioverJB_3_e : doExpand1(bioverJB_3, bmagBasis),
+  b_1_e : doExpand1(b_x, bmagBasis),
+  b_2_e : doExpand1(b_y, bmagBasis),
+  if (no_by or cdim = 1) then (b_2_e : 0),
+  b_3_e : doExpand1(b_z, bmagBasis),
+  jacobTotInv_e : doExpand1(jacobtot_inv, bmagBasis),
+
+  if cdim = 3 then (
+    dualcurlbhatoverB_vec : [dualcurlbhatoverB_1_e, dualcurlbhatoverB_2_e, dualcurlbhatoverB_3_e],
+    bioverJB_vec : [bioverJB_1_e, bioverJB_2_e, bioverJB_3_e]
+  ) else if cdim = 2 then (
+    dualcurlbhatoverB_vec : [dualcurlbhatoverB_1_e, dualcurlbhatoverB_3_e],
+    bioverJB_vec : [bioverJB_1_e, bioverJB_3_e]
+  ) else if cdim = 1 then (
+    dualcurlbhatoverB_vec : [dualcurlbhatoverB_3_e],
+    bioverJB_vec : [bioverJB_3_e]
+  ),
+
+  /* Velocity mapping fields. */
+  [vmap_e,vmapSq_e,vmap_prime_e] : expandVmapFields(varsP),
+
+  /* Redefine vmap_prime to exploit the relationship between it and vmap. */
+  vmap_prime_e : makelist(diff(vmap_e[d],varsP[cdim+d]),d,1,vdim),
+
+  /* Write out the hamiltonian*/
+  hamil_e : q_*phi_e + (1/2)*m_*vmapSq_e[1],
+  if vdim > 1 then ( hamil_e : hamil_e + vmap_e[2]*bmag_e ),
+  hamil_c : calcInnerProdList(varsP, 1, bP, hamil_e),
+  printf(fh, "  double hamil[~a] = {0.}; ~%", numP),
+  replaceList : [wvpar^2=wvparSq, rdvpar2^2=rdvpar2Sq, rdx2^2=rdx2Sq, rdy2^2=rdy2Sq, rdz2^2=rdz2Sq, m_^2=mSq, q_^2=qSq],
+  hamilCvar : eval_string(sconcat("hamil")),
+  writeCExprsNoExpand1(hamilCvar, gcfac(float(expand(subst(replaceList, hamil_c))))),
+  printf(fh, "~%"),
+  flush_output(fh),
+  hamilNoZero_c : makelistNoZeros1(hamil_c, hamilCvar),
+  /* Expand projected Hamiltonian on basis. */
+  hamil_e : hamilNoZero_c . bP,
+
+  /*Expand Jf*/
+  Jf_e : doExpand1(fin,bP),
+
+  /* Calculate expressions for dericatives of the hamiltonian*/
+  vpardim : pdim-1,
+  if vdim = 1 then ( vpardim : pdim ),
+  dH_dvpar : diff(hamil_e,varsP[vpardim])/vmap_prime_e[1], /* We do not multply by rdv2vec because it is inside vmap_prime */
+  if cdim = 3 then (
+    dH_dx : diff(hamil_e*rdx2vec[xidx],varsP[xidx]),
+    dH_dy : diff(hamil_e*rdx2vec[yidx],varsP[yidx]),
+    dH_dz : diff(hamil_e*rdx2vec[zidx],varsP[zidx]),
+    gradH_vec : [dH_dx, dH_dy, dH_dz]
+  ) else if cdim = 2 then (
+    dH_dx : diff(hamil_e*rdx2vec[xidx],varsP[xidx]),
+    dH_dz : diff(hamil_e*rdx2vec[zidx],varsP[zidx]),
+    gradH_vec : [dH_dx, dH_dz]
+  ) else if cdim = 1 then (
+    dH_dz : diff(hamil_e*rdx2vec[zidx],varsP[zidx]),
+    gradH_vec : [dH_dz]
+  ),
+
+  /*Make sure to avoid having hamil[i]^2 or vmap[i]^2 in expressions*/
+  replaceListVpar : [vmap[1]^2=vmap2],
+  printf(fh, "  double vmap2 = vmap[1]*vmap[1]; ~%"),
+  printf(fh, "~%"),
+
+  mvpar_e : dH_dvpar,
+  mvparsq_e : mvpar_e*mvpar_e/m_,
+  isqlist : [],
+  for i : 1 thru numP do (
+    if freeof(hamil[i]^2, expand(mvparsq_e)) = false then (
+      isqlist : append(isqlist,[i])
     )
-    else if dir = vpardim then (
-      alpha_e : alpha_e/vmap_prime_e[1]
+  ),
+
+  replaceListHamil : [],
+  printf(fh, "  double hamil2[~a] = {0.}; ~%", length(isqlist)),
+  for i : 1 thru length(isqlist) do (
+    printf(fh, "  hamil2[~a] = hamil[~a]*hamil[~a]; ~%", i-1, isqlist[i], isqlist[i]),
+    replaceListHamil : append(replaceListHamil, [hamil[isqlist[i]]^2=hamil2[i-1]])
+  ),
+  printf(fh, "~%"),
+
+  /* Expand Apar.*/
+  /* Apar_e : subst(y=0, doExpand1(apar,bC)), */
+  Apar_e : doExpand1(apar,bC),
+
+  /* Expand dBperp/Bmag using product rule method Apar ∇ x bhat + ∇Apar x bhat. */
+  /* Calculate ∇Apar in a vector */
+  if cdim = 3 then (
+    dA_dx : diff(Apar_e*rdx2vec[xidx],varsP[xidx]),
+    dA_dy : diff(Apar_e*rdx2vec[yidx],varsP[yidx]),
+    dA_dz : diff(Apar_e*rdx2vec[zidx],varsP[zidx]),
+    gradA_vec : [dA_dx, dA_dy, dA_dz]
+  ) else if cdim = 2 then (
+    dA_dx : diff(Apar_e*rdx2vec[xidx],varsP[xidx]),
+    dA_dy : 0,
+    dA_dz : diff(Apar_e*rdx2vec[zidx],varsP[zidx]),
+    gradA_vec : [dA_dx, dA_dz]
+  ) else if cdim = 1 then (
+    dA_dx : 0,
+    dA_dy : 0,
+    dA_dz : diff(Apar_e*rdx2vec[zidx],varsP[zidx]),
+    gradA_vec : [dA_dz]
+  ),
+
+  /* Use bioverJB to calculate ∇Apar x b / B in a vector */
+  if cdim = 3 then (
+    gradAxbhatoverB_x : dA_dy*bioverJB_3_e - dA_dz*bioverJB_2_e,
+    gradAxbhatoverB_y : dA_dz*bioverJB_1_e - dA_dx*bioverJB_3_e,
+    gradAxbhatoverB_z : dA_dx*bioverJB_2_e - dA_dy*bioverJB_1_e,
+    gradAxboverB_vec : [gradAxbhatoverB_x, gradAxbhatoverB_y, gradAxbhatoverB_z]
+  ) else if cdim = 2 then (
+    gradAxbhatoverB_x : dA_dy*bioverJB_3_e - dA_dz*bioverJB_2_e,
+    gradAxbhatoverB_z : dA_dx*bioverJB_2_e - dA_dy*bioverJB_1_e,
+    gradAxboverB_vec : [gradAxbhatoverB_x, gradAxbhatoverB_z]
+  ) else if cdim = 1 then (
+    gradAxbhatoverB_z : dA_dx*bioverJB_2_e - dA_dy*bioverJB_1_e,
+    gradAxboverB_vec : [gradAxbhatoverB_z]
+  ),
+  
+  /* Auxiliary variables to improve reading */
+  gradpsi : rdx2vec,
+  dpsidvpar : 1 / vmap_prime_e[1],
+
+  /* Note: no contribution from mu. */
+  for dir : 1 thru cdim+1 do (
+
+    dirLabel : varLabel[dir],
+    alpha_e : 0,
+
+    if no_by = false then (
+      if dir < vpardim then ( /* Config space contributions ( curv drift * dH/dvpar . ∇ψ ) */
+        /* 1/m Apar [(∇ x b)/B] dH/dvpar . ∇ψ */
+        alpha_e : alpha_e + 1/m_ * Apar_e * dualcurlbhatoverB_vec[dir] * dH_dvpar * gradpsi[dir],
+        /* 1/m [∇Apar x b/B] dH/dvpar . ∇ψ */
+        alpha_e : alpha_e + 1/m_ * gradAxboverB_vec[dir] * dH_dvpar * gradpsi[dir]
+
+      ) else ( /* Vpar contribution ( curv drift . ∇H * dψ/dvpar )*/
+        for k : 1 thru cdim do (
+          /* 1/m [Apar . (∇ x b)/B] . ∇H * dψ/dvpar */
+          alpha_e : alpha_e - 1/m_ * Apar_e * dualcurlbhatoverB_vec[k] * gradH_vec[k] * dpsidvpar,
+          /* 1/m [∇Apar x b/B] . ∇H * dψ/dvpar */
+          alpha_e : alpha_e - 1/m_ * gradAxboverB_vec[k] * gradH_vec[k] * dpsidvpar
+
+        )
+      )
     ),
 
     /* Project alpha on basis and write to array. */
@@ -235,7 +506,7 @@ buildGKVolKernel(fh, funcNm, cdim, vdim, basisFun, polyOrder, varsInB, no_by) :=
     alpha_c : subst(dvparSimp, alpha_c),
     alphaLabel : eval_string(sconcat(alpha, dirLabel)),
     clst : [rdx2vec, rdv2vec, m_, q_, wvpar, rdvpar2Sq,
-            makelist(dxv[i-1],i,1,pDim), makelist(vmap[i-1],i,1,2*length(vmap_e[1]))],
+            makelist(dxv[i-1],i,1,pdim), makelist(vmap[i-1],i,1,2*length(vmap_e[1]))],
     writeCExprsCollect1(alphaLabel, alpha_c, clst),
     printf(fh, "~%"),
     flush_output(fh),
@@ -257,3 +528,70 @@ buildGKVolKernel(fh, funcNm, cdim, vdim, basisFun, polyOrder, varsInB, no_by) :=
   printf(fh, "} ~%")
 
 )$
+
+addApardotGKEMVolKernel(fh, funcNm, cdim, vdim, basisFun, polyOrder, varsInB, no_by) := block(
+  [pdim,varsC,bC,varsP,bP,vSub,numC,numP,varLabel,d,rdx2vec,rdv2vec, vmap_e,vpardir,dirLabel,
+   alpha_e,alpha_c,alphaLabel,alphaNoZero_c,volTerm_c, alphaJf_e, Jf_e,vmapSq_e,vmap_prime_e,
+   apardot_e,clst],
+
+  kill(varsC,varsP,bC,bP),
+  pdim : cdim+vdim,
+
+  [varsC,bC,varsP,bP,vSub] : loadGkBasis(basisFun, cdim, vdim, polyOrder),
+  numC : length(bC),  numP : length(bP), 
+
+  varLabel : makelist(string(varsP[d]),d,1,pdim),
+
+  print("Working on ", funcNm),
+  printf(fh, "GKYL_CU_DH double ~a(const double *vmap, const double q_, const double m_, const double *apardot,
+    const double *fin, double* GKYL_RESTRICT out) ~%{ ~%", funcNm),
+  printf(fh, "  // q_,m_: species charge and mass.~%"),
+  printf(fh, "  // apardot: time derivative of parallel component of magnetic vector potential.~%"),
+  printf(fh, "  // fin: Distribution function.~%"),
+  printf(fh, "  // out: output increment.~%"),
+  printf(fh, "~%"),
+
+  rdx2vec : makelist(eval_string(sconcat("rd",varLabel[i],"2")),i,1,cdim),
+  rdv2vec : makelist(eval_string(sconcat("rd",varLabel[i],"2")),i,cdim+1,pdim),
+
+
+  /* Velocity mapping fields. */
+  [vmap_e,vmapSq_e,vmap_prime_e] : expandVmapFields(varsP),
+  vmap_prime_e : makelist(diff(vmap_e[d],varsP[cdim+d]),d,1,vdim),
+
+  /*Expand Jf*/
+  Jf_e : doExpand1(fin,bP),
+
+  /* Expand Apardot */
+  apardot_e : doExpand1(apardot,bC),
+
+  /* Note: only a vpar contribution. */
+  vpardir : cdim+1,
+  dirLabel : varLabel[vpardir],
+
+  alpha_e : -q_/m_ * apardot_e /vmap_prime_e[1],
+
+  /* Project alpha on basis and write to array. */
+  printf(fh, "  double alpha~a[~a] = {0.}; ~%", dirLabel, numP),
+  alpha_c : fullratsimp(calcInnerProdList(varsP, 1, bP, alpha_e)),
+  alphaLabel : eval_string(sconcat(alpha, dirLabel)),
+  clst : [rdx2vec, rdv2vec, m_, q_, wvpar, rdvpar2Sq,
+          makelist(dxv[i-1],i,1,pdim), makelist(vmap[i-1],i,1,2*length(vmap_e[1]))],
+  writeCExprsCollect1(alphaLabel, alpha_c, clst),
+  printf(fh, "~%"),
+  flush_output(fh),
+  alphaNoZero_c : makelistNoZeros1(alpha_c, alphaLabel),
+  alpha_e       : doExpand(alphaNoZero_c, bP),
+
+  alphaJf_e : alpha_e*Jf_e,
+
+  printf(fh, "~%"),
+  volTerm_c : fullratsimp(calcInnerProdList(varsP, 1, diff(bP,varsP[vpardir]), alphaJf_e)),
+  writeCIncrExprsNoExpand(gcfac(float(expand(volTerm_c)))),
+  flush_output(fh),
+  printf(fh, "~%"),
+
+  printf(fh, "  return 0.; ~%"),
+  printf(fh, "} ~%")
+
+)$
\ No newline at end of file
diff --git a/maxima/g0/gk_collisionless/gk_collisionless_flux-surf-conf.mac b/maxima/g0/gk_collisionless/gk_collisionless_flux-surf-conf.mac
index 34e06ea7..ac9faffc 100644
--- a/maxima/g0/gk_collisionless/gk_collisionless_flux-surf-conf.mac
+++ b/maxima/g0/gk_collisionless/gk_collisionless_flux-surf-conf.mac
@@ -6,13 +6,15 @@ load("utilities_gyrokinetic")$
 load("nodal_operations/nodal_functions")$
 fpprec : 24$
 
-buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no_by, edge, mb_bound) := block(
+buildGKFluxConfKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no_by, edge, mb_bound, scheme) := block(
   [pDim,varsC,bC,varsP,bP,vSub,numC,numP,surfVar,varLabel,dirLabel,surfIntVars,surf_cvars,surf_vvars,
    surfNodes,nodeVars,bSurf,basisNodal,surfConfigNodes,numSurfNodes,numSurfConfigNodes,numVelNodes,
-   numMuNodes,numVparNodes,d,rdx2vec,rdv2vec,rdSurfVar2,bmagBasis,phi_e,bmagSurf_e,vmap_e,vmapSq_e,
+   numMuNodes,numVparNodes,d,rdx2vec,rdv2vec,rdSurfVar2,bmagBasis,phi_e,apar_e,
+   bmagSurf_e,vmap_e,vmapSq_e,
    vmap_prime_e,evPoint,hamil_e,hamil_c,replaceList,hamilNoZero_c,JfL_e,JfR_e,JfL_c,JfR_c,
    jacobgeo_rat_surfR_e,jacobgeo_rat_surfL_e,JfL_nodes,JfR_nodes,vmap_prime_nodes,vpardim,
-   dH_dz_nodes,mvpar_nodes,di3,i,j,j0index,j1index,vparindex,vpar0index,pOrderCFL
+   dH_dz_nodes,mvpar_nodes,i,j,j0index,j1index,vparindex,vpar0index,pOrderCFL,
+   surfIntVarsC,bSurfC,hamilCvar,aparCvar,apar_nodes,apar_c,dA_dx_nodes
    ],
 
   kill(varsC,varsP,bC,bP),
@@ -60,7 +62,8 @@ buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
     const double *w, const double *dxv,
     const double *vmap, const double *vmapSq, const double q_, const double m_,
     const struct gkyl_dg_surf_geom *dgs, const struct gkyl_gk_dg_surf_geom *gkdgs, 
-    const double *bmag, const double *jacobgeo_rat_surfL, const double *jacobgeo_rat_surfR, const double *phi,
+    const double *bmag, const double *jacobgeo_rat_surfL, const double *jacobgeo_rat_surfR, 
+    const double *phi, const double *apar, const double *apardot,
     const double *JfL, const double *JfR, double* GKYL_RESTRICT flux_surf) ~%{ ~%", funcNm),
   printf(fh, "  // w[NDIM]: cell-center.~%"),
   printf(fh, "  // dxv[NDIM]: cell length.~%"),
@@ -73,6 +76,8 @@ buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
   printf(fh, "  // jacobgeo_rat_surfL: Ratio of surface conf-space Jacobians in left cell.~%"),
   printf(fh, "  // jacobgeo_rat_surfR: Ratio of surface conf-space Jacobians in right cell.~%"),
   printf(fh, "  // phi: electrostatic potential.~%"),
+  printf(fh, "  // apar: parallel component of vector potential.~%"),
+  printf(fh, "  // apardot: time derivative of parallel component of vector potential.~%"),
   printf(fh, "  // JfL: distribution times total jacobian in left cell.~%"),
   printf(fh, "  // JfR: distribution times total jacobian in right cell.~%"),
   printf(fh, "  // flux_surf: output surface phase space flux in each direction (cdim + 1 components).~%"),
@@ -120,7 +125,7 @@ buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
   flush_output(fh),
   hamilNoZero_c : makelistNoZeros1(hamil_c, hamilCvar),
   /* Expand projected Hamiltonian on basis. */
-  hamil_e : hamilNoZero_c . bSurf,
+  hamil_e : doExpand(hamilNoZero_c, bSurf),
 
   /* fl and fr */
   JfL_e : doExpand1(JfL, bP),
@@ -173,16 +178,28 @@ buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
     )
   ),
 
-  mvpar_nodes : [],
+  dHdvpar_nodes : [],
   for i : 1 thru numVparNodes do (
-    mvpar_nodes : append(mvpar_nodes, [dH_dz_nodes[vpardim][i]])
+    dHdvpar_nodes : append(dHdvpar_nodes, [dH_dz_nodes[vpardim][i]])
   ),
 
-  if surfDir = cdim then(
-    di3 : true
-  )
-   else (
-    di3 : false
+  /* Expand Aparallel. */
+  apar_e : doExpand1(apar,bC),
+  apar_c : calcInnerProdList(surfIntVars, 1, bSurf, subst(surfVar=evPoint,apar_e)),
+  printf(fh, "  double apar_surf[~a] = {0.}; ~%", length(bSurf)),
+  aparSurfCvar : eval_string(sconcat("apar_surf")),
+  writeCExprsNoExpand1(aparSurfCvar, gcfac(float(expand(apar_c)))),
+  printf(fh, "~%"),
+  flush_output(fh),
+  apar_c : makelistNoZeros1(apar_c, aparSurfCvar),
+  /* Expand projected Apar on basis. */
+  apar_e : doExpand(apar_c, bSurf),
+  /* Eval Aparallel at nodes. */
+  apar_nodes : float(evAtNodes(apar_e,surfNodes,surfIntVars)),
+  /* Compute gradient of Aparallel */
+  dA_dx_nodes : makelist(0, i, 1, cdim),
+  for i : 1 thru cdim do (
+    dA_dx_nodes[i] : float(evAtNodes(diff(apar_e*rdx2vec[i],varsP[i]),surfNodes,surfIntVars))
   ),
 
   /* Now calculate flux at all quadrature nodes */
@@ -202,14 +219,14 @@ buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
   printf(fh, "  double Jfavg_quad = 0.0; ~%"),
   printf(fh, "  double Jfjump_quad = 0.0; ~%"),
 
-  printf(fh, "  double mvpar_quad[3] = {0.0}; ~%"),
+  printf(fh, "  double dHdvpar_quad[3] = {0.0}; ~%"),
   for i : 1 thru numVparNodes do (
-    printf(fh, "  mvpar_quad[~a] = ~a; ~%", i-1, mvpar_nodes[i])
+    printf(fh, "  dHdvpar_quad[~a] = ~a; ~%", i-1, dHdvpar_nodes[i])
   ),
 
   printf(fh, "  double mvparsq_quad[3] = {0.0}; ~%"),
   for i : 1 thru numVparNodes do (
-    printf(fh, "  mvparsq_quad[~a] = mvpar_quad[~a]*mvpar_quad[~a]/m_; ~%", i-1, i-1,i-1)
+    printf(fh, "  mvparsq_quad[~a] = dHdvpar_quad[~a]*dHdvpar_quad[~a]/m_; ~%", i-1, i-1,i-1)
   ),
   printf(fh, "~%"),
 
@@ -229,64 +246,91 @@ buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
       vparindex : mod(j-1, numVparNodes) + 1,
       vpar0index : mod(j-1, numVparNodes),
       printf(fh, "~%"),
-      if no_by = true then (
-        if di3 = true then (
-          printf(fh, "  alpha_quad = (mvpar_quad[~a]*B3_quad/(m_*bmag_quad))*area_elem_quad/Jc_quad; ~%", vpar0index)
-        )
-        else (
-          printf(fh, "  alpha_quad = 0.0; ~%")
-        )
+
+      /* Electrostatic term */
+      if surfDir = cdim then (
+        /* Parallel streaming term */
+        printf(fh, "  alpha_quad = (dHdvpar_quad[~a]*B3_quad/(m_*bmag_quad)); ~%", vpar0index)
+      ) else (
+        printf(fh, "  alpha_quad = 0.0; ~%")
       ),
+
       if no_by = false then (
-        /*printf(fh, "  alpha_quad += mvparsq_quad[~a]*normcurlbhat_quad/(bmag_quad*q_) ;~%", vpar0index),*/
+        printf(fh, "  alpha_quad += mvparsq_quad[~a]*normcurlbhat_quad/(bmag_quad*q_); ~%", vpar0index),
+        /* + 1/q b/B x ∇H . ∇ψ*/
         if cdim = 3 then (
            if surfDir = 1 then(
-             printf(fh, "  alpha_quad = (mvparsq_quad[~a]*normcurlbhat_quad/(bmag_quad*q_) + 1/(q_*bmag_quad*area_elem_quad) * (bhat_quad[1]*(~a) - bhat_quad[2]*(~a)))*area_elem_quad/Jc_quad; ~%", vpar0index, dH_dz_nodes[3][j1index], dH_dz_nodes[2][j1index])
+             printf(fh, "  alpha_quad += 1/(q_*bmag_quad*area_elem_quad) * (bhat_quad[1]*(~a) - bhat_quad[2]*(~a)); ~%", dH_dz_nodes[3][j1index], dH_dz_nodes[2][j1index])
            ),
            if surfDir = 2 then(
-             printf(fh, "  alpha_quad = (mvparsq_quad[~a]*normcurlbhat_quad/(bmag_quad*q_) + 1/(q_*bmag_quad*area_elem_quad) * (bhat_quad[2]*(~a) - bhat_quad[0]*(~a)))*area_elem_quad/Jc_quad; ~%", vpar0index, dH_dz_nodes[1][j1index], dH_dz_nodes[3][j1index])
+             printf(fh, "  alpha_quad += 1/(q_*bmag_quad*area_elem_quad) * (bhat_quad[2]*(~a) - bhat_quad[0]*(~a)); ~%", dH_dz_nodes[1][j1index], dH_dz_nodes[3][j1index])
            ),
            if surfDir = 3 then(
-             printf(fh, "  alpha_quad = (mvpar_quad[~a]*B3_quad/(m_*bmag_quad) + mvparsq_quad[~a]*normcurlbhat_quad/(bmag_quad*q_) + 1/(q_*bmag_quad*area_elem_quad) * (bhat_quad[0]*(~a) - bhat_quad[1]*(~a)))*area_elem_quad/Jc_quad; ~%", vpar0index, vpar0index, dH_dz_nodes[2][j1index], dH_dz_nodes[1][j1index])
+             printf(fh, "  alpha_quad += 1/(q_*bmag_quad*area_elem_quad) * (bhat_quad[0]*(~a) - bhat_quad[1]*(~a)); ~%", dH_dz_nodes[2][j1index], dH_dz_nodes[1][j1index])
            )
         ),
         if cdim = 2 then (
            if surfDir = 1 then(
-             printf(fh, "  alpha_quad = (mvparsq_quad[~a]*normcurlbhat_quad/(bmag_quad*q_) + 1/(q_*bmag_quad*area_elem_quad) * bhat_quad[1]*(~a))*area_elem_quad/Jc_quad; ~%", vpar0index, dH_dz_nodes[2][j1index])
+             printf(fh, "  alpha_quad += 1/(q_*bmag_quad*area_elem_quad) * bhat_quad[1]*(~a); ~%", dH_dz_nodes[2][j1index])
            ),
            if surfDir = 2 then(
-             printf(fh, "  alpha_quad = (mvpar_quad[~a]*B3_quad/(m_*bmag_quad) + mvparsq_quad[~a]*normcurlbhat_quad/(bmag_quad*q_) + 1/(q_*bmag_quad*area_elem_quad) * -bhat_quad[1]*(~a))*area_elem_quad/Jc_quad;~%", vpar0index, vpar0index, dH_dz_nodes[1][j1index])
+             printf(fh, "  alpha_quad -= 1/(q_*bmag_quad*area_elem_quad) * bhat_quad[1]*(~a); ~%", dH_dz_nodes[1][j1index])
            )
         ),
-        if cdim = 1 then (
-          printf(fh, "  alpha_quad = (mvpar_quad[~a]*B3_quad/(m_*bmag_quad))*area_elem_quad/Jc_quad; ~%", vpar0index)
+
+        /* Electromagnetic Apar contribution using ∇ x (A b) = A ∇ x b + ∇A x b */
+        /* No contribution for cdim = 1 */
+        if cdim > 1 then (
+          /* + 1/m A (∇ x b)/B . ∇ψ dH/dvpar */
+          printf(fh, "  alpha_quad += 1/m_ * (~a) * normcurlbhat_quad/bmag_quad * dHdvpar_quad[~a]; ~%", apar_nodes[j1index], vpar0index)
+        ),
+        /* + 1/m ∇A x b/B . ∇ψ dH/dvpar */
+        if cdim = 3 then (
+          if surfDir = 1 then(
+            printf(fh, "  alpha_quad += 1/(m_*bmag_quad*area_elem_quad) * ((~a) * bhat_quad[2] - (~a) * bhat_quad[1]) * dHdvpar_quad[~a]; ~%", dA_dx_nodes[2][j1index], dA_dx_nodes[3][j1index], vpar0index)
+          ),
+          if surfDir = 2 then(
+            printf(fh, "  alpha_quad += 1/(m_*bmag_quad*area_elem_quad) * ((~a) * bhat_quad[0] - (~a) * bhat_quad[2]) * dHdvpar_quad[~a]; ~%", dA_dx_nodes[3][j1index], dA_dx_nodes[1][j1index], vpar0index)
+          ),
+          if surfDir = 3 then(
+            printf(fh, "  alpha_quad += 1/(m_*bmag_quad*area_elem_quad) * ((~a) * bhat_quad[1] - (~a) * bhat_quad[0]) * dHdvpar_quad[~a]; ~%", dA_dx_nodes[1][j1index], dA_dx_nodes[2][j1index], vpar0index)
+          )
+        ),
+        /* in 2D, we take the first and last component of the 3D cross product, setting d/dy = 0.*/
+        if cdim = 2 then (
+          if surfDir = 1 then(
+            printf(fh, "  alpha_quad -= 1/(m_*bmag_quad*area_elem_quad) * (~a) * bhat_quad[1] * dHdvpar_quad[~a]; ~%", dA_dx_nodes[2][j1index], vpar0index)
+          ),
+          if surfDir = 2 then(
+            printf(fh, "  alpha_quad += 1/(m_*bmag_quad*area_elem_quad) * (~a) * bhat_quad[1] * dHdvpar_quad[~a]; ~%", dA_dx_nodes[1][j1index], vpar0index)
+          )
         )
+        /* Nothing in 1D */
       ),
 
+      /* Multiply by |e^i|, note: area_elem_quad = J_c |e^i| */
+      printf(fh, "  alpha_quad = alpha_quad * area_elem_quad/Jc_quad; ~%"),
+
       printf(fh, "~%"),
-      /*printf(fh, "  alpha_quad = alpha_quad*area_elem_quad/Jc_quad; ~%"),*/
       printf(fh, "  cfl = fmax(fabs(alpha_quad), fabs(cfl)); ~%"),
       printf(fh, "  JfL_quad = ~a; ~%", JfL_nodes[j1index]),
       printf(fh, "  JfR_quad = ~a; ~%", JfR_nodes[j1index]),
       printf(fh, "  Jfavg_quad = (JfL_quad + JfR_quad)/2.0; ~%"),
       printf(fh, "  Jfjump_quad = (JfR_quad - JfL_quad)/2.0; ~%"),
-      printf(fh, "  flux_surf_nodal[~a] = alpha_quad*Jfavg_quad - fabs(alpha_quad)*Jfjump_quad; ~%", j0index)
+      if scheme = "upwind" then (
+        printf(fh, "  flux_surf_nodal[~a] = alpha_quad*Jfavg_quad - fabs(alpha_quad)*Jfjump_quad; ~%", j0index)
+      ) else if scheme = "central" then (
+        printf(fh, "  flux_surf_nodal[~a] = alpha_quad*Jfavg_quad; ~%", j0index)
+      )       else if scheme = "downwind" then (
+        printf(fh, "  flux_surf_nodal[~a] = alpha_quad*Jfavg_quad + fabs(alpha_quad)*Jfjump_quad; ~%", j0index)
+      ) else (
+        /* Stop the script if an invalid scheme is provided. */
+        error("Invalid flux scheme provided. Options are: upwind, central, downwind.")
+      )
     ),
     printf(fh, "~%")
   ),
 
-  /* Do the quad nodal to modal ops directly here*/ 
-  /*printf(fh, "~%"),
-  printf(fh, "  double *fmodal = &flux_surf[~a]; ~%", length(bSurf)*(surfDir-1)),
-  flux_surf_nodal_e : doExpand1(flux_surf_nodal,basisNodal),
-  fmodproj_e : fullratsimp(calcInnerProdList(surfIntVars, 1, bSurf, flux_surf_nodal_e)),
-
-  for i : 1 thru length(fmodproj_e) do (
-    printf(fh, "  fmodal[~a] = ~a; ~%", i-1, float(expand(fmodproj_e[i])))
-  ),
-
-  printf(fh, "~%"),*/
-
   /*Calculate the cfl*/
   pOrderCFL : polyOrder,
   printf(fh, "~%"),
@@ -296,4 +340,4 @@ buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
   flush_output(fh),
   printf(fh, "} ~%")
 
-)$
+)$
\ No newline at end of file
diff --git a/maxima/g0/gk_collisionless/gk_collisionless_flux-surf-vpar.mac b/maxima/g0/gk_collisionless/gk_collisionless_flux-surf-vpar.mac
index a172abb2..409eab38 100644
--- a/maxima/g0/gk_collisionless/gk_collisionless_flux-surf-vpar.mac
+++ b/maxima/g0/gk_collisionless/gk_collisionless_flux-surf-vpar.mac
@@ -5,13 +5,15 @@ load("scifac")$
 load("utilities_gyrokinetic")$
 fpprec : 24$
 
-buildGKFluxVparESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no_by, edge) := block(
+buildGKFluxVparKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no_by, edge, scheme) := block(
   [pDim,varsC,bC,varsP,bP,vSub,numC,numP,surfVar,varLabel,dirLabel,surfIntVars,surf_cvars,surf_vvars,
    surfIntVarsC,bSurfC,surfNodes,nodeVars,bSurf,basisNodal,configNodes,numSurfNodes,numConfigNodes,
    numVelNodes,tempVars,tempBasis,NSurfIndexing,numNodesIndexing,d,rdx2vec,rdv2vec,rdSurfVar2,
    bmagBasis,phi_e,bmag_e,vmap_e,vmapSq_e,vmap_prime_e,evPoint,hamil_e,hamil_c,replaceList,
    hamilCvar,hamilNoZero_c,JfL_e,JfR_e,JfL_c,JfR_c,JfL_nodes,JfR_nodes,vmap_prime_nodes,vpardim,
-   dH_dz_nodes,i,j,j0index,j1index,pOrderCFL,vprimeStr
+   dH_dz_nodes,i,i0index,i1index,j,j0index,j1index,pOrderCFL,vprimeStr,apar_e,apar_nodes,dA_dx_nodes,k,
+   apardot_e,apardot_nodes,mvpar_j1,dH_dx_e,dH_dy_e,dH_dz_e,dA_dx_e,dA_dy_e,dA_dz_e,
+   gradHxgradA_nodes, apargradH_nodes
   ],
 
   kill(varsC,varsP,bC,bP),
@@ -57,18 +59,17 @@ buildGKFluxVparESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
     NSurfIndexing    : length(tempBasis),
     numNodesIndexing : length(tempBasis)
   ) else (
-    NSurfIndexing    : NSurf,
-    numNodesIndexing : numNodes
+    error("Only p=1 is currently supported for the surfvpar kernels.")
   ),
 
   print("Working on ", funcNm),
   printf(fh, "GKYL_CU_DH double ~a(
-    const double *w, const double *dxv, 
-    const double *vmap_prime_l, const double *vmap_prime_r,
-    const double *vmap, const double *vmapSq, const double q_, const double m_,
-    const struct gkyl_dg_vol_geom *dgv, const struct gkyl_gk_dg_vol_geom *gkdgv,
-    const double *bmag, const double *phi, const double *JfL, const double *JfR,
-    double* GKYL_RESTRICT flux_surf) ~%{ ~%", funcNm),
+  const double *w, const double *dxv, 
+  const double *vmap_prime_l, const double *vmap_prime_r,
+  const double *vmap, const double *vmapSq, const double q_, const double m_,
+  const struct gkyl_dg_vol_geom *dgv, const struct gkyl_gk_dg_vol_geom *gkdgv,
+  const double *bmag, const double *phi, const double *apar, const double *apardot, const double *JfL, const double *JfR,
+  double* GKYL_RESTRICT flux_surf) ~%{ ~%", funcNm),
   printf(fh, "  // w[NDIM]: cell-center.~%"),
   printf(fh, "  // dxv[NDIM]: cell length.~%"),
   printf(fh, "  // vmap_prime_l,vmap_prime_r: velocity space mapping derivative in left and right cells.~%"),
@@ -79,6 +80,8 @@ buildGKFluxVparESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
   printf(fh, "  // gkdgv: gyrokinetic volume DG geometry.~%"),
   printf(fh, "  // bmag: magnetic field amplitude.~%"),
   printf(fh, "  // phi: electrostatic potential.~%"),
+  printf(fh, "  // apar: parallel component of vector potential.~%"),
+  printf(fh, "  // apardot: time derivative of parallel component of vector potential.~%"),
   printf(fh, "  // JfL: distribution times total jacobian in left cell.~%"),
   printf(fh, "  // JfR: distribution times total jacobian in right cell.~%"),
   printf(fh, "  // flux_surf: output surface phase space flux in each direction (cdim + 1 components).~%"),
@@ -127,9 +130,18 @@ buildGKFluxVparESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
   flush_output(fh),
   hamilNoZero_c : makelistNoZeros1(hamil_c, hamilCvar),
   /* Expand projected Hamiltonian on basis. */
-  hamil_e : hamilNoZero_c . bP,
+  hamil_e : doExpand(hamilNoZero_c,bP),
   /*hamil_e : subst(surfVar=evPoint,hamil_e),*/
 
+  /* Expand Apar */
+  apar_e : doExpand1(apar, bC),
+  apar_nodes : float(evAtNodes(apar_e,surfNodes,surfIntVars)),
+  /* Compute gradient of Aparallel */
+  dA_dx_nodes : makelist(0, i, 1, cdim),  
+  for i : 1 thru cdim do (
+    dA_dx_nodes[i] : float(evAtNodes(diff(apar_e*rdx2vec[i],varsC[i]),surfNodes,surfIntVars))
+  ),
+
   /*fl and fr */
   JfL_e : doExpand1(JfL, bP),
   JfR_e : doExpand1(JfR, bP),
@@ -157,76 +169,137 @@ buildGKFluxVparESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
     )
   ),
 
-  /* Now calculate apha at all quadrature nodes */
+  if cdim = 3 then (
+    dH_dx_e : diff(hamil_e*rdx2vec[1],varsP[1]),
+    dH_dy_e : diff(hamil_e*rdx2vec[2],varsP[2]),
+    dH_dz_e : diff(hamil_e*rdx2vec[3],varsP[3]),
+    dA_dx_e : diff(apar_e*rdx2vec[1],varsP[1]),
+    dA_dy_e : diff(apar_e*rdx2vec[2],varsP[2]),
+    dA_dz_e : diff(apar_e*rdx2vec[3],varsP[3])
+  ) else if cdim = 2 then (
+    dH_dx_e : diff(hamil_e*rdx2vec[1],varsP[1]),
+    dH_dy_e : 0,
+    dH_dz_e : diff(hamil_e*rdx2vec[2],varsP[2]),
+    dA_dx_e : diff(apar_e*rdx2vec[1],varsP[1]),
+    dA_dy_e : 0,
+    dA_dz_e : diff(apar_e*rdx2vec[2],varsP[2])
+  ) else if cdim = 1 then (
+    dH_dx_e : 0,
+    dH_dy_e : 0,
+    dH_dz_e : diff(hamil_e*rdx2vec[1],varsP[1]),
+    dA_dx_e : 0,
+    dA_dy_e : 0,
+    dA_dz_e : diff(apar_e*rdx2vec[1],varsP[1])
+  ),
+  gradHxgradA_nodes : [
+    float(evAtNodes(dH_dy_e*dA_dz_e - dH_dz_e*dA_dy_e,surfNodes,surfIntVars)),
+    float(evAtNodes(dH_dz_e*dA_dx_e - dH_dx_e*dA_dz_e,surfNodes,surfIntVars)),
+    float(evAtNodes(dH_dx_e*dA_dy_e - dH_dy_e*dA_dx_e,surfNodes,surfIntVars))
+  ],
+
+  apargradH_nodes : makelist(0, i, 1, cdim),
+  for i : 1 thru cdim do (
+     apargradH_nodes[i] : float(evAtNodes(apar_e*diff(hamil_e*rdx2vec[i],varsP[i]),surfNodes,surfIntVars))
+  ),
+  
+  /* Now calculate alpha at all quadrature nodes */
   /*printf(fh, "  double flux_surf_nodal[~a]= {0.0}; ~%", numSurfNodes),*/
   printf(fh, "  double *flux_surf_nodal = &flux_surf[~a]; ~%", NSurfIndexing*(surfDir-1)),
   printf(fh, "  double cfl = 0.0; ~%"),
   printf(fh, "  double bmag_quad = 0.0; ~%"),
   printf(fh, "  double B3_quad = 0.0; ~%"),
-  printf(fh, "  double Jc_quad = 0.0; ~%"),
   printf(fh, "  double dualcurlbhat_quad[3] = {0.0}; ~%"),
-
   printf(fh, "  double alpha_quad = 0.0; ~%"),
   printf(fh, "  double JfL_quad = 0.0; ~%"),
   printf(fh, "  double JfR_quad = 0.0; ~%"),
   printf(fh, "  double Jfavg_quad = 0.0; ~%"),
   printf(fh, "  double Jfjump_quad = 0.0; ~%"),
+  printf(fh, "  double bioverJB_quad[3] = {0.0}; ~%"),
   printf(fh, "~%"),
 
   for i : 1 thru numConfigNodes do (
-    printf(fh, "  bmag_quad = gkdgv[~a].bmag; ~%", i-1), 
-    printf(fh, "  B3_quad = gkdgv[~a].B3; ~%", i-1), 
-    printf(fh, "  Jc_quad = dgv[~a].Jc; ~%", i-1), 
-    printf(fh, "  dualcurlbhat_quad[0] = gkdgv[~a].dualcurlbhat.x[0]; ~%", i-1), 
-    printf(fh, "  dualcurlbhat_quad[1] = gkdgv[~a].dualcurlbhat.x[1]; ~%", i-1), 
-    printf(fh, "  dualcurlbhat_quad[2] = gkdgv[~a].dualcurlbhat.x[2]; ~%", i-1), 
-    printf(fh, "~%"),
+    i0index : i-1,
+    i1index : i,
+    printf(fh, "  bmag_quad = gkdgv[~a].bmag; ~%", i0index), 
+    printf(fh, "  B3_quad = gkdgv[~a].B3; ~%", i0index), 
+    printf(fh, "  dualcurlbhat_quad[0] = gkdgv[~a].dualcurlbhat.x[0]; ~%", i0index), 
+    printf(fh, "  dualcurlbhat_quad[1] = gkdgv[~a].dualcurlbhat.x[1]; ~%", i0index), 
+    printf(fh, "  dualcurlbhat_quad[2] = gkdgv[~a].dualcurlbhat.x[2]; ~%", i0index), 
+
+    printf(fh, "  bioverJB_quad[0] = gkdgv[~a].bioverJB.x[0]; ~%", i0index), 
+    printf(fh, "  bioverJB_quad[1] = gkdgv[~a].bioverJB.x[1]; ~%", i0index), 
+    printf(fh, "  bioverJB_quad[2] = gkdgv[~a].bioverJB.x[2]; ~%", i0index),
+
     for j : 1 thru numVelNodes do (
       j0index : j-1+(i-1)*numVelNodes,
       j1index : j+(i-1)*numVelNodes,
+      mvpar_j1 :  dH_dz_nodes[vpardim][j1index]/vmap_prime_nodes[j1index],
       printf(fh, "~%"),
-      if no_by = true then (
-        printf(fh, "  alpha_quad = -(~a)/m_/bmag_quad * B3_quad ;~%", dH_dz_nodes[cdim][j1index])
-      ),
+
+      /* Start ES term */
+      printf(fh, "  alpha_quad = -(~a)/m_/bmag_quad * B3_quad; ~%", dH_dz_nodes[cdim][j1index]),
+
       if no_by = false then (
-        printf(fh, "  alpha_quad = -(~a)/m_/bmag_quad * B3_quad ", dH_dz_nodes[cdim][j1index]),
-         if cdim = 3 then (
-           for k : 1 thru cdim do (
-             printf(fh, "-(~a)/m_/bmag_quad * 1/q_*dualcurlbhat_quad[~a]*(~a)", dH_dz_nodes[k][j1index], k-1, dH_dz_nodes[vpardim][j1index]/vmap_prime_nodes[j1index])
-           )
-         ),
-         if cdim = 2 then (
-             printf(fh, "-(~a)/m_/bmag_quad * 1/q_*dualcurlbhat_quad[~a]*(~a)", dH_dz_nodes[1][j1index], 0, dH_dz_nodes[vpardim][j1index]/vmap_prime_nodes[j1index]),
-             printf(fh, "-(~a)/m_/bmag_quad * 1/q_*dualcurlbhat_quad[~a]*(~a)", dH_dz_nodes[2][j1index], 2, dH_dz_nodes[vpardim][j1index]/vmap_prime_nodes[j1index])
-         ),
-         if cdim = 1 then (
-             printf(fh, "-(~a)/m_/bmag_quad * 1/q_*dualcurlbhat_quad[~a]*(~a)", dH_dz_nodes[1][j1index], 2, dH_dz_nodes[vpardim][j1index]/vmap_prime_nodes[j1index])
-         ),
-        printf(fh, ";~%")
+        if cdim = 3 then (
+          /* Finish ES term */
+          for k : 1 thru cdim do (
+            printf(fh, "  alpha_quad += -(~a)/m_/bmag_quad * 1/q_*dualcurlbhat_quad[~a]*(~a); ~%", dH_dz_nodes[k][j1index], k-1, mvpar_j1)
+          ),
+          /* EM term using ∇ x (A b) . ∇H  = A . ∇ x b . ∇H + ∇A x b . ∇H = A . ∇ x b . ∇H + ∇H x ∇A . b */
+          for k : 1 thru cdim do (
+            /* - A ∇ x b . ∇H dψ/dvpar = - A (∇ x b)_i dH/dxi dψ/dvpar */
+            printf(fh, "  alpha_quad += -1.0/m_/bmag_quad * (~a)*dualcurlbhat_quad[~a]; ~%", apargradH_nodes[k][j1index], k-1),
+            /* - ∇H x ∇A . b dψ/dvpar  = - b_i/B 1/J eps_ijk dH/dx_j dA/dx_k */
+            printf(fh, "  alpha_quad += -bioverJB_quad[~a]/m_ *(~a); ~%", k-1, gradHxgradA_nodes[k][j1index])
+          )
+        ),
+        if cdim = 2 then (
+          /* Finish ES term */
+          printf(fh, "  alpha_quad += -(~a)/m_/bmag_quad * 1/q_*dualcurlbhat_quad[~a]*(~a); ~%", dH_dz_nodes[1][j1index], 0, mvpar_j1),
+          printf(fh, "  alpha_quad += -(~a)/m_/bmag_quad * 1/q_*dualcurlbhat_quad[~a]*(~a); ~%", dH_dz_nodes[2][j1index], 2, mvpar_j1),
+          /* EM - A * ∇ x b . ∇H dψ/dvpar*/
+          printf(fh, "  alpha_quad += -1.0/m_/bmag_quad * (~a)*dualcurlbhat_quad[~a]; ~%", apargradH_nodes[1][j1index], 0),
+          printf(fh, "  alpha_quad += -1.0/m_/bmag_quad * (~a)*dualcurlbhat_quad[~a]; ~%", apargradH_nodes[2][j1index], 2),
+          /* EM - ∇H x ∇A . b dψ/dvpar */
+          /* Only b_2 term is non zero in 2x2v */
+          printf(fh, "  alpha_quad += -bioverJB_quad[~a]/m_ *(~a); ~%", 2, gradHxgradA_nodes[2][j1index])
+        ),
+        if cdim = 1 then (
+          /* Finish ES term */
+          printf(fh, "  alpha_quad += -(~a)/m_/bmag_quad * 1/q_*dualcurlbhat_quad[~a]*(~a); ~%", dH_dz_nodes[1][j1index], 2, mvpar_j1),
+          /* EM - A * ∇ x b . ∇H dψ/dvpar*/
+          printf(fh, "  alpha_quad += -1.0/m_/bmag_quad * (~a)*dualcurlbhat_quad[~a]; ~%", apargradH_nodes[1][j1index], 2)
+          /* EM - ∇H x ∇A . b dψ/dvpar */
+          /* none */
+        )
       ),
 
+      /* Add EM (apardot) contribution */
+      apardot_e : doExpand1(apardot, bC),
+      apardot_nodes : float(evAtNodes(apardot_e,configNodes,surf_cvars)),
+      printf(fh, "  alpha_quad += -q_/m_*(~a); ~%", apardot_nodes[i1index]),
+
       printf(fh, "~%"),
       printf(fh, "  cfl = fmax(fabs(alpha_quad), fabs(cfl)) ;~%", j0index),
       printf(fh, "  JfL_quad =  (~a)/~a;~%", JfL_nodes[j1index], vmap_prime_l[surfDir-cdim-1]),
       printf(fh, "  JfR_quad =  (~a)/~a;~%", JfR_nodes[j1index], vmap_prime_r[surfDir-cdim-1]),
       printf(fh, "  Jfavg_quad = (JfL_quad + JfR_quad)/2.0 ;~%"),
       printf(fh, "  Jfjump_quad = (JfR_quad - JfL_quad)/2.0 ;~%"),
-      printf(fh, "  flux_surf_nodal[~a] = alpha_quad*Jfavg_quad - fabs(alpha_quad)*Jfjump_quad ;~%", j0index)
+      if scheme = "upwind" then (
+        printf(fh, "  flux_surf_nodal[~a] = alpha_quad*Jfavg_quad - fabs(alpha_quad)*Jfjump_quad ;~%", j0index)
+      ) else if scheme = "central" then (
+        printf(fh, "  flux_surf_nodal[~a] = alpha_quad*Jfavg_quad; ~%", j0index)
+      ) else if scheme = "downwind" then (
+        printf(fh, "  flux_surf_nodal[~a] = alpha_quad*Jfavg_quad + fabs(alpha_quad)*Jfjump_quad; ~%", j0index)
+      ) else (
+        /* Stop the script if an invalid scheme is provided. */
+        error("Invalid flux scheme provided. Options are: upwind, central, downwind.")
+      )
     ),
     printf(fh, "~%")
   ),
 
-  /* Do the quad nodal to modal ops directly here*/ 
-  /*printf(fh, "~%"),
-  printf(fh, "  double *fmodal = &flux_surf[~a]; ~%", NSurfIndexing*(surfDir-1)),
-  flux_surf_nodal_e : doExpand1(flux_surf_nodal,basisNodal),
-  fmodproj_e : fullratsimp(calcInnerProdList(surfIntVars, 1, bSurf, flux_surf_nodal_e)),
-
-  for i : 1 thru length(fmodproj_e) do (
-    printf(fh, "  fmodal[~a] = ~a; ~%", i-1, float(expand(fmodproj_e[i])))
-  ),
-
-  printf(fh, "~%"),*/
+  printf(fh, "~%"),
   /*Calculate the cfl*/
   pOrderCFL : polyOrder,
   if polyOrder=1 then ( pOrderCFL : 2 ),
@@ -239,4 +312,4 @@ buildGKFluxVparESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
   flush_output(fh),
   printf(fh, "} ~%")
 
-)$
+)$
\ No newline at end of file
diff --git a/maxima/g0/gk_collisionless/ms-dg_gyrokinetic-header.mac b/maxima/g0/gk_collisionless/ms-dg_gyrokinetic-header.mac
index 08e41563..06c9a6de 100644
--- a/maxima/g0/gk_collisionless/ms-dg_gyrokinetic-header.mac
+++ b/maxima/g0/gk_collisionless/ms-dg_gyrokinetic-header.mac
@@ -2,6 +2,9 @@
 
 /* ...... USER INPUTS........ */
 
+/* Output directory for generated files. */
+outputDir : "~/max-out/"$
+
 /* Serendipity basis. */
 maxPolyOrder_Ser : 2$
 minCdim_Ser : 1$
@@ -86,10 +89,31 @@ printPrototypes() := block([],
         )
       )   
     )
+  ),
+  /* EM terms */
+  for bInd : 1 thru length(bName) do (
+    for c : minCdim[bInd] thru maxCdim[bInd] do (
+      for gkV : 1 thru length(gkVdims[c]) do (
+        v : gkVdims[c][gkV],
+  
+        maxPolyOrderB : maxPolyOrder[bInd],
+        maxPolyOrderB : 1, /* Only declare p=1 kernels for 3x2v */
+        for polyOrder : 1 thru maxPolyOrderB do (
+          printf(fh, "GKYL_CU_DH double dg_gyrokinetic_add_apar_vol_~ax~av_~a_p~a(const double *w, const double *dxv, 
+            const double *vmap, const double *vmapSq, const double q_, const double m_, 
+            const double *bmag, const double *jacobtot_inv, const double *dualcurlbhatoverB, const double *bioverJB,
+            const double *b_i, const double *phi, const double *apar,
+            const double *fin, double* GKYL_RESTRICT out); ~%", c, v, bName[bInd], polyOrder),
+          printf(fh, "GKYL_CU_DH double dg_gyrokinetic_add_apardot_vol_~ax~av_~a_p~a(const double *vmap, const double q_, const double m_,
+            const double *apardot, const double *fin, double* GKYL_RESTRICT out); ~%", c, v, bName[bInd], polyOrder),
+          printf(fh, "~%")
+        )
+      )   
+    )
   )
 )$
 
-fh : openw("~/max-out/gkyl_dg_gyrokinetic_kernels.h")$
+fh : openw(sconcat(outputDir, "gkyl_dg_gyrokinetic_kernels.h"))$
 printf(fh, "#pragma once~%")$
 printf(fh, "~%")$
 printf(fh, "#include <math.h>~%")$
diff --git a/maxima/g0/gk_collisionless/ms-dg_gyrokinetic-surf.mac b/maxima/g0/gk_collisionless/ms-dg_gyrokinetic-surf.mac
index 18f5db41..f8dc8b49 100644
--- a/maxima/g0/gk_collisionless/ms-dg_gyrokinetic-surf.mac
+++ b/maxima/g0/gk_collisionless/ms-dg_gyrokinetic-surf.mac
@@ -9,6 +9,9 @@ load("gk_collisionless/dg_gk-surf")$
 
 /* ...... USER INPUTS........ */
 
+/* Output directory for generated files. */
+outputDir : "~/max-out/"$
+
 /* Serendipity basis. */
 minPolyOrder_Ser : 1$
 maxPolyOrder_Ser : 1$
@@ -60,7 +63,7 @@ for bInd : 1 thru length(bName) do (
       for polyOrder : minPolyOrder[bInd] thru maxPolyOrderB do (
         for dir : 1 thru c do (
           /* Advection in configuration space.*/
-          fname : sconcat("~/max-out/dg_gyrokinetic_surf",clabels[dir],"_", c, "x", v, "v_", bName[bInd], "_p",polyOrder, ".c"),
+          fname : sconcat(outputDir,"dg_gyrokinetic_surf",clabels[dir],"_", c, "x", v, "v_", bName[bInd], "_p",polyOrder, ".c"),
           disp(printf(false,"Creating surface file: ~a",fname)),
   
           fh : openw(fname),
@@ -70,7 +73,7 @@ for bInd : 1 thru length(bName) do (
           close(fh),
 
           /* Advection in configuration space in the skin cell (for boundary flux operations) .*/
-          fname : sconcat("~/max-out/dg_gyrokinetic_boundary_surf",clabels[dir],"_", c, "x", v, "v_", bName[bInd], "_p",polyOrder, ".c"),
+          fname : sconcat(outputDir,"dg_gyrokinetic_boundary_surf",clabels[dir],"_", c, "x", v, "v_", bName[bInd], "_p",polyOrder, ".c"),
           disp(printf(false,"Creating surface file: ~a",fname)),
 
           fh : openw(fname),
@@ -81,7 +84,7 @@ for bInd : 1 thru length(bName) do (
         ),
 
         /* Advection in velocity space.*/
-        fname : sconcat("~/max-out/dg_gyrokinetic_surf",vlabels[1],"_", c, "x", v, "v_", bName[bInd], "_p",polyOrder, ".c"),
+        fname : sconcat(outputDir,"dg_gyrokinetic_surf",vlabels[1],"_", c, "x", v, "v_", bName[bInd], "_p",polyOrder, ".c"),
         disp(printf(false,"Creating surface file: ~a",fname)),
   
         fh : openw(fname),
@@ -91,7 +94,7 @@ for bInd : 1 thru length(bName) do (
         close(fh),
 
         /* Advection in velocity space in the skin cell along vpar (for zero-flux BCs).*/
-        fname : sconcat("~/max-out/dg_gyrokinetic_boundary_surf",vlabels[1],"_", c, "x", v, "v_", bName[bInd], "_p",polyOrder, ".c"),
+        fname : sconcat(outputDir,"dg_gyrokinetic_boundary_surf",vlabels[1],"_", c, "x", v, "v_", bName[bInd], "_p",polyOrder, ".c"),
         disp(printf(false,"Creating surface file: ~a",fname)),
   
         fh : openw(fname),
diff --git a/maxima/g0/gk_collisionless/ms-dg_gyrokinetic-vol.mac b/maxima/g0/gk_collisionless/ms-dg_gyrokinetic-vol.mac
index ef0b8a40..965586ee 100644
--- a/maxima/g0/gk_collisionless/ms-dg_gyrokinetic-vol.mac
+++ b/maxima/g0/gk_collisionless/ms-dg_gyrokinetic-vol.mac
@@ -7,6 +7,9 @@ load("gk_collisionless/dg_gk-vol")$
 
 /* ...... USER INPUTS........ */
 
+/* Output directory for generated files. */
+outputDir : "~/max-out/"$
+
 /* Serendipity basis. */
 minPolyOrder_Ser : 1$
 maxPolyOrder_Ser : 1$
@@ -44,19 +47,35 @@ for bInd : 1 thru length(bName) do (
       maxPolyOrderB : maxPolyOrder[bInd],
       if (c=3) then maxPolyOrderB : 1, /* Only generate p=1 kernels for 3x2v */
       for polyOrder : minPolyOrder[bInd] thru maxPolyOrderB do (
-        fname : sconcat("~/max-out/dg_gyrokinetic_vol_", c, "x", v, "v_", bName[bInd], "_p", polyOrder, ".c"),
+        fname : sconcat(outputDir,"dg_gyrokinetic_vol_", c, "x", v, "v_", bName[bInd], "_p", polyOrder, ".c"),
         disp(printf(false,"Creating volume file: ~a",fname)),
   
         fh : openw(fname),
         printf(fh, "#include <gkyl_dg_gyrokinetic_kernels.h> ~%"),
-
         funcName : sconcat("dg_gyrokinetic_vol_", c, "x", v, "v_", bName[bInd], "_p", polyOrder),
         buildGKVolKernel(fh, funcName, c, v, bName[bInd], polyOrder, bVarsList, false),
         close(fh),
 
+        /* Add em kernels */
+        fname : sconcat(outputDir,"dg_gyrokinetic_add_apar_vol_", c, "x", v, "v_", bName[bInd], "_p", polyOrder, ".c"),
+        disp(printf(false,"Creating volume file (add em): ~a",fname)),
+        fh : openw(fname),
+        printf(fh, "#include <gkyl_dg_gyrokinetic_kernels.h> ~%"),
+        funcName : sconcat("dg_gyrokinetic_add_apar_vol_", c, "x", v, "v_", bName[bInd], "_p", polyOrder),
+        addAparGKEMVolKernel(fh, funcName, c, v, bName[bInd], polyOrder, bVarsList, false),
+        close(fh),
+
+        fname : sconcat(outputDir,"dg_gyrokinetic_add_apardot_vol_", c, "x", v, "v_", bName[bInd], "_p", polyOrder, ".c"),
+        disp(printf(false,"Creating volume file (add em): ~a",fname)),
+        fh : openw(fname),
+        printf(fh, "#include <gkyl_dg_gyrokinetic_kernels.h> ~%"),
+        funcName : sconcat("dg_gyrokinetic_add_apardot_vol_", c, "x", v, "v_", bName[bInd], "_p", polyOrder),
+        addApardotGKEMVolKernel(fh, funcName, c, v, bName[bInd], polyOrder, bVarsList, false),
+        close(fh),
+
         /* if cdim > 1, also generate a set of kernels for the case where there is no toroidal field (by = 0) */
         if (c > 1) then (
-          fname : sconcat("~/max-out/dg_gyrokinetic_no_by_vol_", c, "x", v, "v_", bName[bInd], "_p", polyOrder, ".c"),
+          fname : sconcat(outputDir,"dg_gyrokinetic_no_by_vol_", c, "x", v, "v_", bName[bInd], "_p", polyOrder, ".c"),
           disp(printf(false,"Creating volume file (no by): ~a",fname)),
     
           fh : openw(fname),
diff --git a/maxima/g0/gk_collisionless/ms-gk_collisionless_flux-header.mac b/maxima/g0/gk_collisionless/ms-gk_collisionless_flux-header.mac
index ebef045a..9b66d599 100644
--- a/maxima/g0/gk_collisionless/ms-gk_collisionless_flux-header.mac
+++ b/maxima/g0/gk_collisionless/ms-gk_collisionless_flux-header.mac
@@ -2,8 +2,11 @@
 
 /* ...... USER INPUTS........ */
 
+/* Output directory for generated files. */
+outputDir : "~/max-out/"$
+
 /* Serendipity basis. */
-maxPolyOrder_Ser : 2$
+maxPolyOrder_Ser : 1$
 minCdim_Ser : 1$
 minVdim_Ser : 1$
 maxCdim_Ser : 3$
@@ -43,9 +46,16 @@ byStr : ["", "no_by_"]$
 mb_bcOpt : [[false,true],[false,true],[false,true]]$
 mb_bcStr : ["", "multib_boundary_"]$
 
+/* Options for writing kernels at surface and edge surfaces. */
+edgeBool : [false, true]$
+edgeOpt : ["surf", "edge_surf"]$
+
 printPrototypes() := block([],
+
   for bInd : 1 thru length(bName) do (
+
     for c : minCdim[bInd] thru maxCdim[bInd] do (
+
       for gkV : 1 thru length(gkVdims[c]) do (
         v : gkVdims[c][gkV],
   
@@ -63,50 +73,63 @@ printPrototypes() := block([],
 
               for surfDir : 1 thru c do (
                 dirlabel : varsC[surfDir],
-                extraargs : "const struct gkyl_dg_surf_geom *dgs, const struct gkyl_gk_dg_surf_geom *gkdgs, ",
-                vprimeargs : "",
-
-                printf(fh, "GKYL_CU_DH double gk_collisionless_flux_~a~asurf~a_~ax~av_~a_p~a(
-                  const double *w, const double *dxv,
-                  ~a
-                  const double *vmap, const double *vmapSq, const double q_, const double m_, 
-                  ~a
-                  const double *bmag, const double *jacobgeo_rat_surfL, const double *jacobgeo_rat_surfR,
-                  const double *phi, const double *JfL, const double *JfR, 
-                  double* GKYL_RESTRICT flux_surf); ~%", no_byStr, mb_boundStr, dirlabel, c, v, bName[bInd], polyOrder, vprimeargs, extraargs),  
-
-                printf(fh, "GKYL_CU_DH double gk_collisionless_flux_~a~aedge_surf~a_~ax~av_~a_p~a(
-                  const double *w, const double *dxv,
-                  ~a
-                  const double *vmap, const double *vmapSq, const double q_, const double m_,
-                  ~a
-                  const double *bmag, const double *jacobgeo_rat_surfL, const double *jacobgeo_rat_surfR,
-                  const double *phi, const double *JfL, const double *JfR, 
-                  double* GKYL_RESTRICT flux_surf); ~%", no_byStr, mb_boundStr, dirlabel, c, v, bName[bInd], polyOrder, vprimeargs, extraargs)
+
+                for edgeI : 1 thru 2 do (
+                  edge : edgeBool[edgeI],
+                  edgeStr : edgeOpt[edgeI],
+                
+                  funcName : sconcat("gk_collisionless_flux_",no_byStr,mb_boundStr,edgeStr,dirlabel,"_",c,"x",v,"v_",bName[bInd],"_p",polyOrder),
+                  printf(fh, "GKYL_CU_DH double ~a(
+    const double *w, const double *dxv,
+    const double *vmap, const double *vmapSq, const double q_, const double m_,
+    const struct gkyl_dg_surf_geom *dgs, const struct gkyl_gk_dg_surf_geom *gkdgs,
+    const double *bmag, const double *jacobgeo_rat_surfL, const double *jacobgeo_rat_surfR,
+    const double *phi, const double *apar, const double *apardot, const double *JfL, const double *JfR,
+    double* GKYL_RESTRICT flux_surf); ~%", funcName)
+                )
               )
             ),
 
             dirlabel : varsV[1],
-            extraargs : "const struct gkyl_dg_vol_geom *dgv, const struct gkyl_gk_dg_vol_geom *gkdgv, ",
-            vprimeargs : "const double *vmap_prime_l, const double *vmap_prime_r, ",
-
-            printf(fh, "GKYL_CU_DH double gk_collisionless_flux_~asurf~a_~ax~av_~a_p~a(
-              const double *w, const double *dxv,
-              ~a
-              const double *vmap, const double *vmapSq, const double q_, const double m_, 
-              ~a
-              const double *bmag, const double *phi, const double *JfL, const double *JfR, 
-              double* GKYL_RESTRICT flux_surf); ~%", no_byStr, dirlabel, c, v, bName[bInd], polyOrder, vprimeargs, extraargs)
-          ),
-
-          printf(fh, "~%")
+            edgeStr : edgeOpt[1],
+            edge : false,
+            funcName : sconcat("gk_collisionless_flux_",no_byStr,edgeStr,dirlabel,"_",c,"x",v,"v_",bName[bInd],"_p",polyOrder),
+            printf(fh, "GKYL_CU_DH double ~a(
+  const double *w, const double *dxv,
+  const double *vmap_prime_l, const double *vmap_prime_r,
+  const double *vmap, const double *vmapSq, const double q_, const double m_,
+  const struct gkyl_dg_vol_geom *dgv, const struct gkyl_gk_dg_vol_geom *gkdgv,
+  const double *bmag, const double *phi, const double *apar, const double *apardot, const double *JfL, const double *JfR,
+  double* GKYL_RESTRICT flux_surf); ~%", funcName),
+
+            printf(fh, "~%")
+          )
         )
       )   
     )
-  )
+  ),
+
+  printf(fh,"GKYL_CU_DH double gk_collisionless_flux_surfconf_none(
+  const double *w, const double *dxv,
+  const double *vmap, const double *vmapSq, const double q_, const double m_,
+  const struct gkyl_dg_surf_geom *dgs, const struct gkyl_gk_dg_surf_geom *gkdgs,
+  const double *bmag, const double *jacobgeo_rat_surfL, const double *jacobgeo_rat_surfR,
+  const double *phi, const double *apar, const double *apardot,
+  const double *JfL, const double *JfR, double* GKYL_RESTRICT flux_surf);~%"),
+
+  printf(fh, "~%"),
+
+  printf(fh,"GKYL_CU_DH double gk_collisionless_flux_surfvpar_none(
+  const double *w, const double *dxv,
+  const double *vmap_prime_l, const double *vmap_prime_r,
+  const double *vmap, const double *vmapSq, const double q_, const double m_, 
+  const struct gkyl_dg_vol_geom *dgv, const struct gkyl_gk_dg_vol_geom *gkdgv, 
+  const double *bmag, const double *phi, const double *apar, const double *apardot, 
+  const double *JfL, const double *JfR, double* GKYL_RESTRICT flux_surf);~%")
+
 )$
 
-fh : openw("~/max-out/gkyl_gk_collisionless_flux_kernels.h")$
+fh : openw(sconcat(outputDir,"gkyl_gk_collisionless_flux_kernels.h"))$
 printf(fh, "#pragma once~%")$
 printf(fh, "~%")$
 printf(fh, "#include <math.h>~%")$
diff --git a/maxima/g0/gk_collisionless/ms-gk_collisionless_flux.mac b/maxima/g0/gk_collisionless/ms-gk_collisionless_flux.mac
index 25581b85..5b208422 100644
--- a/maxima/g0/gk_collisionless/ms-gk_collisionless_flux.mac
+++ b/maxima/g0/gk_collisionless/ms-gk_collisionless_flux.mac
@@ -8,6 +8,9 @@ load("gk_collisionless/gk_collisionless_flux-surf-vpar")$
 
 /* ...... USER INPUTS........ */
 
+/* Output directory for generated files. */
+outputDir : "~/max-out/"$
+
 /* Serendipity basis. */
 minPolyOrder_Ser : 1$
 maxPolyOrder_Ser : 1$
@@ -45,61 +48,102 @@ byStr : ["", "no_by_"]$
 mb_bcOpt : [[false,true],[false,true],[false,true]]$
 mb_bcStr : ["", "multib_boundary_"]$
 
+/* Options for writing kernels at surface and edge surfaces. */
+edgeBool : [false, true]$
+edgeOpt : ["surf", "edge_surf"]$
+
 /* Generate kernels of selected types. */
 for bInd : 1 thru length(bName) do (
+  bStr : bName[bInd],
+  /* Loop over configuration space dimensions. */
   for c : minCdim[bInd] thru maxCdim[bInd] do (
+    /* Loop over velocity space dimensions. */
     for gkV : 1 thru length(gkVdims[c]) do (
       v : gkVdims[c][gkV],
-
       maxPolyOrderB : maxPolyOrder[bInd],
       if (c=3) then maxPolyOrderB : 1, /* Only generate p=1 kernels for 3x2v */
-
+      /* Loop over polynomial order. */
       for polyOrder : minPolyOrder[bInd] thru maxPolyOrderB do (
-        for byI : 1 thru length(byOpt[c]) do (
-          no_by : byOpt[c][byI],
-          no_byStr : byStr[byI],
 
+      /* With/without toroidal field loop. */
+      for byI : 1 thru length(byOpt[c]) do (
+        no_by : byOpt[c][byI],
+        no_byStr : byStr[byI],
+
+          /* Singleblock/multiblock loop. */
           for mbI : 1 thru length(mb_bcOpt[c]) do (
             mb_bound : mb_bcOpt[c][mbI],
             mb_boundStr : mb_bcStr[mbI],
 
-            /* Surface flux in direction dir in configuration space.*/
+            /* Configuration space surface direction loop. */
             for dir : 1 thru c do (
-
-              fname : sconcat("~/max-out/gk_collisionless_flux_",no_byStr,mb_boundStr,"surf",clabels[dir],"_", c, "x", v, "v_", bName[bInd], "_p", polyOrder, ".c"),
-              disp(printf(false,"Creating flux surf~a ~a ~a file: ~a",clabels[dir],no_byStr,mb_boundStr,fname)),
-      
-              fh : openw(fname),
-              printf(fh, "#include <gkyl_gk_collisionless_flux_kernels.h> ~%"),
-  
-              funcName : sconcat("gk_collisionless_flux_",no_byStr,mb_boundStr,"surf",clabels[dir],"_", c, "x", v, "v_", bName[bInd], "_p", polyOrder),
-              buildGKFluxConfESKernel(dir, fh, funcName, c, v, bName[bInd], polyOrder, no_by, false, mb_bound),
-              close(fh),
-  
-              fname : sconcat("~/max-out/gk_collisionless_flux_",no_byStr,mb_boundStr,"edge_surf",clabels[dir],"_", c, "x", v, "v_", bName[bInd], "_p", polyOrder, ".c"),
-              disp(printf(false,"Creating flux edge surf~a ~a ~a file: ~a",clabels[dir],no_byStr,mb_boundStr,fname)),
-      
-              fh : openw(fname),
-              printf(fh, "#include <gkyl_gk_collisionless_flux_kernels.h> ~%"),
-  
-              funcName : sconcat("gk_collisionless_flux_",no_byStr,mb_boundStr,"edge_surf",clabels[dir],"_", c, "x", v, "v_", bName[bInd], "_p", polyOrder),
-              buildGKFluxConfESKernel(dir, fh, funcName, c, v, bName[bInd], polyOrder, no_by, true, mb_bound),
-              close(fh)
-            )
-          ),
-
-          /* Surface flux in vparallel direction.*/
-          fname : sconcat("~/max-out/gk_collisionless_flux_",no_byStr,"surf",vlabels[1],"_", c, "x", v, "v_", bName[bInd], "_p", polyOrder, ".c"),
-          disp(printf(false,"Creating flux surfvpar ~a file: ~a",no_byStr,fname)),
-  
-          fh : openw(fname),
-          printf(fh, "#include <gkyl_gk_collisionless_flux_kernels.h> ~%"),
-
-          funcName : sconcat("gk_collisionless_flux_",no_byStr,"surf",vlabels[1],"_", c, "x", v, "v_", bName[bInd], "_p", polyOrder),
-          buildGKFluxVparESKernel(c+1, fh, funcName, c, v, bName[bInd], polyOrder, no_by, false),
-          close(fh)
+              dirStr : clabels[dir],
+
+              /* Surface/edge surface loop. */
+              for edgeI : 1 thru length(edgeBool) do (
+                edge : edgeBool[edgeI],
+                edgeStr : edgeOpt[edgeI],
+
+                scheme : "upwind",
+                fname : sconcat(outputDir,"gk_collisionless_flux_",no_byStr,mb_boundStr,edgeStr,dirStr,"_", c, "x", v, "v_", bStr, "_p", polyOrder, ".c"),
+                disp(printf(false,"Creating ~a flux surf~a ~a ~a file: ~a",edgeStr,dirStr,no_byStr,mb_boundStr,fname)),
+                fh : openw(fname),
+                printf(fh, "#include <gkyl_gk_collisionless_flux_kernels.h> ~%"),
+                funcName : sconcat("gk_collisionless_flux_",no_byStr,mb_boundStr,edgeStr,dirStr,"_", c, "x", v, "v_", bStr, "_p", polyOrder),
+                buildGKFluxConfKernel(dir, fh, funcName, c, v, bStr, polyOrder, no_by, edge, mb_bound, scheme),
+                close(fh)
+              )
+            ),
+
+            /* Surface flux in vparallel direction (no edge).*/
+            dirStr : vlabels[1],
+            edge : false,
+            edgeStr : edgeOpt[1],
+
+            scheme : "upwind",
+            fname : sconcat(outputDir,"gk_collisionless_flux_",no_byStr,edgeStr,dirStr,"_", c, "x", v, "v_", bStr, "_p", polyOrder, ".c"),
+            disp(printf(false,"Creating flux surfvpar ~a file: ~a",no_byStr,fname)),
+            fh : openw(fname),
+            printf(fh, "#include <gkyl_gk_collisionless_flux_kernels.h> ~%"),
+            funcName : sconcat("gk_collisionless_flux_",no_byStr,edgeStr,dirStr,"_", c, "x", v, "v_", bStr, "_p", polyOrder),
+            buildGKFluxVparKernel(c+1, fh, funcName, c, v, bStr, polyOrder, no_by, edge, scheme),
+            close(fh)
+          )
         )
       )
     )
   )
 )$
+
+/* Generate the return zero kernel */
+fname : sconcat(outputDir,"gk_collisionless_flux_surfconf_none.c")$
+disp(printf(false,"Creating return zero kernel file: ~a",fname))$
+fh : openw(fname)$
+printf(fh, "#include <gkyl_gk_collisionless_flux_kernels.h> ~%")$
+printf(fh, "GKYL_CU_DH double gk_collisionless_flux_surfconf_none(~%")$
+printf(fh, "    const double *w, const double *dxv,~%")$
+printf(fh, "    const double *vmap, const double *vmapSq, const double q_, const double m_,~%")$
+printf(fh, "    const struct gkyl_dg_surf_geom *dgs, const struct gkyl_gk_dg_surf_geom *gkdgs, ~%")$
+printf(fh, "    const double *bmag, const double *jacobgeo_rat_surfL, const double *jacobgeo_rat_surfR, ~%")$
+printf(fh, "    const double *phi, const double *apar, const double *apardot, ~%")$
+printf(fh, "    const double *JfL, const double *JfR, double* GKYL_RESTRICT flux_surf) ~%")$
+printf(fh, "{ ~%")$
+printf(fh, "  return 0.0; ~%")$
+printf(fh, "}~%")$
+close(fh)$
+
+fname : sconcat(outputDir,"gk_collisionless_flux_surfvpar_none.c")$
+disp(printf(false,"Creating return zero kernel file: ~a",fname))$
+fh : openw(fname)$
+printf(fh, "#include <gkyl_gk_collisionless_flux_kernels.h> ~%")$
+printf(fh, "GKYL_CU_DH double gk_collisionless_flux_surfvpar_none(~%")$
+printf(fh, "    const double *w, const double *dxv,~%")$
+printf(fh, "    const double *vmap_prime_l, const double *vmap_prime_r,~%")$
+printf(fh, "    const double *vmap, const double *vmapSq, const double q_, const double m_, ~%")$
+printf(fh, "    const struct gkyl_dg_vol_geom *dgv, const struct gkyl_gk_dg_vol_geom *gkdgv, ~%")$
+printf(fh, "    const double *bmag, const double *phi, const double *apar, const double *apardot, ~%")$
+printf(fh, "    const double *JfL, const double *JfR, double* GKYL_RESTRICT flux_surf) ~%")$
+printf(fh, "{ ~%")$
+printf(fh, "  return 0.0; ~%")$
+printf(fh, "}~%")$
+close(fh)$