Commit e69ad1fa9d473cc484eaced026e06070d14d21aa

Authored by dumoda01
1 parent fa31f137
Exists in master and in 1 other branch snow

narr2gotm.m produit les fichiers momentum_flux. Reste a produire les heat_flux...

Showing 1 changed file with 68 additions and 48 deletions   Show diff stats
scripts/matlab/narr2gotm.m
... ... @@ -24,24 +24,26 @@ end
24 24  
25 25 % Time scale given in Matlab serial date number
26 26 ncfile = ['uwnd.narr.',name,'.',yrstr,'.nc'];
27   -t = nc_varget(ncfile,'TIME');
  27 +t1 = nc_varget(ncfile,'TIME');
28 28 tunits = nc_attget(ncfile,'TIME','units');
29 29 init_date = datenum([year 01 01 00 00 00]);
30   -tlength = length(t);
  30 +tlength = length(t1);
31 31  
32   -delta_t = t(2)-t(1);
  32 +delta_t = t1(2)-t1(1);
33 33 if strfind(tunits,'seconds');
34   - time = t./86400;
  34 + time = t1./86400;
35 35 time = time - time(1) + init_date;
36 36 timevec = datevec(time);
37 37 end
38 38 if strfind(tunits,'hours');
39   - time = t./24;
  39 + time = t1./24;
40 40 time = time - time(1) + init_date;
  41 + timevec = datevec(time);
41 42 end
42 43 if strfind(tunits,'days');
43   - time = t;
  44 + time = t1;
44 45 time = time - time(1) + init_date;
  46 + timevec = datevec(time);
45 47 end
46 48  
47 49 timevec = datevec(time);
... ... @@ -88,7 +90,6 @@ uwnd = nc_varget(ncfile,'UWND');
88 90 ncfile = ['vwnd.narr.',name,'.',yrstr,'.nc'];
89 91 vwnd = nc_varget(ncfile,'VWND');
90 92  
91   -
92 93 % Wind stress calculation
93 94 % cf. Stewart, R.H. (2002) Introduction to Physical Oceanography
94 95 rho = 1.3; % density of air (kg m^-3)
... ... @@ -108,11 +109,7 @@ for l = 1:tlength
108 109 end
109 110 end
110 111  
111   -% C10 = 2.18e-3 .* C10;
112   -% C10(U >=3 && U < 6) = (0.29 + 3.1./U + 7.7./U.^2).*1e-3;
113   -% C10(U >= 6) = (0.60 + 0.070.*U).*1e-3; % Yelland and Taylor (1996)
114   -
115   - % C10 = (0.75 + 0.065.*U).*1e-3; % Garratt (1977)
  112 +% C10 = (0.75 + 0.065.*U).*1e-3; % Garratt (1977)
