MITgcm_contrib/gmaze_pv/compute_JFz.m

%
% [] = compute_JFz(SNAPSHOT)
%
% Here we compute the PV flux due to frictionnal forces as
% JFz = ( TAUx * dSIGMATHETA/dy - TAUy * dSIGMATHETA/dx ) / RHO / EKL
%
% where:
%  TAU is the surface 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_TAUY>.<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>/JFz.<netcdf_domain>.<netcdf_suff>
% 
% with netcdf_* as global variables
%
% 06/27/06
% gmaze@mit.edu

function compute_JFz(snapshot)

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


% NETCDF file name:
filST  = netcdf_SIGMATHETA;
filTx  = netcdf_TAUX;
filTy  = netcdf_TAUY;
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,filTy,'.',netcdf_domain,'.',ext);
ncTy    = 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 ;
nz = length(STdpt);

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

% Meridional gradient of SIGMATHETA:
if toshow, disp('grad'), end
% Assuming the grid is regular, dy is independent of x:
[dy b] = meshgrid( m_lldist([1 1]*STlon(1),STlat(1:ny+1) ), STdpt ) ; clear b
for ix = 1 : nx
  if toshow, disp(strcat(num2str(ix),'/',num2str(nx))), end
  STup  = squeeze(ncST{4}(:,2:ny+1,ix));
  STdw  = squeeze(ncST{4}(:,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
iz    = 1;
JFz_a = TAUx .* squeeze(dSIGMATHETAdy(iz,:,:)) ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Second term
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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

% Zonal gradient of SIGMATHETA
if toshow, disp('grad'), end
for iy = 1 : ny
  if toshow, disp(strcat(num2str(iy),'/',num2str(ny))), end
  [dx b] = meshgrid( m_lldist(STlon(1:nx+1),[1 1]*STlat(iy)), STdpt ) ; clear b
  STup    = squeeze(ncST{4}(:,iy,2:nx+1));
  STdw    = squeeze(ncST{4}(:,iy,1:nx));
  dSTdx   = ( STup - STdw ) ./ dx;
  % Change horizontal grid point definition to fit with SIGMATHETA:
  dSTdx   = ( dSTdx(:,1:nx-1) + dSTdx(:,2:nx) )./2;
  dSIGMATHETAdx(:,iy,:) = dSTdx;
end %for iy

% Make TAUy having same limits:
TAUy  = ncTy{4}(1,:,2:nx);

% Compute second term: TAUy * dSIGMATHETA/dx
iz    = 1;
JFz_b = TAUy .* squeeze(dSIGMATHETAdx(iz,:,:)) ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Finish ...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Then make all terms having same limits:
nx = length(STlon) ;
ny = length(STlat) ;
nz = length(STdpt) ;
JFz_a   = squeeze(JFz_a(:,2:nx-1));
JFz_b   = squeeze(JFz_b(2:ny-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 - JFz_b)./delta_e./rho;


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

% General informations: 
netfil     = 'JFz';
units      = 'kg/m3/s2';
ncid       = 'JFz';
longname   = 'Vertical PV flux due to frictional forces';
uniquename = 'JFz';

% 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;

nc=close(nc);


1	gmaze	1.1	%
2			% [] = 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	gmaze	1.2	% EKL is the Ekman layer depth (m, positive)
12	gmaze	1.1	%
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 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