MITgcm_contrib/gmaze_pv/compute_JBz.m

%
% [JBz] = compute_JBz(SNAPSHOT)
%
% Here we compute the PV flux due to diabatic processes as
% JFz = - alpha * f * Qnet / MLD / Cw
% where:
%  alpha = 2.5*E-4 1/K is the thermal expansion coefficient
%  f = 2*OMEGA*sin(LAT) is the Coriolis parameter
%  Qnet is the net surface heat flux (W/m^2), positive downward
%  MLD is the mixed layer depth (m, positive)
%  Cw = 4187 J/kg/K is the specific heat of seawater
%
% Files names are:
% INPUT:
% ./netcdf-files/<SNAPSHOT>/<netcdf_Qnet>.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/<netcdf_MLD>.<netcdf_domain>.<netcdf_suff>
% OUTPUT:
% ./netcdf-files/<SNAPSHOT>/JBz.<netcdf_domain>.<netcdf_suff>
% 
% with: netcdf_* as global variables
%
% 06/27/06
% gmaze@mit.edu

function varargout = compute_JBz(snapshot)

global sla toshow
global netcdf_suff netcdf_domain
global netcdf_Qnet netcdf_MLD
pv_checkpath


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

% Load files:
ferfile = strcat(pathname,sla,snapshot,sla,netcdf_Qnet,'.',netcdf_domain,'.',ext);
ncQ     = netcdf(ferfile,'nowrite');
[Qlon Qlat Qdpt] = coordfromnc(ncQ);

ferfile = strcat(pathname,sla,snapshot,sla,netcdf_MLD,'.',netcdf_domain,'.',ext);
ncMLD   = netcdf(ferfile,'nowrite');
[MLDlon MLDlat MLDdpt] = coordfromnc(ncMLD);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% surface PV flux
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define axis:
nx = length(Qlon) ;
ny = length(Qlat) ;
nz = length(Qdpt) ;


% Planetary vorticity:
f     = 2*(2*pi/86400)*sin(Qlat*pi/180);
[a f] = meshgrid(Qlon,f); clear a c


% Net surface heat flux:
Qnet = ncQ{4}(:,:,:);


% Mixed layer Depth:
MLD = ncMLD{4}(:,:,:);


% Coefficient:
alpha = 2.5*10^(-4); % Surface average value
Cw    = 4187;             
coef  = - alpha / Cw;


% JBz:
JBz = zeros(nz,ny,nx).*NaN;
JBz(1,:,:) = coef*f.*Qnet./MLD;


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

% General informations: 
netfil     = 'JBz';
units      = 'kg/m3/s2';
ncid       = 'JBz';
longname   = 'Vertical PV flux due to diabatic processes';
uniquename = 'JBz';

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

% Define axis:
nx = length(Qlon) ;
ny = length(Qlat) ;
nz = 1 ;

nc('X') = nx;
nc('Y') = ny;
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'}(:)         = Qlon;
 
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'}(:)         = Qlat;
 
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'}(:)         = Qdpt(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}(:,:,:)        = JBz;

nc=close(nc);
close(ncQ);
close(ncMLD);


% Output:
output = struct('JBz',JBz,'lat',Qlat,'lon',Qlon);
switch nargout
 case 1
  varargout(1) = {output};
