MITgcm_contrib/gmaze_pv/compute_QEk.m

%
% [QEk] = compute_QEk(SNAPSHOT)
%
% Here we compute the lateral heat flux induced by Ekman currents
% from JFz, the PV flux induced by frictional forces:
% QEk = - Cw * EKL * JFz / alpha / f
% where:
%  Cw = 4187 J/kg/K is the specific heat of seawater
%  EKL is the Ekman layer depth (m)
%  JFz is the PV flux (kg/m3/s2)
%  alpha = 2.5*E-4 1/K is the thermal expansion coefficient
%  f = 2*OMEGA*sin(LAT) is the Coriolis parameter
%
% This allows a direct comparison with the net surface heat flux Qnet
% which forces the surface Pv flux due to diabatic processes.
%   
% Remind that:
% JFz = ( TAUx * dSIGMATHETA/dy - TAUy * dSIGMATHETA/dx ) / RHO / EKL
%
% Files names are:
% INPUT:
% ./netcdf-files/<SNAPSHOT>/<netcdf_JFz>.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/<netcdf_EKL>.<netcdf_domain>.<netcdf_suff>
% OUPUT:
% ./netcdf-files/<SNAPSHOT>/QEk.<netcdf_domain>.<netcdf_suff>
%
% with netcdf_* as global variables
%
% 06/27/06
% gmaze@mit.edu

function varargout = compute_QEk(snapshot)

global sla toshow
global netcdf_suff netcdf_domain
global netcdf_JFz netcdf_EKL
pv_checkpath


% NETCDF file name:
filJFz  = netcdf_JFz;
filEKL  = netcdf_EKL;

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

% Load files:
ferfile = strcat(pathname,sla,snapshot,sla,filJFz,'.',netcdf_domain,'.',ext);
ncJFz   = netcdf(ferfile,'nowrite');
JFz     = ncJFz{4}(1,:,:);
[JFzlon JFzlat JFzdpt] = coordfromnc(ncJFz);

ferfile = strcat(pathname,sla,snapshot,sla,filEKL,'.',netcdf_domain,'.',ext);
ncEKL   = netcdf(ferfile,'nowrite');
EKL     = ncEKL{4}(1,:,:);
[EKLlon EKLlat EKLdpt] = coordfromnc(ncEKL);

% Make them having same limits:
% (JFz is defined with first/last points removed from the EKL grid)
nx = length(JFzlon) ;
ny = length(JFzlat) ;
nz = length(JFzdpt) ;
EKL = squeeze(EKL(2:ny+1,2:nx+1));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Dim:
if toshow, disp('dim'), end
nx = length(JFzlon) ;
ny = length(JFzlat) ;
nz = length(JFzdpt) ;

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

% Planetary vorticity:
f = 2*(2*pi/86400)*sin(JFzlat*pi/180);
[a f]=meshgrid(JFzlon,f); clear a c

% Coefficient:
Cw = 4187;
al = 2.5*10^(-4); % Average surface value of alpha
coef = - Cw / al;

% Compute flux:
QEk = coef.* EKL .* JFz ./ f;


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

% General informations: 
netfil     = 'QEk';
units      = 'W/m2';
ncid       = 'QEk';
longname   = 'Lateral heat flux induced by Ekman currents';
uniquename = 'QEk';

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

% Define axis:
nx = length(JFzlon) ;
ny = length(JFzlat) ;
nz = 1 ;

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'}(:)         = JFzlon;
 
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'}(:)         = JFzlat;
 
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'}(:)         = JFzdpt(1);
 
% 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}(:,:,:)        = QEk;

nc=close(nc);


% Output:
output = struct('QEk',QEk,'lat',JFzlat,'lon',JFzlon);
switch nargout
 case 1
  varargout(1) = {output};