116 113 dir = atan2(vwnd,uwnd); % wind direction (rad)
117 114 tau = rho.*C10.*U.*U; % wind stress amplitude (N m^-2)
118 115 taux = tau.*cos(dir); % zonal wind stress (N m^-2)
... ... @@ -132,27 +129,31 @@ for n = 1:nloc
132 129 array(:,7) = squeeze(taux(:,iloc(n),jloc(n)));
133 130 array(:,8) = squeeze(tauy(:,iloc(n),jloc(n)));
134 131  
135   - if n < 10
136   - memstr = ['00',num2str(n)];
137   - elseif n < 100
138   - memstr = ['0',num2str(n)];
139   - else
140   - memstr = num2str(n);
  132 + for m = 1:3
  133 +
  134 + mem = n + (m - 1).*nloc + m - 1;
  135 +
  136 + if mem < 10
  137 + memstr = ['00',num2str(mem)];
  138 + elseif n < 100
  139 + memstr = ['0',num2str(mem)];
  140 + else
  141 + memstr = num2str(mem);
  142 + end
  143 +
  144 + fid = fopen(['narr_momentumflux_',yrstr,'_',memstr,'.dat'], 'wt');
  145 + fprintf(fid, '%4.0f-%02.0f-%02.0f %02.0f:%02.0f:%02.0f\t%f\t%f\n', ...
  146 + array');
  147 + fclose(fid);
  148 +
  149 + tau = sqrt(array(:,7).^2 + array(:,8).^2);
  150 + plot(time,tau,'Color',[0.7 0.7 0.7]);
  151 + hold on
  152 +
141 153 end
142   -
143   - fid = fopen(['narr_momentumflux_',yrstr,'_',memstr,'.dat'], 'wt');
144   - fprintf(fid, '%4.0f-%02.0f-%02.0f %02.0f:%02.0f:%02.0f\t%f\t%f\n', ...
145   - array');
146   - fclose(fid);
147   -
148   - tau = sqrt(array(:,7).^2 + array(:,8).^2);
149   - plot(time,tau,'Color',[0.7 0.7 0.7]);
150   - hold on
151   -
152 154 end
153 155  
154   -% Average file
155   -
  156 +% Average file 000
156 157 array = nan .* ones(tlength,8);
157 158 array(:,1:6) = timevec;
158 159 array(:,7) = mean(mean(taux(:,8:12,7:9),3),2);
... ... @@ -165,12 +166,9 @@ fclose(fid);
165 166  
166 167 tau = sqrt(array(:,7).^2 + array(:,8).^2);
167 168 plot(time,tau,'-k');
168   -
169 169 hold off
170 170  
171 171 %% Incoming shortwave radiative flux and surface heat flux
172   -% from NARR 3-hourly reanalysis
173   -% ----------------------------------------------------------------------
174 172  
175 173 ncfile = ['dswrf.narr.',name,'.',yrstr,'.nc'];
176 174 dswrf = nc_varget(ncfile,'DSWRF'); % W/m^2
... ... @@ -180,30 +178,57 @@ dswrf(dswrf &lt; 0.0) = 0.0;
180 178 ncfile = ['uswrf.narr.',name,'.',yrstr,'.nc'];
181 179 uswrf = nc_varget(ncfile,'USWRF'); % W/m^2
182 180  
183   -ncfile = ['dlwrf].narr.',name,'.',yrstr,'.nc'];
  181 +ncfile = ['dlwrf.narr.',name,'.',yrstr,'.nc'];
184 182 dlwrf = nc_varget(ncfile,'DLWRF'); % W/m^2
185 183  
186   -ncfile = ['ulwrf].narr.',name,'.',yrstr,'.nc'];
  184 +ncfile = ['ulwrf.narr.',name,'.',yrstr,'.nc'];
187 185 ulwrf = nc_varget(ncfile,'ULWRF'); % W/m^2
188 186  
189   -ncfile = ['shtfl].narr.',name,'.',yrstr,'.nc'];
  187 +ncfile = ['shtfl.narr.',name,'.',yrstr,'.nc'];
190 188 shtfl = nc_varget(ncfile,'SHTFL'); % W/m^2
191 189  
192   -ncfile = ['lhtfl].narr.',name,'.',yrstr,'.nc'];
  190 +ncfile = ['lhtfl.narr.',name,'.',yrstr,'.nc'];
193 191 lhtfl = nc_varget(ncfile,'LHTFL'); % W/m^2
194 192  
195 193 % Transmitted shortwave radiative flux
196 194 % that include the extinction due to the ice cover
197 195 % Canadian Ice Service
198   -eval(['load cis_icec_',yrstr])
  196 +% eval(['load cis_icec_',yrstr])
199 197  
200 198 % ECMWF
201 199 ncfile = ['ci.',yrstr,'.nc'];
202   -icec = nc_varget(ncfile,'CI');
  200 +ci = nc_varget(ncfile,'ci');
  201 +
  202 +t2 = nc_varget(ncfile,'time');
  203 +t2 = double(t2);
  204 +tunits = nc_attget(ncfile,'time','units');
  205 +init_date = datenum([year 01 01 00 00 00]);
  206 +tlength = length(t2);
  207 +
  208 +if strfind(tunits,'seconds');
  209 + time2 = t2./86400;
  210 + time2 = time2 - time2(1) + init_date;
  211 + timevec = datevec(time2);
  212 +end
  213 +if strfind(tunits,'hours');
  214 + time2 = t2./24;
  215 + time2 = time2 - time2(1) + init_date;
  216 + timevec = datevec(time2);
  217 +end
  218 +if strfind(tunits,'days');
  219 + time2 = t2;
  220 + time2 = time2 - time2(1) + init_date;
  221 + timevec = datevec(time2);
  222 +end
203 223  
204 224 % Indices over the Amundsen Gulf
205   -iecmwf = [3 4 5];
206   -jecmwf = [2 2 2];
  225 +iecmwf = [2 2 2];
  226 +jecmwf = [3 4 5];
  227 +necmwf = length(iecmwf);
  228 +
  229 +for m = 1:necmwf
  230 + icec(:,m) = squeeze(ci(:,iecmwf(m),jecmwf(m)));
  231 +end
207 232  
208 233 swr = (dswrf - uswrf).*(1.0-icec);
209 234 for ii = 1:length(swr)
... ... @@ -211,17 +236,12 @@ for ii = 1:length(swr)
211 236 swr(ii) = 0.0;
212 237 end
213 238 end
214   -% No ice effect
215   -% swr = (dswrf - uswrf);
216 239  
217 240 % Total heat flux (minus the shortwave)
218 241 % positive = towards the ocean surface
219 242 % qt = dlwrf - ulwrf + shtfl + lhtfl; (????)
220   -% that include the insulation due to the ice cover
221   -% (Data from D. Barber group)
222   -qt = (dlwrf - ulwrf + shtfl + lhtfl) .* (1.0-icec);
223   -% No ice effect
224   -% qt = dlwrf - ulwrf + shtfl + lhtfl;
  243 +
  244 +qt_001 = (dlwrf - ulwrf + shtfl + lhtfl) .* (1.0-icec_001);
225 245  
226 246 figure; plot(time,qt)
227 247 figure; plot(time,swr)
... ...