/[MITgcm]/MITgcm_contrib/gmaze_pv/compute_density.m
ViewVC logotype

Diff of /MITgcm_contrib/gmaze_pv/compute_density.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by gmaze, Thu Jun 22 18:13:30 2006 UTC revision 1.3 by gmaze, Tue Jan 30 22:10:10 2007 UTC
# Line 1  Line 1 
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:
13    % INPUT:
14    % ./netcdf-files/<SNAPSHOT>/<netcdf_THETA>.<netcdf_domain>.<netcdf_suff>
15    % ./netcdf-files/<SNAPSHOT>/<netcdf_SALTanom>.<netcdf_domain>.<netcdf_suff>
16    % OUPUT:
17    % ./netcdf-files/<SNAPSHOT>/RHO.<netcdf_domain>.<netcdf_suff>
18    %
19  % 06/21/2006  % 06/21/2006
20  % gmaze@mit.edu  % gmaze@mit.edu
21  %  %
# Line 18  function compute_density(snapshot) Line 28  function compute_density(snapshot)
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    
# Line 63  nz = length(dpt); Line 74  nz = length(dpt);
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);
# Line 130  nc{ncid}.longname      = ncchar(longname Line 150  nc{ncid}.longname      = ncchar(longname
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
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22