MITgcm_contrib/gmaze_pv/A_compute_potential_density.m

%
% [ST,LON,LAT,DPT] = A_compute_potential_density(SNAPSHOT)
%
% For a time snapshot, this program computes the 
% 3D potential density from potential temperature and salinity.
% THETA and SALTanom are supposed to be defined on the same 
% domain and grid.
% SALTanom is by default a salinity anomaly vs 35.
% If not, (is absolute value) set the global variable is_SALTanom to 0
% 
% Files names are:
% INPUT:
% ./netcdf-files/<SNAPSHOT>/<netcdf_THETA>.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/<netcdf_SALTanom>.<netcdf_domain>.<netcdf_suff>
% OUPUT:
% ./netcdf-files/<SNAPSHOT>/SIGMATHETA.<netcdf_domain>.<netcdf_suff>
%
% 06/07/2006
% gmaze@mit.edu
%

  
function A_compute_potential_density(snapshot)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global sla netcdf_THETA netcdf_SALTanom netcdf_domain netcdf_suff
pv_checkpath


%% THETA and SALTanom files name:
filTHETA = strcat(netcdf_THETA   ,'.',netcdf_domain);
filSALTa = strcat(netcdf_SALTanom,'.',netcdf_domain);

%% Path and extension to find them:
pathname = strcat('netcdf-files',sla,snapshot);
ext      = strcat('.',netcdf_suff);

%% Load netcdf files:
ferfile = strcat(pathname,sla,filTHETA,ext);
ncTHETA = netcdf(ferfile,'nowrite');
THETAvariables = var(ncTHETA);

ferfile = strcat(pathname,sla,filSALTa,ext);
ncSALTa = netcdf(ferfile,'nowrite');
SALTavariables = var(ncSALTa);

global is_SALTanom
if exist('is_SALTanom')
  if is_SALTanom == 1
    bS = 35;
  else
    bS = 0;
  end
end

%% Gridding:
% Don't care about the grid here !
% SALTanom and THETA are normaly defined on the same grid
% So we compute sigma_theta on it.

%% Flags:
global toshow % Turn to 1 to follow the computing process


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Now we compute potential density
%% The routine used is densjmd95.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Axis (usual netcdf files):
if toshow,disp('Dim');end
[lon lat dpt] = coordfromnc(ncTHETA);
nx = length(lon);
ny = length(lat);
nz = length(dpt);

% Pre-allocate:
if toshow,disp('Pre-allocate');end
SIGMATHETA = zeros(nz,ny,nx);

% Then compute potential density SIGMATHETA:
for iz = 1 : nz
  if toshow,disp(strcat('Compute potential density at level:',num2str(iz),'/',num2str(nz)));end
  
  S = SALTavariables{4}(iz,:,:) + bS; % Evantualy move the anom to an absolute field
  T = THETAvariables{4}(iz,:,:);
  SIGMATHETA(iz,:,:) = densjmd95(S,T,zeros(ny,nx)) - 1000;
   
  % Eventualy make a plot of the field:
  if 0 & iz==1
    clf;pcolor(squeeze(SIGMATHETA(iz,:,:)));
    shading interp;caxis([20 30]);colorbar
    drawnow
    %M(iz)=getframe; % To make a video
  end %if1
end %for iz


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Record output:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% General informations: 
netfil     = strcat('SIGMATHETA','.',netcdf_domain,'.',netcdf_suff);
units      = 'kg/m^3-1000';
ncid       = 'ST';
longname   = 'Potential Density';
uniquename = 'potential_density';

% Open output file:
nc = netcdf(strcat(pathname,sla,netfil),'clobber');

% Define axis:
nc('X') = nx;
nc('Y') = ny;
nc('Z') = nz;

nc{'X'} = 'X';
nc{'Y'} = 'Y';
nc{'Z'} = 'Z';

