--- MITgcm_contrib/gmaze_pv/compute_density.m 2006/07/10 15:09:01 1.2 +++ MITgcm_contrib/gmaze_pv/compute_density.m 2007/01/30 22:10:10 1.3 @@ -1,10 +1,13 @@ % -% [] = compute_density(SNAPSHOT) +% [RHO,LON,LAT,DPT] = 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: @@ -25,6 +28,7 @@ %% 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,28 @@ nc{ncid}.uniquename = ncchar(uniquename); nc{ncid}(:,:,:) = RHO; -nc=close(nc); + +% Close files: +close(ncTHETA); +close(ncSALTa); +close(nc); + + +% Output: +switch nargout + case 1 + varargout(1) = RHO; + case 2 + varargout(1) = RHO; + varargout(2) = lon; + case 3 + varargout(1) = RHO; + varargout(2) = lon; + varargout(3) = lat; + case 4 + varargout(1) = RHO; + varargout(2) = lon; + varargout(3) = lat; + varargout(4) = dpt; +end