end
1	%
2	% [JBz] = compute_JBz(SNAPSHOT)
3	%
4	% Here we compute the PV flux due to diabatic processes as
5	% JFz = - alpha * f * Qnet / MLD / Cw
6	% where:
7	% alpha = 2.5*E-4 1/K is the thermal expansion coefficient
8	% f = 2OMEGAsin(LAT) is the Coriolis parameter
9	% Qnet is the net surface heat flux (W/m^2), positive downward
10	% MLD is the mixed layer depth (m, positive)
11	% Cw = 4187 J/kg/K is the specific heat of seawater
12	%
13	% Files names are:
14	% INPUT:
15	% ./netcdf-files/<SNAPSHOT>/<netcdf_Qnet>.<netcdf_domain>.<netcdf_suff>
16	% ./netcdf-files/<SNAPSHOT>/<netcdf_MLD>.<netcdf_domain>.<netcdf_suff>
17	% OUTPUT:
18	% ./netcdf-files/<SNAPSHOT>/JBz.<netcdf_domain>.<netcdf_suff>
19	%
20	% with: netcdf_* as global variables
21	%
22	% 06/27/06
23	% gmaze@mit.edu
24
25	function varargout = compute_JBz(snapshot)
26
27	global sla toshow
28	global netcdf_suff netcdf_domain
29	global netcdf_Qnet netcdf_MLD
30	pv_checkpath
31
32
33	% Path and extension to find netcdf-files:
34	pathname = strcat('netcdf-files',sla);
35	ext = netcdf_suff;
36
37	% Load files:
38	ferfile = strcat(pathname,sla,snapshot,sla,netcdf_Qnet,'.',netcdf_domain,'.',ext);
39	ncQ = netcdf(ferfile,'nowrite');
40	[Qlon Qlat Qdpt] = coordfromnc(ncQ);
41
42	ferfile = strcat(pathname,sla,snapshot,sla,netcdf_MLD,'.',netcdf_domain,'.',ext);
43	ncMLD = netcdf(ferfile,'nowrite');
44	[MLDlon MLDlat MLDdpt] = coordfromnc(ncMLD);
45
46
47
48	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49	% surface PV flux
50	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51
52	% Define axis:
53	nx = length(Qlon) ;
54	ny = length(Qlat) ;
55	nz = length(Qdpt) ;
56
57
58	% Planetary vorticity:
59	f = 2(2pi/86400)sin(Qlatpi/180);
60	[a f] = meshgrid(Qlon,f); clear a c
61
62
63	% Net surface heat flux:
64	Qnet = ncQ{4}(:,:,:);
65
66
67	% Mixed layer Depth:
68	MLD = ncMLD{4}(:,:,:);
69
70
71	% Coefficient:
72	alpha = 2.5*10^(-4); % Surface average value
73	Cw = 4187;
74	coef = - alpha / Cw;
75
76
77	% JBz:
78	JBz = zeros(nz,ny,nx).*NaN;
79	JBz(1,:,:) = coeff.Qnet./MLD;
80
81
82
83	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
84	% Record
85	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
86	if toshow, disp('record'), end
87
88	% General informations:
89	netfil = 'JBz';
90	units = 'kg/m3/s2';
91	ncid = 'JBz';
92	longname = 'Vertical PV flux due to diabatic processes';
93	uniquename = 'JBz';
94
95	% Open output file:
96	nc = netcdf(strcat(pathname,sla,snapshot,sla,netfil,'.',netcdf_domain,'.',ext),'clobber');
97
98	% Define axis:
99	nx = length(Qlon) ;
100	ny = length(Qlat) ;
101	nz = 1 ;
102
103	nc('X') = nx;
104	nc('Y') = ny;
105	nc('Z') = nz;
106
107	nc{'X'} = ncfloat('X');
108	nc{'X'}.uniquename = ncchar('X');
109	nc{'X'}.long_name = ncchar('longitude');
110	nc{'X'}.gridtype = nclong(0);
111	nc{'X'}.units = ncchar('degrees_east');
112	nc{'X'}(:) = Qlon;
113
114	nc{'Y'} = ncfloat('Y');
115	nc{'Y'}.uniquename = ncchar('Y');
116	nc{'Y'}.long_name = ncchar('latitude');
117	nc{'Y'}.gridtype = nclong(0);
118	nc{'Y'}.units = ncchar('degrees_north');
119	nc{'Y'}(:) = Qlat;
120
121	nc{'Z'} = ncfloat('Z');
122	nc{'Z'}.uniquename = ncchar('Z');
123	nc{'Z'}.long_name = ncchar('depth');
124	nc{'Z'}.gridtype = nclong(0);
125	nc{'Z'}.units = ncchar('m');
126	nc{'Z'}(:) = Qdpt(1);
127
128	% And main field:
129	nc{ncid} = ncfloat('Z', 'Y', 'X');
130	nc{ncid}.units = ncchar(units);
131	nc{ncid}.missing_value = ncfloat(NaN);
132	nc{ncid}.FillValue_ = ncfloat(NaN);
133	nc{ncid}.longname = ncchar(longname);
134	nc{ncid}.uniquename = ncchar(uniquename);
135	nc{ncid}(:,:,:) = JBz;
136
137	nc=close(nc);
138	close(ncQ);
139	close(ncMLD);
140
141
142
143	% Output:
144	output = struct('JBz',JBz,'lat',Qlat,'lon',Qlon);
145	switch nargout
146	case 1
147	varargout(1) = {output};
148	end