--- MITgcm_contrib/gmaze_pv/compute_density.m	2006/06/22 18:13:30	1.1
+++ MITgcm_contrib/gmaze_pv/compute_density.m	2007/09/19 14:45:59	1.5
@@ -1,23 +1,34 @@
 %
-% [] = 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:
+% ./netcdf-files/<SNAPSHOT>/<netcdf_THETA>.<netcdf_domain>.<netcdf_suff>
+% ./netcdf-files/<SNAPSHOT>/<netcdf_SALTanom>.<netcdf_domain>.<netcdf_suff>
+% OUPUT:
+% ./netcdf-files/<SNAPSHOT>/RHO.<netcdf_domain>.<netcdf_suff>
+% 
 % 06/21/2006
 % gmaze@mit.edu
 %
 
   
-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
 
 
@@ -27,6 +38,7 @@
 
 %% Path and extension to find them:
 pathname = strcat('netcdf-files',sla,snapshot);
+%pathname = '.';
 ext      = strcat('.',netcdf_suff);
 
 %% Load netcdf files:
@@ -63,11 +75,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);
@@ -130,5 +151,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