Commit fa31f13743478bfe8eec1b9e7c81b000c73b90bd

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

En developpement ...

Showing 1 changed file with 120 additions and 67 deletions   Show diff stats
scripts/matlab/narr2gotm.m
... ... @@ -23,10 +23,28 @@ else
23 23 end
24 24  
25 25 % Time scale given in Matlab serial date number
  26 +ncfile = ['uwnd.narr.',name,'.',yrstr,'.nc'];
  27 +t = nc_varget(ncfile,'TIME');
  28 +tunits = nc_attget(ncfile,'TIME','units');
26 29 init_date = datenum([year 01 01 00 00 00]);
27   -time = [init_date:1/24:init_date+nod-1/24]';
  30 +tlength = length(t);
  31 +
  32 +delta_t = t(2)-t(1);
  33 +if strfind(tunits,'seconds');
  34 + time = t./86400;
  35 + time = time - time(1) + init_date;
  36 + timevec = datevec(time);
  37 +end
  38 +if strfind(tunits,'hours');
  39 + time = t./24;
  40 + time = time - time(1) + init_date;
  41 +end
  42 +if strfind(tunits,'days');
  43 + time = t;
  44 + time = time - time(1) + init_date;
  45 +end
  46 +
28 47 timevec = datevec(time);
29   -tlength = length(time);
30 48  
31 49 %% Horizontal coordinates
32 50  
... ... @@ -37,31 +55,29 @@ lat = nc_varget(ncfile,'LAT');
37 55 %% Choosing members
38 56  
39 57 % Indices of chosen locations
40   -i_loc = [5 5 6 6 6 7 7 7 8 8 8 9];
41   -j_loc = [4 5 3 4 5 3 4 5 3 4 5 3];
42   -nloc = length(i_loc);
  58 +iloc = [ 7 7 7 7 ...
  59 + 8 8 8 8 8 ...
  60 + 9 9 9 9 9 9 ...
  61 + 10 10 10 10 10 10 ...
  62 + 11 11 11 11 11 11 ...
  63 + 12 12 12 12 12 12 ...
  64 + 13 13 13 13 13];
  65 +
  66 +jloc = [ 8 9 10 11 ...
  67 + 7 8 9 10 11 ...
  68 + 6 7 8 9 10 11 ...
  69 + 6 7 8 9 10 11 ...
  70 + 5 6 7 8 9 10 ...
  71 + 5 6 7 8 9 10 ...
  72 + 5 6 7 8 9];
  73 +nloc = length(iloc);
43 74  
44 75 figure
45   -gray = [0.55 0.55 0.55];
46   -ax = worldmap([67 75],[-135 -115]);
47   -tightmap
48   -gridm('-')
49   -mlabel('on')
50   -plabel('on')
51   -geoshow('landareas.shp', 'FaceColor',gray)
52   -setm(gca,'FontName','Times','FontSize',12,'FontWeight','bold', ...
53   - 'MLabelLocation',[-135 -125], ...
54   - 'PLabelLocation',[70 75], ...
55   - 'PLineLocation',[70 75], ...
56   - 'PLineLimit', [-135 -115], ...
57   - 'MLineLocation',[-135 -130 -120], ...
58   - 'MLineLimit', [67 75], ...
59   - 'FLineWidth',1, ...
60   - 'GColor', [0.4 0.4 0.4])
61   -geoshow('worldlakes.shp', 'FaceColor', 'white');
62   -plotm(lat,lon,'ok')
63   -
64   -return
  76 +plot(lon,lat,'ok'); hold on
  77 +for n = 1:nloc
  78 + plot(lon(iloc(n),jloc(n)),lat(iloc(n),jloc(n)),'xk');
  79 +end
  80 +hold off
