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 = 2*OMEGA*sin(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*(2*pi/86400)*sin(JFzlat*pi/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 |