MITgcm_contrib/gmaze_pv/compute_EKL.m

%
% [EKL] = compute_EKL(SNAPSHOT)
%
% Here we compute the Ekmal Layer Depth as:
% EKL = 0.7 sqrt( |TAU|/RHO )/f 
%
% where:
%  TAU is the amplitude of the 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>
% ./netcdf-files/<SNAPSHOT>/<netcdf_TAUY>.<netcdf_domain>.<netcdf_suff>
% OUTPUT
% ./netcdf-files/<SNAPSHOT>/<netcdf_EKL>.<netcdf_domain>.<netcdf_suff>
% 
% with netcdf_* as global variables
% netcdf_EKL = 'EKL' by default
%
% 08/16/06
% gmaze@mit.edu

function varargout = compute_EKL(snapshot)

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


% NETCDF file name:
filTx  = netcdf_TAUX;
filTy  = netcdf_TAUY;
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,filTy,'.',netcdf_domain,'.',ext);
ncTy    = netcdf(ferfile,'nowrite');
Ty      = ncTy{4}(1,:,:);
[Tylon Tylat Tydpt] = coordfromnc(ncTy);

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);
nz = 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 + Ty.^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_EKL')
  netfil = netcdf_EKL;
else
  netfil = 'EKL';
end
units      = 'm';
ncid       = 'EKL';
longname   = 'Ekman Layer Depth';
uniquename = 'EKL';

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

nc=close(nc);


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