1 |
% |
2 |
% [JFzx] = compute_JFzx(SNAPSHOT) |
3 |
% |
4 |
% Here we compute the PV flux due to the zonal frictionnal force as |
5 |
% JFzx = ( TAUx * dSIGMATHETA/dy ) / RHO / EKL |
6 |
% |
7 |
% where: |
8 |
% TAUx is the surface zonal 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_RHO>.<netcdf_domain>.<netcdf_suff> |
18 |
% ./netcdf-files/<SNAPSHOT>/<netcdf_EKL>.<netcdf_domain>.<netcdf_suff> |
19 |
% OUTPUT: |
20 |
% ./netcdf-files/<SNAPSHOT>/JFzx.<netcdf_domain>.<netcdf_suff> |
21 |
% |
22 |
% with netcdf_* as global variables |
23 |
% |
24 |
% 06/04/12 |
25 |
% gmaze@mit.edu |
26 |
|
27 |
function varargout = compute_JFzx(snapshot) |
28 |
|
29 |
global sla toshow |
30 |
global netcdf_suff netcdf_domain |
31 |
global netcdf_TAUX netcdf_SIGMATHETA netcdf_EKL netcdf_RHO |
32 |
pv_checkpath |
33 |
|
34 |
|
35 |
% NETCDF file name: |
36 |
filST = netcdf_SIGMATHETA; |
37 |
filTx = netcdf_TAUX; |
38 |
filRHO = netcdf_RHO; |
39 |
filH = netcdf_EKL; |
40 |
|
41 |
% Path and extension to find them: |
42 |
pathname = strcat('netcdf-files',sla); |
43 |
ext = netcdf_suff; |
44 |
|
45 |
% Load files: |
46 |
ferfile = strcat(pathname,sla,snapshot,sla,filST,'.',netcdf_domain,'.',ext); |
47 |
ncST = netcdf(ferfile,'nowrite'); |
48 |
[STlon STlat STdpt] = coordfromnc(ncST); |
49 |
|
50 |
ferfile = strcat(pathname,sla,snapshot,sla,filTx,'.',netcdf_domain,'.',ext); |
51 |
ncTx = netcdf(ferfile,'nowrite'); |
52 |
|
53 |
ferfile = strcat(pathname,sla,snapshot,sla,filRHO,'.',netcdf_domain,'.',ext); |
54 |
ncRHO = netcdf(ferfile,'nowrite'); |
55 |
RHO = ncRHO{4}(1,:,:); |
56 |
|
57 |
ferfile = strcat(pathname,sla,snapshot,sla,filH,'.',netcdf_domain,'.',ext); |
58 |
ncH = netcdf(ferfile,'nowrite'); |
59 |
EKL = ncH{4}(1,:,:); |
60 |
|
61 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
62 |
% First term |
63 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
64 |
|
65 |
% Dim: |
66 |
if toshow, disp('dim'), end |
67 |
nx = length(STlon) ; |
68 |
ny = length(STlat) - 1 ; |
69 |
|
70 |
% Pre-allocate: |
71 |
if toshow, disp('pre-allocate'), end |
72 |
dSIGMATHETAdy = zeros(ny-1,nx).*NaN; |
73 |
dy = zeros(1,ny).*NaN; |
74 |
STup = zeros(1,ny); |
75 |
STdw = zeros(1,ny); |
76 |
|
77 |
% Meridional gradient of SIGMATHETA: |
78 |
if toshow, disp('grad'), end |
79 |
% Assuming the grid is regular, dy is independent of x: |
80 |
dy = m_lldist([1 1]*STlon(1),STlat(1:ny+1) ) ; |
81 |
for ix = 1 : nx |
82 |
if toshow, disp(strcat(num2str(ix),'/',num2str(nx))), end |
83 |
STup = squeeze(ncST{4}(1,2:ny+1,ix)); |
84 |
STdw = squeeze(ncST{4}(1,1:ny,ix)); |
85 |
dSTdy = ( STup - STdw ) ./ dy; |
86 |
% Change horizontal grid point definition to fit with SIGMATHETA: |
87 |
dSTdy = ( dSTdy(1:ny-1) + dSTdy(2:ny) )./2; |
88 |
dSIGMATHETAdy(:,ix) = dSTdy; |
89 |
end %for iy |
90 |
|
91 |
% Make TAUx having same limits: |
92 |
TAUx = ncTx{4}(1,2:ny,:); |
93 |
|
94 |
% Compute first term: TAUx * dSIGMATHETA/dy |
95 |
JFz_a = TAUx .* dSIGMATHETAdy ; |
96 |
|
97 |
|
98 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
99 |
% Finish ... |
100 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
101 |
% Then make all terms having same limits: |
102 |
nx = length(STlon) ; |
103 |
ny = length(STlat) ; |
104 |
JFz_a = squeeze(JFz_a(:,2:nx-1)); |
105 |
delta_e = squeeze(EKL(2:ny-1,2:nx-1)); |
106 |
rho = squeeze(RHO(2:ny-1,2:nx-1)); |
107 |
|
108 |
% and finish: |
109 |
JFz = JFz_a./delta_e./rho; |
110 |
|
111 |
|
112 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
113 |
% Record |
114 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
115 |
if toshow, disp('record'), end |
116 |
|
117 |
% General informations: |
118 |
netfil = 'JFzx'; |
119 |
units = 'kg/m3/s2'; |
120 |
ncid = 'JFzx'; |
121 |
longname = 'Vertical PV flux due to the zonal frictional force'; |
122 |
uniquename = 'JFzx'; |
123 |
|
124 |
% Open output file: |
125 |
nc = netcdf(strcat(pathname,sla,snapshot,sla,netfil,'.',netcdf_domain,'.',ext),'clobber'); |
126 |
|
127 |
% Define axis: |
128 |
nx = length(STlon) ; |
129 |
ny = length(STlat) ; |
130 |
nz = 1 ; |
131 |
|
132 |
nc('X') = nx-2; |
133 |
nc('Y') = ny-2; |
134 |
nc('Z') = nz; |
135 |
|
136 |
nc{'X'} = ncfloat('X'); |
137 |
nc{'X'}.uniquename = ncchar('X'); |
138 |
nc{'X'}.long_name = ncchar('longitude'); |
139 |
nc{'X'}.gridtype = nclong(0); |
140 |
nc{'X'}.units = ncchar('degrees_east'); |
141 |
nc{'X'}(:) = STlon(2:nx-1); |
142 |
|
143 |
nc{'Y'} = ncfloat('Y'); |
144 |
nc{'Y'}.uniquename = ncchar('Y'); |
145 |
nc{'Y'}.long_name = ncchar('latitude'); |
146 |
nc{'Y'}.gridtype = nclong(0); |
147 |
nc{'Y'}.units = ncchar('degrees_north'); |
148 |
nc{'Y'}(:) = STlat(2:ny-1); |
149 |
|
150 |
nc{'Z'} = ncfloat('Z'); |
151 |
nc{'Z'}.uniquename = ncchar('Z'); |
152 |
nc{'Z'}.long_name = ncchar('depth'); |
153 |
nc{'Z'}.gridtype = nclong(0); |
154 |
nc{'Z'}.units = ncchar('m'); |
155 |
nc{'Z'}(:) = STdpt(1); |
156 |
|
157 |
% And main field: |
158 |
nc{ncid} = ncfloat('Z', 'Y', 'X'); |
159 |
nc{ncid}.units = ncchar(units); |
160 |
nc{ncid}.missing_value = ncfloat(NaN); |
161 |
nc{ncid}.FillValue_ = ncfloat(NaN); |
162 |
nc{ncid}.longname = ncchar(longname); |
163 |
nc{ncid}.uniquename = ncchar(uniquename); |
164 |
nc{ncid}(:,:,:) = JFz; |
165 |
|
166 |
|
167 |
|
168 |
%%% Close files: |
169 |
close(ncST); |
170 |
close(ncTx); |
171 |
close(ncRHO); |
172 |
close(ncH); |
173 |
close(nc); |
174 |
|
175 |
% Output: |
176 |
output = struct('JFzx',JFz,'lat',STlat(2:ny-1),'lon',STlon(2:nx-1)); |
177 |
switch nargout |
178 |
case 1 |
179 |
varargout(1) = {output}; |
180 |
end |