--- MITgcm_contrib/gmaze_pv/compute_JBz.m	2006/07/10 15:09:01	1.1
+++ MITgcm_contrib/gmaze_pv/compute_JBz.m	2007/02/01 17:02:02	1.3
@@ -1,5 +1,5 @@
 %
-% [] = compute_JBz(SNAPSHOT)
+% [JBz] = compute_JBz(SNAPSHOT)
 %
 % Here we compute the PV flux due to diabatic processes as
 % JFz = - alpha * f * Qnet / MLD / Cw
@@ -7,7 +7,7 @@
 %  alpha = 2.5*E-4 1/K is the thermal expansion coefficient
 %  f = 2*OMEGA*sin(LAT) is the Coriolis parameter
 %  Qnet is the net surface heat flux (W/m^2), positive downward
-%  MLD is the mixed layer depth (m)
+%  MLD is the mixed layer depth (m, positive)
 %  Cw = 4187 J/kg/K is the specific heat of seawater
 %
 % Files names are:
@@ -22,7 +22,7 @@
 % 06/27/06
 % gmaze@mit.edu
 
-function compute_JBz(snapshot)
+function varargout = compute_JBz(snapshot)
 
 global sla toshow
 global netcdf_suff netcdf_domain
@@ -123,7 +123,7 @@
 nc{'Z'}.long_name  = ncchar('depth');
 nc{'Z'}.gridtype   = nclong(0);
 nc{'Z'}.units      = ncchar('m');
-nc{'Z'}(:)         = Qdpt;
+nc{'Z'}(:)         = Qdpt(1);
  
 % And main field:
 nc{ncid}               = ncfloat('Z', 'Y', 'X'); 
@@ -135,5 +135,14 @@
 nc{ncid}(:,:,:)        = JBz;
 
 nc=close(nc);
+close(ncQ);
+close(ncMLD);
 
 
+
+% Output:
+output = struct('JBz',JBz,'lat',Qlat,'lon',Qlon);
+switch nargout
+ case 1
+  varargout(1) = {output};
+end