65 81  
66 82  
67 83 %% Momentum flux
... ... @@ -72,24 +88,30 @@ uwnd = nc_varget(ncfile,'UWND');
72 88 ncfile = ['vwnd.narr.',name,'.',yrstr,'.nc'];
73 89 vwnd = nc_varget(ncfile,'VWND');
74 90  
75   -% interpolation des donn???es ??? chaque heure
76   -uwnd = interp(uwnd,3);
77   -vwnd = interp(vwnd,3);
78 91  
79 92 % Wind stress calculation
80 93 % cf. Stewart, R.H. (2002) Introduction to Physical Oceanography
81 94 rho = 1.3; % density of air (kg m^-3)
82 95 U = sqrt(uwnd.^2 + vwnd.^2); % velocity amplitude (m s^-1)
83   -C10 = nan .* ones(length(U),1);
84   -for i = 1:length(U)
85   - if (U(i) >= 3 && U(i) < 6)
86   - C10(i) = (0.29 + 3.1./U(i) + 7.7./U(i).^2).*1e-3;
87   - elseif U(i) >= 6
88   - C10(i) = (0.60 + 0.070.*U(i)).*1e-3; % Yelland and Taylor (1996)
89   - else
90   - C10(i) = 2.18e-3;
  96 +C10 = nan .* ones(tlength,size(U,2),size(U,3));
  97 +for l = 1:tlength
  98 + for i = 1:size(U,2)
  99 + for j = 1:size(U,3)
  100 + if (U(l,i,j) >= 3 && U(l,i,j) < 6)
  101 + C10(l,i,j) = (0.29 + 3.1./U(l,i,j) + 7.7./U(l,i,j).^2).*1e-3;
  102 + elseif U(l,i,j) >= 6
  103 + C10(l,i,j) = (0.60 + 0.070.*U(l,i,j)).*1e-3; % Yelland and Taylor (1996)
  104 + else
  105 + C10(l,i,j) = 2.18e-3;
  106 + end
  107 + end
91 108 end
92 109 end
  110 +
  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 +
93 115 % C10 = (0.75 + 0.065.*U).*1e-3; % Garratt (1977)
94 116 dir = atan2(vwnd,uwnd); % wind direction (rad)
95 117 tau = rho.*C10.*U.*U; % wind stress amplitude (N m^-2)
... ... @@ -99,59 +121,90 @@ tauy = tau.*sin(dir); % meridional wind stress (N m^-2)
99 121 % figure; plot(U,tau,'ok')
100 122 clear U dir tau C10 rho
101 123  
102   -% Creation of the file 'narr_momentumflux.dat'
103   -% file with tx and ty given in N/m^2
104   -momentumflux = nan .* ones(tlength,8);
105   -momentumflux(:,1:6) = timevec;
106   -momentumflux(:,7) = taux;
107   -momentumflux(:,8) = tauy;
  124 +% Creation of nloc files 'narr_momentumflux_????_???.dat', one for each member
  125 +% tx and ty in N/m^2
108 126  
109   -fid = fopen(['narr_hourly_momentumflux_',yrstr,'.dat'], 'wt');
  127 +figure
  128 +for n = 1:nloc
  129 +
  130 + array = nan .* ones(tlength,8);
  131 + array(:,1:6) = timevec;
  132 + array(:,7) = squeeze(taux(:,iloc(n),jloc(n)));
  133 + array(:,8) = squeeze(tauy(:,iloc(n),jloc(n)));
  134 +
  135 + if n < 10
  136 + memstr = ['00',num2str(n)];
  137 + elseif n < 100
  138 + memstr = ['0',num2str(n)];
  139 + else
  140 + memstr = num2str(n);
  141 + 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 +end
  153 +
  154 +% Average file
  155 +
  156 +array = nan .* ones(tlength,8);
  157 +array(:,1:6) = timevec;
  158 +array(:,7) = mean(mean(taux(:,8:12,7:9),3),2);
  159 +array(:,8) = mean(mean(tauy(:,8:12,7:9),3),2);
  160 +
  161 +fid = fopen(['narr_momentumflux_',yrstr,'_000.dat'], 'wt');
