/[MITgcm]/MITgcm_contrib/gmaze_pv/compute_EKLx.m
ViewVC logotype

Annotation of /MITgcm_contrib/gmaze_pv/compute_EKLx.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.4 - (hide annotations) (download)
Wed Sep 19 15:37:38 2007 UTC (17 years, 10 months ago) by gmaze
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +0 -0 lines
General Update

1 gmaze 1.1 %
2 gmaze 1.2 % [EKL] = compute_EKLx(SNAPSHOT)
3 gmaze 1.1 %
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 gmaze 1.2 function varargout = compute_EKLx(snapshot)
27 gmaze 1.1
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*(2*pi/86400)*sin(RHOlat*pi/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 gmaze 1.2
148     % Output:
149     output = struct('EKL',EKL,'lat',RHOlat,'lon',RHOlon);
150     switch nargout
151     case 1
152     varargout(1) = {output};
153     end

  ViewVC Help
Powered by ViewVC 1.1.22