MITgcm_contrib/gmaze_pv/compute_alpha.m

%
% [ALPHA] = compute_alpha(SNAPSHOT)
%
% This function computes the thermal expansion coefficient from
% files of potential temperature THETA and salinity anomaly 
% SALTanom.
% SALTanom is by default a salinity anomaly vs 35PSU.
% If not, (is absolute value) set the global variable is_SALTanom to 0
%
% Files name are:
% INPUT:
% ./netcdf-files/<SNAPSHOT>/<netcdf_THETA>.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/<netcdf_SALTanom>.<netcdf_domain>.<netcdf_suff>
% OUTPUT:
% ./netcdf-files/<SNAPSHOT>/ALPHA.<netcdf_domain>.<netcdf_suff>
%
% with: netcdf_* as global variables
%
% Alpha is computed with the subroutine sw_alpha from package SEAWATER
%
% 06/27/06
% gmaze@mit.edu

function varargout = compute_alpha(snapshot)

global sla toshow
global netcdf_suff netcdf_domain
global netcdf_SALTanom netcdf_THETA
pv_checkpath


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

% Load files:
ferfile = strcat(pathname,sla,snapshot,sla,netcdf_THETA,'.',netcdf_domain,'.',ext);
ncT     = netcdf(ferfile,'nowrite');
[Tlon Tlat Tdpt] = coordfromnc(ncT);

ferfile = strcat(pathname,sla,snapshot,sla,netcdf_SALTanom,'.',netcdf_domain,'.',ext);
ncS   = netcdf(ferfile,'nowrite');
[Slon Slat Sdpt] = coordfromnc(ncS); % but normaly is the same grid as T

% Salinity field ref;
global is_SALTanom
if exist('is_SALTanom')
  if is_SALTanom == 1
    bS = 35;
  else
    bS = 0;
  end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% surface PV flux
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define axis:
nx = length(Tlon) ;
ny = length(Tlat) ;
nz = length(Tdpt) ;


% Pre-allocation:
if toshow,disp('Pre-allocate');end
ALPHA = zeros(nz,ny,nx).*NaN;

% Compute alpha:
for iz = 1 : nz
  if toshow,disp(strcat('Compute alpha for level:',num2str(iz),'/',num2str(nz)));end
  TEMP = ncT{4}(iz,:,:);
  SALT = ncS{4}(iz,:,:) + bS;
  PRES = (0.09998*9.81*Tdpt(iz))*ones(ny,nx);
  ALPHA(iz,:,:) = sw_alpha(SALT,TEMP,PRES,'ptmp');
end %for iz


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Record
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if toshow, disp('record'), end

% General informations: 
netfil     = 'ALPHA';
units      = '1/K';
ncid       = 'ALPHA';
longname   = 'Thermal expansion coefficient';
uniquename = 'ALPHA';

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

% Define axis:
nx = length(Tlon) ;
ny = length(Tlat) ;
nz = length(Tdpt) ;

nc('X') = nx;
nc('Y') = ny;
nc('Z') = nz;
 
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'}(:)         = Tlon;
 
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'}(:)         = Tlat;
 
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'}(:)         = Tdpt;
 
% 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}(:,:,:)        = ALPHA;

nc=close(nc);
close(ncS);
close(ncT);

% Output:
output = struct('ALPHA',ALPHA,'dpt',Tdpt,'lat',Tlat,'lon',Tlon);
switch nargout
 case 1
  varargout(1) = {output};
