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	gmaze	1.1	%
2	gmaze	1.2	% [JFzx] = compute_JFzx(SNAPSHOT)
3	gmaze	1.1	%
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	gmaze	1.2	function varargout = compute_JFzx(snapshot)
28	gmaze	1.1
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	gmaze	1.2
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