MITgcm_contrib/gmaze_pv/compute_EKLx.m

%
% [EKL] = compute_EKLx(SNAPSHOT)
%
% Here we compute the Ekman Layer Depth as:
% EKL = 0.7 sqrt( TAUx/RHO )/f 
%
% where:
%  TAUx is the amplitude of the zonal surface wind-stress (N/m2)
%  RHO is the density of seawater (kg/m3)
%  f is the Coriolis parameter (kg/m3)
%  EKL is the Ekman layer depth (m)
%
% Files names are:
% INPUT:
% ./netcdf-files/<SNAPSHOT>/<netcdf_RHO>.<netcdf_domain>.<netcdf_suff>
% ./netcdf-files/<SNAPSHOT>/<netcdf_TAUX>.<netcdf_domain>.<netcdf_suff>
% OUTPUT
% ./netcdf-files/<SNAPSHOT>/<netcdf_EKLx>.<netcdf_domain>.<netcdf_suff>
% 
% with netcdf_* as global variables
% netcdf_EKLx = 'EKLx' by default
%
% 12/04/06
% gmaze@mit.edu

function varargout = compute_EKLx(snapshot)

global sla toshow
global netcdf_suff netcdf_domain
global netcdf_TAUX netcdf_RHO netcdf_EKLx
pv_checkpath
global EKL Tx Ty TAU RHO f


% NETCDF file name:
filTx  = netcdf_TAUX;
filRHO = netcdf_RHO;

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

% Load files:
ferfile = strcat(pathname,sla,snapshot,sla,filTx,'.',netcdf_domain,'.',ext);
ncTx    = netcdf(ferfile,'nowrite');
Tx      = ncTx{4}(1,:,:);

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get EKL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Dim:
if toshow, disp('dim'), end
nx = length(RHOlon);
ny = length(RHOlat);
ynz = length(RHOdpt);

% Pre-allocate:
if toshow, disp('pre-allocate'), end
EKL = zeros(ny,nx);

% Planetary vorticity:
f = 2*(2*pi/86400)*sin(RHOlat*pi/180);
[a f c]=meshgrid(RHOlon,f,RHOdpt); clear a c
f = permute(f,[3 1 2]);
f = squeeze(f(1,:,:));

% Windstress amplitude:
TAU = sqrt( Tx.^2 );

% Ekman Layer Depth:
EKL = 0.7* sqrt(TAU ./ RHO) ./f;
%EKL = 1.7975 * sqrt( TAU ./ RHO ./ f );

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

% General informations: 
if ~isempty('netcdf_EKLx')
  netfil = netcdf_EKLx;
else
  netfil = 'EKLx';
end
units      = 'm';
ncid       = 'EKLx';
longname   = 'Ekman Layer Depth from TAUx';
uniquename = 'EKLx';

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

% Define axis:
nx = length(RHOlon) ;
ny = length(RHOlat) ;
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'}(:)         = RHOlon;
 
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'}(:)         = RHOlat;
 
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'}(:)         = RHOdpt(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}(:,:,:)        = EKL;


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


% Output:
output = struct('EKL',EKL,'lat',RHOlat,'lon',RHOlon);
switch nargout
 case 1
  varargout(1) = {output};
