| 1 | % | % | 
| 2 | % [] = compute_density(SNAPSHOT) | % [RHO] = 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: | 
| 21 | % | % | 
| 22 |  |  | 
| 23 |  |  | 
| 24 | function compute_density(snapshot) | function varargout = compute_density(snapshot) | 
| 25 |  |  | 
| 26 |  |  | 
| 27 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | 
| 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 |  |  | 
| 38 |  |  | 
| 39 | %% Path and extension to find them: | %% Path and extension to find them: | 
| 40 | pathname = strcat('netcdf-files',sla,snapshot); | pathname = strcat('netcdf-files',sla,snapshot); | 
| 41 |  | %pathname = '.'; | 
| 42 | ext      = strcat('.',netcdf_suff); | ext      = strcat('.',netcdf_suff); | 
| 43 |  |  | 
| 44 | %% Load netcdf files: | %% Load netcdf files: | 
| 75 | if toshow,disp('Pre-allocate');end | if toshow,disp('Pre-allocate');end | 
| 76 | RHO = zeros(nz,ny,nx); | RHO = zeros(nz,ny,nx); | 
| 77 |  |  | 
| 78 |  | global is_SALTanom | 
| 79 |  | if exist('is_SALTanom') | 
| 80 |  | if is_SALTanom == 1 | 
| 81 |  | bS = 35; | 
| 82 |  | else | 
| 83 |  | bS = 0; | 
| 84 |  | end | 
| 85 |  | end | 
| 86 |  |  | 
| 87 | % Then compute density RHO: | % Then compute density RHO: | 
| 88 | for iz = 1 : nz | for iz = 1 : nz | 
| 89 | 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 | 
| 90 |  |  | 
| 91 | S = SALTavariables{4}(iz,:,:) + 35; % Move the anom to an absolute field | S = SALTavariables{4}(iz,:,:) + bS; % Move the anom to an absolute field | 
| 92 | T = THETAvariables{4}(iz,:,:); | T = THETAvariables{4}(iz,:,:); | 
| 93 | P = (0.09998*9.81*dpt(iz))*ones(ny,nx); | P = (0.09998*9.81*dpt(iz))*ones(ny,nx); | 
| 94 | RHO(iz,:,:) = densjmd95(S,T,P); | RHO(iz,:,:) = densjmd95(S,T,P); | 
| 151 | nc{ncid}.uniquename    = ncchar(uniquename); | nc{ncid}.uniquename    = ncchar(uniquename); | 
| 152 | nc{ncid}(:,:,:)        = RHO; | nc{ncid}(:,:,:)        = RHO; | 
| 153 |  |  | 
|  | nc=close(nc); |  | 
| 154 |  |  | 
| 155 |  |  | 
| 156 |  | % Close files: | 
| 157 |  | close(ncTHETA); | 
| 158 |  | close(ncSALTa); | 
| 159 |  | close(nc); | 
| 160 |  |  | 
| 161 |  |  | 
| 162 |  | % Output: | 
| 163 |  | output = struct('RHO',RHO,'dpt',dpt,'lat',lat,'lon',lon); | 
| 164 |  | switch nargout | 
| 165 |  | case 1 | 
| 166 |  | varargout(1) = {output}; | 
| 167 |  | end |