110 162 fprintf(fid, '%4.0f-%02.0f-%02.0f %02.0f:%02.0f:%02.0f\t%f\t%f\n', ...
111   - momentumflux');
  163 + array');
112 164 fclose(fid);
113 165  
  166 +tau = sqrt(array(:,7).^2 + array(:,8).^2);
  167 +plot(time,tau,'-k');
114 168  
115   -% ----------------------------------------------------------------------
116   -% Incoming shortwave radiative flux and surface heat flux
  169 +hold off
  170 +
  171 +%% Incoming shortwave radiative flux and surface heat flux
117 172 % from NARR 3-hourly reanalysis
118 173 % ----------------------------------------------------------------------
119   -ncfile = ['dswrf_',yrstr,'_narr_amdgulf.nc'];
  174 +
  175 +ncfile = ['dswrf.narr.',name,'.',yrstr,'.nc'];
120 176 dswrf = nc_varget(ncfile,'DSWRF'); % W/m^2
121   -% interpolation des donn???es ??? chaque heure
122   -dswrf = interp(dswrf,3);
123   -a = find(dswrf < 0.0);
124   -dswrf(a) = 0.0;
  177 +% Remove negative values
  178 +dswrf(dswrf < 0.0) = 0.0;
125 179  
126   -ncfile = ['uswrf_',yrstr,'_narr_amdgulf.nc'];
  180 +ncfile = ['uswrf.narr.',name,'.',yrstr,'.nc'];
127 181 uswrf = nc_varget(ncfile,'USWRF'); % W/m^2
128   -% interpolation des donn???es ??? chaque heure
129   -uswrf = interp(uswrf,3);
130 182  
131   -ncfile = ['dlwrf_',yrstr,'_narr_amdgulf.nc'];
  183 +ncfile = ['dlwrf].narr.',name,'.',yrstr,'.nc'];
132 184 dlwrf = nc_varget(ncfile,'DLWRF'); % W/m^2
133   -% interpolation des donn???es ??? chaque heure
134   -dlwrf = interp(dlwrf,3);
135 185  
136   -ncfile = ['ulwrf_',yrstr,'_narr_amdgulf.nc'];
  186 +ncfile = ['ulwrf].narr.',name,'.',yrstr,'.nc'];
137 187 ulwrf = nc_varget(ncfile,'ULWRF'); % W/m^2
138   -% interpolation des donn???es ??? chaque heure
139   -ulwrf = interp(ulwrf,3);
140 188  
141   -ncfile = ['shtfl_',yrstr,'_narr_amdgulf.nc'];
  189 +ncfile = ['shtfl].narr.',name,'.',yrstr,'.nc'];
142 190 shtfl = nc_varget(ncfile,'SHTFL'); % W/m^2
143   -% interpolation des donn???es ??? chaque heure
144   -shtfl = interp(shtfl,3);
145 191  
146   -ncfile = ['lhtfl_',yrstr,'_narr_amdgulf.nc'];
  192 +ncfile = ['lhtfl].narr.',name,'.',yrstr,'.nc'];
147 193 lhtfl = nc_varget(ncfile,'LHTFL'); % W/m^2
148   -% interpolation des donn???es ??? chaque heure
149   -lhtfl = interp(lhtfl,3);
150 194  
151 195 % Transmitted shortwave radiative flux
152 196 % that include the extinction due to the ice cover
153   -% (Data from Canadian Ice Service)
  197 +% Canadian Ice Service
154 198 eval(['load cis_icec_',yrstr])
  199 +
  200 +% ECMWF
  201 +ncfile = ['ci.',yrstr,'.nc'];
  202 +icec = nc_varget(ncfile,'CI');
  203 +
  204 +% Indices over the Amundsen Gulf
  205 +iecmwf = [3 4 5];
  206 +jecmwf = [2 2 2];
  207 +
155 208 swr = (dswrf - uswrf).*(1.0-icec);
156 209 for ii = 1:length(swr)
157 210 if swr(ii) < 0.0;
... ...