diff --git a/docs/conf.py b/docs/conf.py
index 12c00f0a..cf305423 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -63,7 +63,7 @@
 # The short X.Y version.
 version = '1.6'
 # The full version, including alpha/beta/rc tags.
-release = '1.6.0'
+release = '1.6.1'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
diff --git a/docs/releases.rst b/docs/releases.rst
index 17d70c57..c1032df3 100644
--- a/docs/releases.rst
+++ b/docs/releases.rst
@@ -10,6 +10,15 @@ Note, backwards compatibility is defined as the ability to return to a previous
 
 The `github release notes <https://github.com/PrincetonUniversity/GPEC/releases>`_ are reproduced below.
 
+GPEC v1.6.1
+===========
+
+Fixes
+-----
+- GPEC - The Clebsch displacements ``xigradpsi``, ``xigradpsi_dpsi`` and ``xigradalpha`` in the netCDF output are labeled unitless instead of m.
+- GPEC - ``Jbgradpsi_rmatch`` is now filled when ``jac_out`` differs from ``jac_type`` or ``tmag_out=0``; it previously held uninitialized values in those runs.
+
+
 GPEC v1.6.0
 ===========
 
diff --git a/gpec/gpout.f b/gpec/gpout.f
index 897b60bd..5d42c2bb 100644
--- a/gpec/gpout.f
+++ b/gpec/gpout.f
@@ -3418,7 +3418,7 @@ SUBROUTINE gpout_xbnormal(egnum,xspmn,spots,interpspot,nspot)
       REAL(r8) :: ximax,rmax,area
 
       REAL(r8), DIMENSION(0:mthsurf) :: delpsi,jacs,dphi
-      COMPLEX(r8), DIMENSION(0:mthsurf) :: xwp_fun,bwp_fun
+      COMPLEX(r8), DIMENSION(0:mthsurf) :: xwp_fun,bwp_fun,bwpr_fun
 
       COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: xmns,ymns,
      $     xnomns,bnomns,bwpmns,bwpmns_rmatch
@@ -3431,7 +3431,7 @@ SUBROUTINE gpout_xbnormal(egnum,xspmn,spots,interpspot,nspot)
       COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: rss,zss,
      $     xnofuns,bnofuns,intbwpmns
 
-      COMPLEX(r8), DIMENSION(mpert) :: interpbwn
+      COMPLEX(r8), DIMENSION(mpert) :: interpbwn,bwpr_mn
       TYPE(cspline_type) :: fsp_sol
 c-----------------------------------------------------------------------
 c     allocation puts memory in heap, avoiding stack overfill
@@ -3532,6 +3532,13 @@ SUBROUTINE gpout_xbnormal(egnum,xspmn,spots,interpspot,nspot)
          IF ((jac_out /= jac_type).OR.(tout==0)) THEN
             CALL gpeq_bcoordsout(bwpmns(istep,:),bno_mn,
      $                           psifac(istep),ji=1)
+            IF(galsol%gal_flag) THEN
+               CALL iscdftb(mfac,mpert,bwpr_fun,mthsurf,bwp_mn_rmatch)
+               bwpr_fun=bwpr_fun/(jacs*delpsi)
+               CALL iscdftf(mfac,mpert,bwpr_fun,mthsurf,bwpr_mn)
+               CALL gpeq_bcoordsout(bwpmns_rmatch(istep,:),bwpr_mn,
+     $                              psifac(istep),ji=1)
+            ENDIF
             IF (msing>0) THEN
                CALL gpeq_bcoordsout(intbwpmns(istep,:),interpbwn,
      $              psifac(istep),ji=0)
@@ -5860,7 +5867,7 @@ SUBROUTINE gpout_xclebsch(egnum,xspmn)
      $      "Contravariant Clebsch angle displacement") )
          ids = (/xp_id,dp_id,xa_id/)
          DO i=1,3
-            CALL check( nf90_put_att(fncid,ids(i),"units","m") )
+            CALL check( nf90_put_att(fncid,ids(i),"units","unitless") )
             CALL check( nf90_put_att(fncid,ids(i),"jacobian",jac_out) )
          ENDDO
          CALL check( nf90_enddef(fncid) )