end
1	%
2	% [EKL] = compute_EKLx(SNAPSHOT)
3	%
4	% Here we compute the Ekman Layer Depth as:
5	% EKL = 0.7 sqrt( TAUx/RHO )/f
6	%
7	% where:
8	% TAUx is the amplitude of the zonal surface wind-stress (N/m2)
9	% RHO is the density of seawater (kg/m3)
10	% f is the Coriolis parameter (kg/m3)
11	% EKL is the Ekman layer depth (m)
12	%
13	% Files names are:
14	% INPUT:
15	% ./netcdf-files/<SNAPSHOT>/<netcdf_RHO>.<netcdf_domain>.<netcdf_suff>
16	% ./netcdf-files/<SNAPSHOT>/<netcdf_TAUX>.<netcdf_domain>.<netcdf_suff>
17	% OUTPUT
18	% ./netcdf-files/<SNAPSHOT>/<netcdf_EKLx>.<netcdf_domain>.<netcdf_suff>
19	%
20	% with netcdf_* as global variables
21	% netcdf_EKLx = 'EKLx' by default
22	%
23	% 12/04/06
24	% gmaze@mit.edu
25
26	function varargout = compute_EKLx(snapshot)
27
28	global sla toshow
29	global netcdf_suff netcdf_domain
30	global netcdf_TAUX netcdf_RHO netcdf_EKLx
31	pv_checkpath
32	global EKL Tx Ty TAU RHO f
33
34
35	% NETCDF file name:
36	filTx = netcdf_TAUX;
37	filRHO = netcdf_RHO;
38
39	% Path and extension to find them:
40	pathname = strcat('netcdf-files',sla);
41	ext = netcdf_suff;
42
43	% Load files:
44	ferfile = strcat(pathname,sla,snapshot,sla,filTx,'.',netcdf_domain,'.',ext);
45	ncTx = netcdf(ferfile,'nowrite');
46	Tx = ncTx{4}(1,:,:);
47
48	ferfile = strcat(pathname,sla,snapshot,sla,filRHO,'.',netcdf_domain,'.',ext);
49	ncRHO = netcdf(ferfile,'nowrite');
50	RHO = ncRHO{4}(1,:,:);
51	[RHOlon RHOlat RHOdpt] = coordfromnc(ncRHO);
52
53
54	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
55	% Get EKL
56	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57
58	% Dim:
59	if toshow, disp('dim'), end
60	nx = length(RHOlon);
61	ny = length(RHOlat);
62	ynz = length(RHOdpt);
63
64	% Pre-allocate:
65	if toshow, disp('pre-allocate'), end
66	EKL = zeros(ny,nx);
67
68	% Planetary vorticity:
69	f = 2(2pi/86400)sin(RHOlatpi/180);
70	[a f c]=meshgrid(RHOlon,f,RHOdpt); clear a c
71	f = permute(f,[3 1 2]);
72	f = squeeze(f(1,:,:));
73
74	% Windstress amplitude:
75	TAU = sqrt( Tx.^2 );
76
77	% Ekman Layer Depth:
78	EKL = 0.7* sqrt(TAU ./ RHO) ./f;
79	%EKL = 1.7975 * sqrt( TAU ./ RHO ./ f );
80
81	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
82	% Record
83	%%%%%%%%%%%%%%%%%%%%%%%%%%%%
84	if toshow, disp('record'), end
85
86	% General informations:
87	if ~isempty('netcdf_EKLx')
88	netfil = netcdf_EKLx;
89	else
90	netfil = 'EKLx';
91	end
92	units = 'm';
93	ncid = 'EKLx';
94	longname = 'Ekman Layer Depth from TAUx';
95	uniquename = 'EKLx';
96
97	% Open output file:
98	nc = netcdf(strcat(pathname,sla,snapshot,sla,netfil,'.',netcdf_domain,'.',ext),'clobber');
99
100	% Define axis:
101	nx = length(RHOlon) ;
102	ny = length(RHOlat) ;
103	nz = 1 ;
104
105	nc('X') = nx;
106	nc('Y') = ny;
107	nc('Z') = nz;
108
109	nc{'X'} = ncfloat('X');
110	nc{'X'}.uniquename = ncchar('X');
111	nc{'X'}.long_name = ncchar('longitude');
112	nc{'X'}.gridtype = nclong(0);
113	nc{'X'}.units = ncchar('degrees_east');
114	nc{'X'}(:) = RHOlon;
115
116	nc{'Y'} = ncfloat('Y');
117	nc{'Y'}.uniquename = ncchar('Y');
118	nc{'Y'}.long_name = ncchar('latitude');
119	nc{'Y'}.gridtype = nclong(0);
120	nc{'Y'}.units = ncchar('degrees_north');
121	nc{'Y'}(:) = RHOlat;
122
123	nc{'Z'} = ncfloat('Z');
124	nc{'Z'}.uniquename = ncchar('Z');
125	nc{'Z'}.long_name = ncchar('depth');
126	nc{'Z'}.gridtype = nclong(0);
127	nc{'Z'}.units = ncchar('m');
128	nc{'Z'}(:) = RHOdpt(1);
129
130	% And main field:
131	nc{ncid} = ncfloat('Z', 'Y', 'X');
132	nc{ncid}.units = ncchar(units);
133	nc{ncid}.missing_value = ncfloat(NaN);
134	nc{ncid}.FillValue_ = ncfloat(NaN);
135	nc{ncid}.longname = ncchar(longname);
136	nc{ncid}.uniquename = ncchar(uniquename);
137	nc{ncid}(:,:,:) = EKL;
138
139
140
141	% Close files:
142	close(ncTx);
143	close(ncRHO);
144	close(nc);
145
146
147
148	% Output:
149	output = struct('EKL',EKL,'lat',RHOlat,'lon',RHOlon);
150	switch nargout
151	case 1
152	varargout(1) = {output};
153	end