--- MITgcm_contrib/gmaze_pv/compute_density.m	2006/07/10 15:09:01	1.2
+++ MITgcm_contrib/gmaze_pv/compute_density.m	2007/02/01 17:02:02	1.4
@@ -1,10 +1,13 @@
 %
-% [] = compute_density(SNAPSHOT)
+% [RHO] = compute_density(SNAPSHOT)
 %
 % For a time snapshot, this program computes the 
-% 3D density from potential temperature and salinity.
+% 3D density from potential temperature and salinity fields.
 % THETA and SALTanom are supposed to be defined on the same 
-% domain and grid.
+% domain and grid. 
+% SALTanom is by default a salinity anomaly vs 35PSU.
+% If not, (is absolute value) set the global variable is_SALTanom to 0
+% 
 %
 % Files names are:
 % INPUT:
@@ -18,13 +21,14 @@
 %
 
   
-function compute_density(snapshot)
+function varargout = compute_density(snapshot)
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Setup
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global sla netcdf_THETA netcdf_SALTanom netcdf_domain netcdf_suff
+global is_SALTanom
 pv_checkpath
 
 
@@ -70,11 +74,20 @@
 if toshow,disp('Pre-allocate');end
 RHO = zeros(nz,ny,nx);
 
+global is_SALTanom
+if exist('is_SALTanom')
+  if is_SALTanom == 1
+    bS = 35;
+  else
+    bS = 0;
+  end
+end
+
 % Then compute density RHO:
 for iz = 1 : nz
   if toshow,disp(strcat('Compute density at level:',num2str(iz),'/',num2str(nz)));end
   
-  S = SALTavariables{4}(iz,:,:) + 35; % Move the anom to an absolute field
+  S = SALTavariables{4}(iz,:,:) + bS; % Move the anom to an absolute field
   T = THETAvariables{4}(iz,:,:);
   P = (0.09998*9.81*dpt(iz))*ones(ny,nx);
   RHO(iz,:,:) = densjmd95(S,T,P);
@@ -137,5 +150,17 @@
 nc{ncid}.uniquename    = ncchar(uniquename);
 nc{ncid}(:,:,:)        = RHO;
 
-nc=close(nc);
 
+
+% Close files:
+close(ncTHETA);
+close(ncSALTa);
+close(nc);
+
+
+% Output:
+output = struct('RHO',RHO,'dpt',dpt,'lat',lat,'lon',lon);
+switch nargout
+ case 1
+  varargout(1) = {output};
+end