nc{'X'}            = ncfloat('X');
nc{'X'}.uniquename = ncchar('X');
nc{'X'}.long_name  = ncchar('longitude');
nc{'X'}.gridtype   = nclong(0);
nc{'X'}.units      = ncchar('degrees_east');
nc{'X'}(:)         = lon;

nc{'Y'}            = ncfloat('Y'); 
nc{'Y'}.uniquename = ncchar('Y');
nc{'Y'}.long_name  = ncchar('latitude');
nc{'Y'}.gridtype   = nclong(0);
nc{'Y'}.units      = ncchar('degrees_north');
nc{'Y'}(:)         = lat;
 
nc{'Z'}            = ncfloat('Z');
nc{'Z'}.uniquename = ncchar('Z');
nc{'Z'}.long_name  = ncchar('depth');
nc{'Z'}.gridtype   = nclong(0);
nc{'Z'}.units      = ncchar('m');
nc{'Z'}(:)         = dpt;

% And main field:
nc{ncid}               = ncfloat('Z', 'Y', 'X'); 
nc{ncid}.units         = ncchar(units);
nc{ncid}.missing_value = ncfloat(NaN);
nc{ncid}.FillValue_    = ncfloat(NaN);
nc{ncid}.longname      = ncchar(longname);
nc{ncid}.uniquename    = ncchar(uniquename);
nc{ncid}(:,:,:)        = SIGMATHETA;

nc=close(nc);


% Output:
switch nargout
 case 1
  varargout(1) = SIGMATHETA;
 case 2
  varargout(1) = SIGMATHETA;
  varargout(2) = lon;
 case 3
  varargout(1) = SIGMATHETA;
  varargout(2) = lon;
  varargout(3) = lat;
 case 4
  varargout(1) = SIGMATHETA;
  varargout(2) = lon;
  varargout(3) = lat;
  varargout(4) = dpt;