end
1	%
2	% [QEk] = compute_QEk(SNAPSHOT)
3	%
4	% Here we compute the lateral heat flux induced by Ekman currents
5	% from JFz, the PV flux induced by frictional forces:
6	% QEk = - Cw * EKL * JFz / alpha / f
7	% where:
8	% Cw = 4187 J/kg/K is the specific heat of seawater
9	% EKL is the Ekman layer depth (m)
10	% JFz is the PV flux (kg/m3/s2)
11	% alpha = 2.5*E-4 1/K is the thermal expansion coefficient
12	% f = 2OMEGAsin(LAT) is the Coriolis parameter
13	%
14	% This allows a direct comparison with the net surface heat flux Qnet
15	% which forces the surface Pv flux due to diabatic processes.
16	%
17	% Remind that:
18	% JFz = ( TAUx * dSIGMATHETA/dy - TAUy * dSIGMATHETA/dx ) / RHO / EKL
19	%
20	% Files names are:
21	% INPUT:
22	% ./netcdf-files/<SNAPSHOT>/<netcdf_JFz>.<netcdf_domain>.<netcdf_suff>
23	% ./netcdf-files/<SNAPSHOT>/<netcdf_EKL>.<netcdf_domain>.<netcdf_suff>
24	% OUPUT:
25	% ./netcdf-files/<SNAPSHOT>/QEk.<netcdf_domain>.<netcdf_suff>
26	%
27	% with netcdf_* as global variables
28	%
29	% 06/27/06
30	% gmaze@mit.edu
31
32	function varargout = compute_QEk(snapshot)
33
34	global sla toshow
35	global netcdf_suff netcdf_domain
36	global netcdf_JFz netcdf_EKL
37	pv_checkpath
38
39
40	% NETCDF file name:
41	filJFz = netcdf_JFz;
42	filEKL = netcdf_EKL;
43
44	% Path and extension to find them:
45	pathname = strcat('netcdf-files',sla);
46	ext = netcdf_suff;
47
48	% Load files:
49	ferfile = strcat(pathname,sla,snapshot,sla,filJFz,'.',netcdf_domain,'.',ext);
50	ncJFz = netcdf(ferfile,'nowrite');
51	JFz = ncJFz{4}(1,:,:);
52	[JFzlon JFzlat JFzdpt] = coordfromnc(ncJFz);
53
54	ferfile = strcat(pathname,sla,snapshot,sla,filEKL,'.',netcdf_domain,'.',ext);
55	ncEKL = netcdf(ferfile,'nowrite');
56	EKL = ncEKL{4}(1,:,:);
57	[EKLlon EKLlat EKLdpt] = coordfromnc(ncEKL);
58
59	% Make them having same limits:
60	% (JFz is defined with first/last points removed from the EKL grid)
61	nx = length(JFzlon) ;
62	ny = length(JFzlat) ;
63	nz = length(JFzdpt) ;
64	EKL = squeeze(EKL(2:ny+1,2:nx+1));
65
66
67	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
68	%
69	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
70
71	% Dim:
72	if toshow, disp('dim'), end
73	nx = length(JFzlon) ;
74	ny = length(JFzlat) ;
75	nz = length(JFzdpt) ;
76
77	% Pre-allocate:
78	if toshow, disp('pre-allocate'), end
79	QEk = zeros(nz,ny,nx).*NaN;
80
81	% Planetary vorticity:
82	f = 2(2pi/86400)sin(JFzlatpi/180);
83	[a f]=meshgrid(JFzlon,f); clear a c
84
85	% Coefficient:
86	Cw = 4187;
87	al = 2.5*10^(-4); % Average surface value of alpha
88	coef = - Cw / al;
89
90	% Compute flux:
91	QEk = coef.* EKL .* JFz ./ f;
92
93
94
95
96
97	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
98	% Record
99	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
100	if toshow, disp('record'), end
101
102	% General informations:
103	netfil = 'QEk';
104	units = 'W/m2';
105	ncid = 'QEk';
106	longname = 'Lateral heat flux induced by Ekman currents';
107	uniquename = 'QEk';
108
109	% Open output file:
110	nc = netcdf(strcat(pathname,sla,snapshot,sla,netfil,'.',netcdf_domain,'.',ext),'clobber');
111
112	% Define axis:
113	nx = length(JFzlon) ;
114	ny = length(JFzlat) ;
115	nz = 1 ;
116
117	nc('X') = nx;
118	nc('Y') = ny;
119	nc('Z') = nz;
120
121	nc{'X'} = ncfloat('X');
122	nc{'X'}.uniquename = ncchar('X');
123	nc{'X'}.long_name = ncchar('longitude');
124	nc{'X'}.gridtype = nclong(0);
125	nc{'X'}.units = ncchar('degrees_east');
126	nc{'X'}(:) = JFzlon;
127
128	nc{'Y'} = ncfloat('Y');
129	nc{'Y'}.uniquename = ncchar('Y');
130	nc{'Y'}.long_name = ncchar('latitude');
131	nc{'Y'}.gridtype = nclong(0);
132	nc{'Y'}.units = ncchar('degrees_north');
133	nc{'Y'}(:) = JFzlat;
134
135	nc{'Z'} = ncfloat('Z');
136	nc{'Z'}.uniquename = ncchar('Z');
137	nc{'Z'}.long_name = ncchar('depth');
138	nc{'Z'}.gridtype = nclong(0);
139	nc{'Z'}.units = ncchar('m');
140	nc{'Z'}(:) = JFzdpt(1);
141
142	% And main field:
143	nc{ncid} = ncfloat('Z', 'Y', 'X');
144	nc{ncid}.units = ncchar(units);
145	nc{ncid}.missing_value = ncfloat(NaN);
146	nc{ncid}.FillValue_ = ncfloat(NaN);
147	nc{ncid}.longname = ncchar(longname);
148	nc{ncid}.uniquename = ncchar(uniquename);
149	nc{ncid}(:,:,:) = QEk;
150
151	nc=close(nc);
152
153
154
155	% Output:
156	output = struct('QEk',QEk,'lat',JFzlat,'lon',JFzlon);
157	switch nargout
158	case 1
159	varargout(1) = {output};
160	end