--- MITgcm_contrib/gmaze_pv/compute_alpha.m	2006/07/10 15:09:01	1.1
+++ MITgcm_contrib/gmaze_pv/compute_alpha.m	2007/01/30 22:10:10	1.2
@@ -1,9 +1,11 @@
 %
-% [] = compute_alpha(SNAPSHOT)
+% [ALPHA,LON,LAT,DPT] = compute_alpha(SNAPSHOT)
 %
 % This function computes the thermal expansion coefficient from
 % files of potential temperature THETA and salinity anomaly 
 % SALTanom.
+% SALTanom is by default a salinity anomaly vs 35PSU.
+% If not, (is absolute value) set the global variable is_SALTanom to 0
 %
 % Files name are:
 % INPUT:
@@ -40,6 +42,15 @@
 ncS   = netcdf(ferfile,'nowrite');
 [Slon Slat Sdpt] = coordfromnc(ncS); % but normaly is the same grid as T
 
+% Salinity field ref;
+global is_SALTanom
+if exist('is_SALTanom')
+  if is_SALTanom == 1
+    bS = 35;
+  else
+    bS = 0;
+  end
+end
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -60,7 +71,7 @@
 for iz = 1 : nz
   if toshow,disp(strcat('Compute alpha for level:',num2str(iz),'/',num2str(nz)));end
   TEMP = ncT{4}(iz,:,:);
-  SALT = ncS{4}(iz,:,:) + 35;
+  SALT = ncS{4}(iz,:,:) + bS;
   PRES = (0.09998*9.81*Tdpt(iz))*ones(ny,nx);
   ALPHA(iz,:,:) = sw_alpha(SALT,TEMP,PRES,'ptmp');
 end %for iz
@@ -123,3 +134,20 @@
 nc=close(nc);
 
 
+% Output:
+switch nargout
+ case 1
+  varargout(1) = ALPHA;
+ case 2
+  varargout(1) = ALPHA;
+  varargout(2) = Tlon;
+ case 3
+  varargout(1) = ALPHA;
+  varargout(2) = Tlon;
+  varargout(3) = Tlat;
+ case 4
+  varargout(1) = ALPHA;
+  varargout(2) = Tlon;
+  varargout(3) = Tlat;
+  varargout(4) = Tdpt;
+end