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	gmaze	1.1	%
2	gmaze	1.2	% [QEk] = compute_QEk(SNAPSHOT)
3	gmaze	1.1	%
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	gmaze	1.2	function varargout = compute_QEk(snapshot)
33	gmaze	1.1
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	gmaze	1.2
155			% Output:
156			output = struct('QEk',QEk,'lat',JFzlat,'lon',JFzlon);
157			switch nargout
158			case 1
159			varargout(1) = {output};
160			end