MITgcm_contrib/gmaze_pv/compute_density.m

%
% [RHO,LON,LAT,DPT] = compute_density(SNAPSHOT)
%
% For a time snapshot, this program computes the 
% 3D density from potential temperature and salinity fields.
% THETA and SALTanom are supposed to be defined on the same 
% domain and grid. 
% SALTanom is by default a salinity anomaly vs 35PSU.
% 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>/RHO.<netcdf_domain>.<netcdf_suff>
% 
% 06/21/2006
% gmaze@mit.edu
%

  
function compute_density(snapshot)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global sla netcdf_THETA netcdf_SALTanom netcdf_domain netcdf_suff
global is_SALTanom
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);

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

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Now we compute the 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
RHO = zeros(nz,ny,nx);

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

% Then compute density RHO:
for iz = 1 : nz
  if toshow,disp(strcat('Compute density at level:',num2str(iz),'/',num2str(nz)));end
  
  S = SALTavariables{4}(iz,:,:) + bS; % Move the anom to an absolute field
  T = THETAvariables{4}(iz,:,:);
  P = (0.09998*9.81*dpt(iz))*ones(ny,nx);
  RHO(iz,:,:) = densjmd95(S,T,P);
  
end %for iz


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

% General informations: 
netfil     = strcat('RHO','.',netcdf_domain,'.',netcdf_suff);
units      = 'kg/m^3';
ncid       = 'RHO';
longname   = 'Density';
uniquename = '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}(:,:,:)        = RHO;


% Close files:
close(ncTHETA);
close(ncSALTa);
close(nc);


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