MITgcm_contrib/gmaze_pv/A_compute_potential_density.m

%
% [] = 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.
%
% 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);

%% 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,:,:) + 35; % 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
    clf;pcolor(squeeze(SIGMATHETA(iz,:,:)));
    shading interp;caxis([10 40]);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);

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