1 |
% |
2 |
% [JFz] = compute_JFz(SNAPSHOT) |
3 |
% |
4 |
% Here we compute the PV flux due to frictionnal forces as |
5 |
% JFz = ( TAUx * dSIGMATHETA/dy - TAUy * dSIGMATHETA/dx ) / RHO / EKL |
6 |
% |
7 |
% where: |
8 |
% TAU is the surface wind-stress (N/m2) |
9 |
% SIGMATHETA is the potential density (kg/m3) |
10 |
% RHO is the density (kg/m3) |
11 |
% EKL is the Ekman layer depth (m, positive) |
12 |
% |
13 |
% Files names are: |
14 |
% INPUT: |
15 |
% ./netcdf-files/<SNAPSHOT>/<netcdf_SIGMATHETA>.<netcdf_domain>.<netcdf_suff> |
16 |
% ./netcdf-files/<SNAPSHOT>/<netcdf_TAUX>.<netcdf_domain>.<netcdf_suff> |
17 |
% ./netcdf-files/<SNAPSHOT>/<netcdf_TAUY>.<netcdf_domain>.<netcdf_suff> |
18 |
% ./netcdf-files/<SNAPSHOT>/<netcdf_RHO>.<netcdf_domain>.<netcdf_suff> |
19 |
% ./netcdf-files/<SNAPSHOT>/<netcdf_EKL>.<netcdf_domain>.<netcdf_suff> |
20 |
% OUTPUT: |
21 |
% ./netcdf-files/<SNAPSHOT>/JFz.<netcdf_domain>.<netcdf_suff> |
22 |
% |
23 |
% with netcdf_* as global variables |
24 |
% |
25 |
% 06/27/06 |
26 |
% gmaze@mit.edu |
27 |
|
28 |
function varargout = compute_JFz(snapshot) |
29 |
|
30 |
global sla toshow |
31 |
global netcdf_suff netcdf_domain |
32 |
global netcdf_TAUX netcdf_TAUY netcdf_SIGMATHETA netcdf_EKL netcdf_RHO |
33 |
pv_checkpath |
34 |
|
35 |
|
36 |
% NETCDF file name: |
37 |
filST = netcdf_SIGMATHETA; |
38 |
filTx = netcdf_TAUX; |
39 |
filTy = netcdf_TAUY; |
40 |
filRHO = netcdf_RHO; |
41 |
filH = netcdf_EKL; |
42 |
|
43 |
% Path and extension to find them: |
44 |
pathname = strcat('netcdf-files',sla); |
45 |
ext = netcdf_suff; |
46 |
|
47 |
% Load files: |
48 |
ferfile = strcat(pathname,sla,snapshot,sla,filST,'.',netcdf_domain,'.',ext); |
49 |
ncST = netcdf(ferfile,'nowrite'); |
50 |
[STlon STlat STdpt] = coordfromnc(ncST); |
51 |
|
52 |
ferfile = strcat(pathname,sla,snapshot,sla,filTx,'.',netcdf_domain,'.',ext); |
53 |
ncTx = netcdf(ferfile,'nowrite'); |
54 |
ferfile = strcat(pathname,sla,snapshot,sla,filTy,'.',netcdf_domain,'.',ext); |
55 |
ncTy = netcdf(ferfile,'nowrite'); |
56 |
|
57 |
ferfile = strcat(pathname,sla,snapshot,sla,filRHO,'.',netcdf_domain,'.',ext); |
58 |
ncRHO = netcdf(ferfile,'nowrite'); |
59 |
RHO = ncRHO{4}(1,:,:); |
60 |
|
61 |
ferfile = strcat(pathname,sla,snapshot,sla,filH,'.',netcdf_domain,'.',ext); |
62 |
ncH = netcdf(ferfile,'nowrite'); |
63 |
EKL = ncH{4}(1,:,:); |
64 |
|
65 |
|
66 |
|
67 |
|
68 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
69 |
% First term |
70 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
71 |
|
72 |
% Dim: |
73 |
if toshow, disp('dim'), end |
74 |
nx = length(STlon) ; |
75 |
ny = length(STlat) - 1 ; |
76 |
nz = length(STdpt); |
77 |
|
78 |
% Pre-allocate: |
79 |
if toshow, disp('pre-allocate'), end |
80 |
dSIGMATHETAdy = zeros(nz,ny-1,nx).*NaN; |
81 |
dy = zeros(1,ny).*NaN; |
82 |
STup = zeros(nz,ny); |
83 |
STdw = zeros(nz,ny); |
84 |
|
85 |
% Meridional gradient of SIGMATHETA: |
86 |
if toshow, disp('grad'), end |
87 |
% Assuming the grid is regular, dy is independent of x: |
88 |
[dy b] = meshgrid( m_lldist([1 1]*STlon(1),STlat(1:ny+1) ), STdpt ) ; clear b |
89 |
for ix = 1 : nx |
90 |
if toshow, disp(strcat(num2str(ix),'/',num2str(nx))), end |
91 |
STup = squeeze(ncST{4}(:,2:ny+1,ix)); |
92 |
STdw = squeeze(ncST{4}(:,1:ny,ix)); |
93 |
dSTdy = ( STup - STdw ) ./ dy; |
94 |
% Change horizontal grid point definition to fit with SIGMATHETA: |
95 |
dSTdy = ( dSTdy(:,1:ny-1) + dSTdy(:,2:ny) )./2; |
96 |
dSIGMATHETAdy(:,:,ix) = dSTdy; |
97 |
end %for iy |
98 |
|
99 |
% Make TAUx having same limits: |
100 |
TAUx = ncTx{4}(1,2:ny,:); |
101 |
|
102 |
% Compute first term: TAUx * dSIGMATHETA/dy |
103 |
iz = 1; |
104 |
JFz_a = TAUx .* squeeze(dSIGMATHETAdy(iz,:,:)) ; |
105 |
|
106 |
|
107 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
108 |
% Second term |
109 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
110 |
|
111 |
% Dim: |
112 |
if toshow, disp('dim'), end |
113 |
nx = length(STlon) - 1; |
114 |
ny = length(STlat) ; |
115 |
nz = length(STdpt) ; |
116 |
|
117 |
% Pre-allocate: |
118 |
if toshow, disp('pre-allocate'), end |
119 |
dSIGMATHETAdx = zeros(nz,ny,nx-1).*NaN; |
120 |
dx = zeros(1,nx).*NaN; |
121 |
STup = zeros(nz,nx); |
122 |
STdw = zeros(nz,nx); |
123 |
|
124 |
% Zonal gradient of SIGMATHETA |
125 |
if toshow, disp('grad'), end |
126 |
for iy = 1 : ny |
127 |
if toshow, disp(strcat(num2str(iy),'/',num2str(ny))), end |
128 |
[dx b] = meshgrid( m_lldist(STlon(1:nx+1),[1 1]*STlat(iy)), STdpt ) ; clear b |
129 |
STup = squeeze(ncST{4}(:,iy,2:nx+1)); |
130 |
STdw = squeeze(ncST{4}(:,iy,1:nx)); |
131 |
dSTdx = ( STup - STdw ) ./ dx; |
132 |
% Change horizontal grid point definition to fit with SIGMATHETA: |
133 |
dSTdx = ( dSTdx(:,1:nx-1) + dSTdx(:,2:nx) )./2; |
134 |
dSIGMATHETAdx(:,iy,:) = dSTdx; |
135 |
end %for iy |
136 |
|
137 |
% Make TAUy having same limits: |
138 |
TAUy = ncTy{4}(1,:,2:nx); |
139 |
|
140 |
% Compute second term: TAUy * dSIGMATHETA/dx |
141 |
iz = 1; |
142 |
JFz_b = TAUy .* squeeze(dSIGMATHETAdx(iz,:,:)) ; |
143 |
|
144 |
|
145 |
|
146 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
147 |
% Finish ... |
148 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
149 |
% Then make all terms having same limits: |
150 |
nx = length(STlon) ; |
151 |
ny = length(STlat) ; |
152 |
nz = length(STdpt) ; |
153 |
JFz_a = squeeze(JFz_a(:,2:nx-1)); |
154 |
JFz_b = squeeze(JFz_b(2:ny-1,:)); |
155 |
delta_e = squeeze(EKL(2:ny-1,2:nx-1)); |
156 |
rho = squeeze(RHO(2:ny-1,2:nx-1)); |
157 |
|
158 |
% and finish: |
159 |
JFz = (JFz_a - JFz_b)./delta_e./rho; |
160 |
|
161 |
|
162 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
163 |
% Record |
164 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
165 |
if toshow, disp('record'), end |
166 |
|
167 |
% General informations: |
168 |
netfil = 'JFz'; |
169 |
units = 'kg/m3/s2'; |
170 |
ncid = 'JFz'; |
171 |
longname = 'Vertical PV flux due to frictional forces'; |
172 |
uniquename = 'JFz'; |
173 |
|
174 |
% Open output file: |
175 |
nc = netcdf(strcat(pathname,sla,snapshot,sla,netfil,'.',netcdf_domain,'.',ext),'clobber'); |
176 |
|
177 |
% Define axis: |
178 |
nx = length(STlon) ; |
179 |
ny = length(STlat) ; |
180 |
nz = 1 ; |
181 |
|
182 |
nc('X') = nx-2; |
183 |
nc('Y') = ny-2; |
184 |
nc('Z') = nz; |
185 |
|
186 |
nc{'X'} = ncfloat('X'); |
187 |
nc{'X'}.uniquename = ncchar('X'); |
188 |
nc{'X'}.long_name = ncchar('longitude'); |
189 |
nc{'X'}.gridtype = nclong(0); |
190 |
nc{'X'}.units = ncchar('degrees_east'); |
191 |
nc{'X'}(:) = STlon(2:nx-1); |
192 |
|
193 |
nc{'Y'} = ncfloat('Y'); |
194 |
nc{'Y'}.uniquename = ncchar('Y'); |
195 |
nc{'Y'}.long_name = ncchar('latitude'); |
196 |
nc{'Y'}.gridtype = nclong(0); |
197 |
nc{'Y'}.units = ncchar('degrees_north'); |
198 |
nc{'Y'}(:) = STlat(2:ny-1); |
199 |
|
200 |
nc{'Z'} = ncfloat('Z'); |
201 |
nc{'Z'}.uniquename = ncchar('Z'); |
202 |
nc{'Z'}.long_name = ncchar('depth'); |
203 |
nc{'Z'}.gridtype = nclong(0); |
204 |
nc{'Z'}.units = ncchar('m'); |
205 |
nc{'Z'}(:) = STdpt(1); |
206 |
|
207 |
% And main field: |
208 |
nc{ncid} = ncfloat('Z', 'Y', 'X'); |
209 |
nc{ncid}.units = ncchar(units); |
210 |
nc{ncid}.missing_value = ncfloat(NaN); |
211 |
nc{ncid}.FillValue_ = ncfloat(NaN); |
212 |
nc{ncid}.longname = ncchar(longname); |
213 |
nc{ncid}.uniquename = ncchar(uniquename); |
214 |
nc{ncid}(:,:,:) = JFz; |
215 |
|
216 |
nc=close(nc); |
217 |
|
218 |
|
219 |
% Output: |
220 |
output = struct('JFz',JFz,'lat',STlat(2:ny-1),'lon',STlon(2:nx-1)); |
221 |
switch nargout |
222 |
case 1 |
223 |
varargout(1) = {output}; |
224 |
end |