end
1	gmaze	1.1	%
2	gmaze	1.3	% [ALPHA] = compute_alpha(SNAPSHOT)
3	gmaze	1.1	%
4			% This function computes the thermal expansion coefficient from
5			% files of potential temperature THETA and salinity anomaly
6			% SALTanom.
7	gmaze	1.2	% SALTanom is by default a salinity anomaly vs 35PSU.
8			% If not, (is absolute value) set the global variable is_SALTanom to 0
9	gmaze	1.1	%
10			% Files name are:
11			% INPUT:
12			% ./netcdf-files/<SNAPSHOT>/<netcdf_THETA>.<netcdf_domain>.<netcdf_suff>
13			% ./netcdf-files/<SNAPSHOT>/<netcdf_SALTanom>.<netcdf_domain>.<netcdf_suff>
14			% OUTPUT:
15			% ./netcdf-files/<SNAPSHOT>/ALPHA.<netcdf_domain>.<netcdf_suff>
16			%
17			% with: netcdf_* as global variables
18			%
19			% Alpha is computed with the subroutine sw_alpha from package SEAWATER
20			%
21			% 06/27/06
22			% gmaze@mit.edu
23
24			function varargout = compute_alpha(snapshot)
25
26			global sla toshow
27			global netcdf_suff netcdf_domain
28			global netcdf_SALTanom netcdf_THETA
29			pv_checkpath
30
31
32			% Path and extension to find netcdf-files:
33			pathname = strcat('netcdf-files',sla);
34			ext = netcdf_suff;
35
36			% Load files:
37			ferfile = strcat(pathname,sla,snapshot,sla,netcdf_THETA,'.',netcdf_domain,'.',ext);
38			ncT = netcdf(ferfile,'nowrite');
39			[Tlon Tlat Tdpt] = coordfromnc(ncT);
40
41			ferfile = strcat(pathname,sla,snapshot,sla,netcdf_SALTanom,'.',netcdf_domain,'.',ext);
42			ncS = netcdf(ferfile,'nowrite');
43			[Slon Slat Sdpt] = coordfromnc(ncS); % but normaly is the same grid as T
44
45	gmaze	1.2	% Salinity field ref;
46			global is_SALTanom
47			if exist('is_SALTanom')
48			if is_SALTanom == 1
49			bS = 35;
50			else
51			bS = 0;
52			end
53			end
54	gmaze	1.1
55
56			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57			% surface PV flux
58			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
59
60			% Define axis:
61			nx = length(Tlon) ;
62			ny = length(Tlat) ;
63			nz = length(Tdpt) ;
64
65
66			% Pre-allocation:
67			if toshow,disp('Pre-allocate');end
68			ALPHA = zeros(nz,ny,nx).*NaN;
69
70			% Compute alpha:
71			for iz = 1 : nz
72			if toshow,disp(strcat('Compute alpha for level:',num2str(iz),'/',num2str(nz)));end
73			TEMP = ncT{4}(iz,:,:);
74	gmaze	1.2	SALT = ncS{4}(iz,:,:) + bS;
75	gmaze	1.1	PRES = (0.099989.81Tdpt(iz))*ones(ny,nx);
76			ALPHA(iz,:,:) = sw_alpha(SALT,TEMP,PRES,'ptmp');
77			end %for iz
78
79
80			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
81			% Record
82			%%%%%%%%%%%%%%%%%%%%%%%%%%%%
83			if toshow, disp('record'), end
84
85			% General informations:
86			netfil = 'ALPHA';
87			units = '1/K';
88			ncid = 'ALPHA';
89			longname = 'Thermal expansion coefficient';
90			uniquename = 'ALPHA';
91
92			% Open output file:
93			nc = netcdf(strcat(pathname,sla,snapshot,sla,netfil,'.',netcdf_domain,'.',ext),'clobber');
94
95			% Define axis:
96			nx = length(Tlon) ;
97			ny = length(Tlat) ;
98			nz = length(Tdpt) ;
99
100			nc('X') = nx;
101			nc('Y') = ny;
102			nc('Z') = nz;
103
104			nc{'X'} = ncfloat('X');
105			nc{'X'}.uniquename = ncchar('X');
106			nc{'X'}.long_name = ncchar('longitude');
107			nc{'X'}.gridtype = nclong(0);
108			nc{'X'}.units = ncchar('degrees_east');
109			nc{'X'}(:) = Tlon;
110
111			nc{'Y'} = ncfloat('Y');
112			nc{'Y'}.uniquename = ncchar('Y');
113			nc{'Y'}.long_name = ncchar('latitude');
114			nc{'Y'}.gridtype = nclong(0);
115			nc{'Y'}.units = ncchar('degrees_north');
116			nc{'Y'}(:) = Tlat;
117
118			nc{'Z'} = ncfloat('Z');
119			nc{'Z'}.uniquename = ncchar('Z');
120			nc{'Z'}.long_name = ncchar('depth');
121			nc{'Z'}.gridtype = nclong(0);
122			nc{'Z'}.units = ncchar('m');
123			nc{'Z'}(:) = Tdpt;
124
125			% And main field:
126			nc{ncid} = ncfloat('Z', 'Y', 'X');
127			nc{ncid}.units = ncchar(units);
128			nc{ncid}.missing_value = ncfloat(NaN);
129			nc{ncid}.FillValue_ = ncfloat(NaN);
130			nc{ncid}.longname = ncchar(longname);
131			nc{ncid}.uniquename = ncchar(uniquename);
132			nc{ncid}(:,:,:) = ALPHA;
133
134			nc=close(nc);
135	gmaze	1.3	close(ncS);
136			close(ncT);
137	gmaze	1.1
138	gmaze	1.2	% Output:
139	gmaze	1.3	output = struct('ALPHA',ALPHA,'dpt',Tdpt,'lat',Tlat,'lon',Tlon);
140	gmaze	1.2	switch nargout
141			case 1
142	gmaze	1.3	varargout(1) = {output};
143	gmaze	1.2	end