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

Diff of /MITgcm_contrib/gmaze_pv/A_compute_potential_density.m

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

revision 1.3 by gmaze, Mon Jul 10 15:09:00 2006 UTC revision 1.4 by gmaze, Tue Jan 30 22:10:10 2007 UTC
# Line 1  Line 1 
1  %  %
2  % [] = A_compute_potential_density(SNAPSHOT)  % [ST,LON,LAT,DPT] = A_compute_potential_density(SNAPSHOT)
3  %  %
4  % For a time snapshot, this program computes the  % For a time snapshot, this program computes the
5  % 3D potential density from potential temperature and salinity.  % 3D potential density from potential temperature and salinity.
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 35.
9    % If not, (is absolute value) set the global variable is_SALTanom to 0
10    %
11  % Files names are:  % Files names are:
12  % INPUT:  % INPUT:
13  % ./netcdf-files/<SNAPSHOT>/<netcdf_THETA>.<netcdf_domain>.<netcdf_suff>  % ./netcdf-files/<SNAPSHOT>/<netcdf_THETA>.<netcdf_domain>.<netcdf_suff>
# Line 45  ferfile = strcat(pathname,sla,filSALTa,e Line 47  ferfile = strcat(pathname,sla,filSALTa,e
47  ncSALTa = netcdf(ferfile,'nowrite');  ncSALTa = netcdf(ferfile,'nowrite');
48  SALTavariables = var(ncSALTa);  SALTavariables = var(ncSALTa);
49    
50    global is_SALTanom
51    if exist('is_SALTanom')
52      if is_SALTanom == 1
53        bS = 35;
54      else
55        bS = 0;
56      end
57    end
58    
59  %% Gridding:  %% Gridding:
60  % Don't care about the grid here !  % Don't care about the grid here !
61  % SALTanom and THETA are normaly defined on the same grid  % SALTanom and THETA are normaly defined on the same grid
# Line 74  SIGMATHETA = zeros(nz,ny,nx); Line 85  SIGMATHETA = zeros(nz,ny,nx);
85  for iz = 1 : nz  for iz = 1 : nz
86    if toshow,disp(strcat('Compute potential density at level:',num2str(iz),'/',num2str(nz)));end    if toshow,disp(strcat('Compute potential density at level:',num2str(iz),'/',num2str(nz)));end
87        
88    S = SALTavariables{4}(iz,:,:) + 35; % Move the anom to an absolute field    S = SALTavariables{4}(iz,:,:) + bS; % Evantualy move the anom to an absolute field
89    T = THETAvariables{4}(iz,:,:);    T = THETAvariables{4}(iz,:,:);
90    SIGMATHETA(iz,:,:) = densjmd95(S,T,zeros(ny,nx)) - 1000;    SIGMATHETA(iz,:,:) = densjmd95(S,T,zeros(ny,nx)) - 1000;
91        
92    % Eventualy make a plot of the field:    % Eventualy make a plot of the field:
93    if 0    if 0 & iz==1
94      clf;pcolor(squeeze(SIGMATHETA(iz,:,:)));      clf;pcolor(squeeze(SIGMATHETA(iz,:,:)));
95      shading interp;caxis([10 40]);colorbar      shading interp;caxis([20 30]);colorbar
96      drawnow      drawnow
97      M(iz)=getframe; % To make a video      %M(iz)=getframe; % To make a video
98    end %if1    end %if1
99  end %for iz  end %for iz
100    
# Line 145  nc{ncid}(:,:,:)        = SIGMATHETA; Line 156  nc{ncid}(:,:,:)        = SIGMATHETA;
156    
157  nc=close(nc);  nc=close(nc);
158    
159    
160    % Output:
161    switch nargout
162     case 1
163      varargout(1) = SIGMATHETA;
164     case 2
165      varargout(1) = SIGMATHETA;
166      varargout(2) = lon;
167     case 3
168      varargout(1) = SIGMATHETA;
169      varargout(2) = lon;
170      varargout(3) = lat;
171     case 4
172      varargout(1) = SIGMATHETA;
173      varargout(2) = lon;
174      varargout(3) = lat;
175      varargout(4) = dpt;
176    end
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22