end
1	gmaze	1.1	%
2	gmaze	1.4	% [ST,LON,LAT,DPT] = A_compute_potential_density(SNAPSHOT)
3	gmaze	1.1	%
4			% For a time snapshot, this program computes the
5			% 3D potential density from potential temperature and salinity.
6			% THETA and SALTanom are supposed to be defined on the same
7			% domain and grid.
8	gmaze	1.4	% 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	gmaze	1.3	% Files names are:
12			% INPUT:
13			% ./netcdf-files/<SNAPSHOT>/<netcdf_THETA>.<netcdf_domain>.<netcdf_suff>
14			% ./netcdf-files/<SNAPSHOT>/<netcdf_SALTanom>.<netcdf_domain>.<netcdf_suff>
15			% OUPUT:
16			% ./netcdf-files/<SNAPSHOT>/SIGMATHETA.<netcdf_domain>.<netcdf_suff>
17			%
18	gmaze	1.1	% 06/07/2006
19			% gmaze@mit.edu
20			%
21
22
23			function A_compute_potential_density(snapshot)
24
25
26			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
27			%% Setup
28			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
29			global sla netcdf_THETA netcdf_SALTanom netcdf_domain netcdf_suff
30			pv_checkpath
31
32
33			%% THETA and SALTanom files name:
34			filTHETA = strcat(netcdf_THETA ,'.',netcdf_domain);
35			filSALTa = strcat(netcdf_SALTanom,'.',netcdf_domain);
36
37			%% Path and extension to find them:
38			pathname = strcat('netcdf-files',sla,snapshot);
39			ext = strcat('.',netcdf_suff);
40
41			%% Load netcdf files:
42			ferfile = strcat(pathname,sla,filTHETA,ext);
43			ncTHETA = netcdf(ferfile,'nowrite');
44			THETAvariables = var(ncTHETA);
45
46			ferfile = strcat(pathname,sla,filSALTa,ext);
47			ncSALTa = netcdf(ferfile,'nowrite');
48			SALTavariables = var(ncSALTa);
49
50	gmaze	1.4	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	gmaze	1.1	%% Gridding:
60			% Don't care about the grid here !
61			% SALTanom and THETA are normaly defined on the same grid
62			% So we compute sigma_theta on it.
63
64			%% Flags:
65			global toshow % Turn to 1 to follow the computing process
66
67
68			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
69			%% Now we compute potential density
70			%% The routine used is densjmd95.m
71			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
72
73			% Axis (usual netcdf files):
74			if toshow,disp('Dim');end
75			[lon lat dpt] = coordfromnc(ncTHETA);
76			nx = length(lon);
77			ny = length(lat);
78			nz = length(dpt);
79
80			% Pre-allocate:
81			if toshow,disp('Pre-allocate');end
82			SIGMATHETA = zeros(nz,ny,nx);
83
84			% Then compute potential density SIGMATHETA:
85			for iz = 1 : nz
86			if toshow,disp(strcat('Compute potential density at level:',num2str(iz),'/',num2str(nz)));end
87
88	gmaze	1.4	S = SALTavariables{4}(iz,:,:) + bS; % Evantualy move the anom to an absolute field
89	gmaze	1.1	T = THETAvariables{4}(iz,:,:);
90			SIGMATHETA(iz,:,:) = densjmd95(S,T,zeros(ny,nx)) - 1000;
91	gmaze	1.2
92	gmaze	1.1	% Eventualy make a plot of the field:
93	gmaze	1.4	if 0 & iz==1
94	gmaze	1.1	clf;pcolor(squeeze(SIGMATHETA(iz,:,:)));
95	gmaze	1.4	shading interp;caxis([20 30]);colorbar
96	gmaze	1.1	drawnow
97	gmaze	1.4	%M(iz)=getframe; % To make a video
98	gmaze	1.1	end %if1
99			end %for iz
100
101
102
103
104			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
105			%% Record output:
106			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
107
108			% General informations:
109			netfil = strcat('SIGMATHETA','.',netcdf_domain,'.',netcdf_suff);
110			units = 'kg/m^3-1000';
111			ncid = 'ST';
112			longname = 'Potential Density';
113			uniquename = 'potential_density';
114
115			% Open output file:
116			nc = netcdf(strcat(pathname,sla,netfil),'clobber');
117
118			% Define axis:
119			nc('X') = nx;
120			nc('Y') = ny;
121			nc('Z') = nz;
122
123			nc{'X'} = 'X';
124			nc{'Y'} = 'Y';
125			nc{'Z'} = 'Z';
126
127			nc{'X'} = ncfloat('X');
128			nc{'X'}.uniquename = ncchar('X');
129			nc{'X'}.long_name = ncchar('longitude');
130			nc{'X'}.gridtype = nclong(0);
131			nc{'X'}.units = ncchar('degrees_east');
132			nc{'X'}(:) = lon;
133
134			nc{'Y'} = ncfloat('Y');
135			nc{'Y'}.uniquename = ncchar('Y');
136			nc{'Y'}.long_name = ncchar('latitude');
137			nc{'Y'}.gridtype = nclong(0);
138			nc{'Y'}.units = ncchar('degrees_north');
139			nc{'Y'}(:) = lat;
140
141			nc{'Z'} = ncfloat('Z');
142			nc{'Z'}.uniquename = ncchar('Z');
143			nc{'Z'}.long_name = ncchar('depth');
144			nc{'Z'}.gridtype = nclong(0);
145			nc{'Z'}.units = ncchar('m');
146			nc{'Z'}(:) = dpt;
147
148			% And main field:
149			nc{ncid} = ncfloat('Z', 'Y', 'X');
150			nc{ncid}.units = ncchar(units);
151			nc{ncid}.missing_value = ncfloat(NaN);
152			nc{ncid}.FillValue_ = ncfloat(NaN);
153			nc{ncid}.longname = ncchar(longname);
154			nc{ncid}.uniquename = ncchar(uniquename);
155			nc{ncid}(:,:,:) = SIGMATHETA;
156
157			nc=close(nc);
158
159	gmaze	1.4
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