diff --git a/src/physics/cam/gw_drag.F90 b/src/physics/cam/gw_drag.F90
index 9a3a651060..f9e08e5a86 100644
--- a/src/physics/cam/gw_drag.F90
+++ b/src/physics/cam/gw_drag.F90
@@ -1549,7 +1549,7 @@ subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat)
   flx_heat = 0._r8
 
   if ( use_gw_nlgw ) then
-    call gw_nlgw_dp_ml(state1,ptend)
+    call gw_nlgw_dp_ml(state1,ptend,lchnk)
   end if
 
   if (use_gw_convect_dp) then
diff --git a/src/physics/cam/gw_nlgw.F90 b/src/physics/cam/gw_nlgw.F90
index 647b771e68..d9eb02e437 100644
--- a/src/physics/cam/gw_nlgw.F90
+++ b/src/physics/cam/gw_nlgw.F90
@@ -12,6 +12,7 @@ module gw_nlgw
 use cam_logfile,    only: iulog
 use physconst,      only: cappa, pi
 use interpolate_data, only: lininterp
+use cam_history,    only: outfld, addfld
 
 use ftorch
 
@@ -104,10 +105,11 @@ module gw_nlgw
 
 !==========================================================================
 
-subroutine gw_nlgw_dp_ml(state_in, ptend)
+subroutine gw_nlgw_dp_ml(state_in, ptend, lchnk)
 
   ! inputs
   type(physics_state), intent(in) :: state_in
+  integer,             intent(in)    :: lchnk
   ! outputs
   type(physics_ptend), intent(inout) :: ptend
 
@@ -173,6 +175,13 @@ subroutine gw_nlgw_dp_ml(state_in, ptend)
   call flux_to_forcing(uflux, utgw)
   call flux_to_forcing(vflux, vtgw)
 
+  ! Write UTGW and VTGW to file
+  call outfld('UTGW_NL', utgw, ncol, lchnk)
+  call outfld('VTGW_NL', vtgw, ncol, lchnk)
+
+  call outfld('UFLUX_NL', uflux, ncol, lchnk)
+  call outfld('VFLUX_NL', vflux, ncol, lchnk)
+
   ! update the tendencies
   ptend%u(:ncol,:pver) = ptend%u(:ncol,:pver) + utgw(:ncol,:pver)
   ptend%v(:ncol,:pver) = ptend%v(:ncol,:pver) + vtgw(:ncol,:pver)
@@ -219,6 +228,11 @@ subroutine gw_nlgw_dp_init(model_path)
      write(iulog,*)'nlgw model loaded from: ', model_path
   endif
 
+  call addfld('UTGW_NL', (/ 'lev' /), 'A', 'm/s2', 'Nonlinear GW zonal wind tendency')
+  call addfld('VTGW_NL', (/ 'lev' /), 'A', 'm/s2', 'Nonlinear GW meridional wind tendency')
+  call addfld('UFLUX_NL', (/ 'lev' /), 'A', 'm/s', 'Nonlinear GW zonal wind flux')
+  call addfld('VFLUX_NL', (/ 'lev' /), 'A', 'm/s', 'Nonlinear GW meridional wind flux')
+
 end subroutine gw_nlgw_dp_init