1 |
% |
% |
2 |
% [] = compute_density(SNAPSHOT) |
% [RHO,LON,LAT,DPT] = compute_density(SNAPSHOT) |
3 |
% |
% |
4 |
% For a time snapshot, this program computes the |
% For a time snapshot, this program computes the |
5 |
% 3D density from potential temperature and salinity. |
% 3D density from potential temperature and salinity fields. |
6 |
% THETA and SALTanom are supposed to be defined on the same |
% THETA and SALTanom are supposed to be defined on the same |
7 |
% domain and grid. |
% domain and grid. |
8 |
|
% SALTanom is by default a salinity anomaly vs 35PSU. |
9 |
|
% If not, (is absolute value) set the global variable is_SALTanom to 0 |
10 |
|
% |
11 |
% |
% |
12 |
% Files names are: |
% Files names are: |
13 |
% INPUT: |
% INPUT: |
28 |
%% Setup |
%% Setup |
29 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
30 |
global sla netcdf_THETA netcdf_SALTanom netcdf_domain netcdf_suff |
global sla netcdf_THETA netcdf_SALTanom netcdf_domain netcdf_suff |
31 |
|
global is_SALTanom |
32 |
pv_checkpath |
pv_checkpath |
33 |
|
|
34 |
|
|
74 |
if toshow,disp('Pre-allocate');end |
if toshow,disp('Pre-allocate');end |
75 |
RHO = zeros(nz,ny,nx); |
RHO = zeros(nz,ny,nx); |
76 |
|
|
77 |
|
global is_SALTanom |
78 |
|
if exist('is_SALTanom') |
79 |
|
if is_SALTanom == 1 |
80 |
|
bS = 35; |
81 |
|
else |
82 |
|
bS = 0; |
83 |
|
end |
84 |
|
end |
85 |
|
|
86 |
% Then compute density RHO: |
% Then compute density RHO: |
87 |
for iz = 1 : nz |
for iz = 1 : nz |
88 |
if toshow,disp(strcat('Compute density at level:',num2str(iz),'/',num2str(nz)));end |
if toshow,disp(strcat('Compute density at level:',num2str(iz),'/',num2str(nz)));end |
89 |
|
|
90 |
S = SALTavariables{4}(iz,:,:) + 35; % Move the anom to an absolute field |
S = SALTavariables{4}(iz,:,:) + bS; % Move the anom to an absolute field |
91 |
T = THETAvariables{4}(iz,:,:); |
T = THETAvariables{4}(iz,:,:); |
92 |
P = (0.09998*9.81*dpt(iz))*ones(ny,nx); |
P = (0.09998*9.81*dpt(iz))*ones(ny,nx); |
93 |
RHO(iz,:,:) = densjmd95(S,T,P); |
RHO(iz,:,:) = densjmd95(S,T,P); |
150 |
nc{ncid}.uniquename = ncchar(uniquename); |
nc{ncid}.uniquename = ncchar(uniquename); |
151 |
nc{ncid}(:,:,:) = RHO; |
nc{ncid}(:,:,:) = RHO; |
152 |
|
|
|
nc=close(nc); |
|
153 |
|
|
154 |
|
|
155 |
|
% Close files: |
156 |
|
close(ncTHETA); |
157 |
|
close(ncSALTa); |
158 |
|
close(nc); |
159 |
|
|
160 |
|
|
161 |
|
% Output: |
162 |
|
switch nargout |
163 |
|
case 1 |
164 |
|
varargout(1) = RHO; |
165 |
|
case 2 |
166 |
|
varargout(1) = RHO; |
167 |
|
varargout(2) = lon; |
168 |
|
case 3 |
169 |
|
varargout(1) = RHO; |
170 |
|
varargout(2) = lon; |
171 |
|
varargout(3) = lat; |
172 |
|
case 4 |
173 |
|
varargout(1) = RHO; |
174 |
|
varargout(2) = lon; |
175 |
|
varargout(3) = lat; |
176 |
|
varargout(4) = dpt; |
177 |
|
end |