MITgcm_contrib/gmaze_pv/compute_density.m

%
% [] = COMPUTE_DENSITY(SNAPSHOT)
%
% For a time snapshot, this program computes the 
% 3D density from potential temperature and salinity.
% THETA and SALTanom are supposed to be defined on the same 
% domain and grid.
%
% 06/21/2006
% gmaze@mit.edu
%

  
function compute_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);

%% 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);

% 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,:,:) + 35; % 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;

nc=close(nc);

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