MITgcm_contrib/gmaze_pv/compute_alpha.m

%
% [] = compute_alpha(SNAPSHOT)
%
% This function computes the thermal expansion coefficient from
% files of potential temperature THETA and salinity anomaly 
% SALTanom.
%
% 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


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


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