MITgcm_contrib/gmaze_pv/compute_JFzx.m

%
% [JFzx] = compute_JFzx(SNAPSHOT)
%
% Here we compute the PV flux due to the zonal frictionnal force as
% JFzx = ( TAUx * dSIGMATHETA/dy ) / RHO / EKL
%
% where:
%  TAUx is the surface zonal wind-stress (N/m2)
%  SIGMATHETA is the potential density (kg/m3)
%  RHO is the density (kg/m3)
%  EKL is the Ekman layer depth (m, positive)
%
% Files names are:
% INPUT:
% ./netcdf-files/<SNAPSHOT>/<netcdf_SIGMATHETA>.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/<netcdf_TAUX>.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/<netcdf_RHO>.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/<netcdf_EKL>.<netcdf_domain>.<netcdf_suff>
% OUTPUT:
% ./netcdf-files/<SNAPSHOT>/JFzx.<netcdf_domain>.<netcdf_suff>
% 
% with netcdf_* as global variables
%
% 06/04/12
% gmaze@mit.edu

function varargout = compute_JFzx(snapshot)

global sla toshow
global netcdf_suff netcdf_domain
global netcdf_TAUX netcdf_SIGMATHETA netcdf_EKL netcdf_RHO
pv_checkpath


% NETCDF file name:
filST  = netcdf_SIGMATHETA;
filTx  = netcdf_TAUX;
filRHO = netcdf_RHO;
filH   = netcdf_EKL;

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

% Load files:
ferfile = strcat(pathname,sla,snapshot,sla,filST,'.',netcdf_domain,'.',ext);
ncST     = netcdf(ferfile,'nowrite');
[STlon STlat STdpt] = coordfromnc(ncST);

ferfile = strcat(pathname,sla,snapshot,sla,filTx,'.',netcdf_domain,'.',ext);
ncTx    = netcdf(ferfile,'nowrite');

ferfile = strcat(pathname,sla,snapshot,sla,filRHO,'.',netcdf_domain,'.',ext);
ncRHO   = netcdf(ferfile,'nowrite');
RHO     = ncRHO{4}(1,:,:);

ferfile = strcat(pathname,sla,snapshot,sla,filH,'.',netcdf_domain,'.',ext);
ncH     = netcdf(ferfile,'nowrite');
EKL     = ncH{4}(1,:,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First term
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Dim:
if toshow, disp('dim'), end
nx = length(STlon) ;
ny = length(STlat) - 1 ;

% Pre-allocate:
if toshow, disp('pre-allocate'), end
dSIGMATHETAdy = zeros(ny-1,nx).*NaN;
dy        = zeros(1,ny).*NaN;
STup      = zeros(1,ny);
STdw      = zeros(1,ny);

% Meridional gradient of SIGMATHETA:
if toshow, disp('grad'), end
% Assuming the grid is regular, dy is independent of x:
dy = m_lldist([1 1]*STlon(1),STlat(1:ny+1) ) ; 
for ix = 1 : nx
  if toshow, disp(strcat(num2str(ix),'/',num2str(nx))), end
  STup  = squeeze(ncST{4}(1,2:ny+1,ix));
  STdw  = squeeze(ncST{4}(1,1:ny,ix));
  dSTdy = ( STup - STdw ) ./ dy;
  % Change horizontal grid point definition to fit with SIGMATHETA:
  dSTdy = ( dSTdy(1:ny-1) + dSTdy(2:ny) )./2; 
  dSIGMATHETAdy(:,ix) = dSTdy;
end %for iy

% Make TAUx having same limits:
TAUx = ncTx{4}(1,2:ny,:);

% Compute first term: TAUx * dSIGMATHETA/dy
JFz_a = TAUx .* dSIGMATHETAdy ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Finish ...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Then make all terms having same limits:
nx = length(STlon) ;
ny = length(STlat) ;
JFz_a   = squeeze(JFz_a(:,2:nx-1));
delta_e = squeeze(EKL(2:ny-1,2:nx-1));
rho     = squeeze(RHO(2:ny-1,2:nx-1));

% and finish:
JFz = JFz_a./delta_e./rho;


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

% General informations: 
netfil     = 'JFzx';
units      = 'kg/m3/s2';
ncid       = 'JFzx';
longname   = 'Vertical PV flux due to the zonal frictional force';
uniquename = 'JFzx';

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

% Define axis:
nx = length(STlon) ;
ny = length(STlat) ;
nz = 1 ;

nc('X') = nx-2;
nc('Y') = ny-2;
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'}(:)         = STlon(2:nx-1);
 
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'}(:)         = STlat(2:ny-1);
 
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'}(:)         = STdpt(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}(:,:,:)        = JFz;


%%% Close files:
close(ncST);
close(ncTx);
close(ncRHO);
close(ncH);
close(nc);

% Output:
output = struct('JFzx',JFz,'lat',STlat(2:ny-1),'lon',STlon(2:nx-1));
switch nargout
 case 1
  varargout(1) = {output};
end
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