[170] | 1 | close all |
---|
| 2 | clear all |
---|
| 3 | |
---|
[191] | 4 | |
---|
| 5 | addpath('../matlab_toolbox'); |
---|
[170] | 6 | osave = true; |
---|
| 7 | |
---|
| 8 | indir='../exp_nemed/'; |
---|
| 9 | outdir = '../doc/ocean_modelling/fig/'; |
---|
| 10 | load ../../data/coastlines-split-4326/east_med.mat |
---|
| 11 | |
---|
| 12 | |
---|
| 13 | robs=load([ indir 'obs.dat']); |
---|
| 14 | rsim=load([ indir 'rfloat_nodiv.dat']); |
---|
[222] | 15 | %rsim=load([ indir 'rfloat_div2.dat']); |
---|
[170] | 16 | uvb=load([ indir 'uv_back_aviso.dat']); |
---|
| 17 | |
---|
| 18 | uvr=load([indir 'uv_tot_nodiv.dat']); |
---|
[222] | 19 | %uvr=load([indir 'uv_tot_div2.dat']); |
---|
[170] | 20 | meshg=load([indir 'meshgrid_aviso.dat']); |
---|
| 21 | mask = load([indir 'mask_aviso.dat']); |
---|
| 22 | |
---|
| 23 | iech=1; |
---|
| 24 | jech=1; |
---|
| 25 | %LIM=[34.5 36 32.5 35]; |
---|
| 26 | LON = [33.5 35.6]; |
---|
| 27 | LAT = [32.5 35]; |
---|
| 28 | LIM= [LON LAT]; |
---|
| 29 | |
---|
| 30 | %nlon=87; |
---|
| 31 | %nlat=58; |
---|
| 32 | lon=unique(meshg(:,1)); |
---|
| 33 | lat=unique(meshg(:,2)); |
---|
| 34 | nlon=length(lon); |
---|
| 35 | nlat=length(lat); |
---|
| 36 | |
---|
| 37 | [Lon,Lat]=meshgrid(lon',lat); |
---|
| 38 | Mask = reshape(mask,[nlon,nlat]); |
---|
| 39 | Mask = Mask'; |
---|
| 40 | |
---|
| 41 | nfloat=length(unique(robs(:,2))); |
---|
| 42 | |
---|
| 43 | %% To change unit |
---|
| 44 | %change grid from unitless to m/s |
---|
| 45 | |
---|
| 46 | R_earth=6371229; %in meters |
---|
| 47 | delta_x=zeros(nlat,nlon); |
---|
| 48 | delta_y=zeros(nlat,nlon); |
---|
| 49 | |
---|
| 50 | for ilon=1:nlon-1; |
---|
| 51 | for jlat=1:nlat-1; |
---|
| 52 | |
---|
| 53 | delta_x(jlat,ilon)=R_earth*(2*pi/360)*(Lon(jlat,ilon+1)-Lon(jlat,ilon))... |
---|
| 54 | .*cos(2*pi*Lat(jlat,ilon)/360); |
---|
| 55 | delta_y(jlat,ilon)=R_earth*(2*pi/360)*(Lat(jlat+1,ilon)-Lat(jlat,ilon)); |
---|
| 56 | |
---|
| 57 | end |
---|
| 58 | end |
---|
| 59 | |
---|
| 60 | |
---|
| 61 | |
---|
| 62 | delta_x(1:nlat,nlon)=delta_x(1:nlat,nlon); |
---|
| 63 | delta_y(1:nlat,nlon)=delta_y(1:nlat,nlon); |
---|
| 64 | |
---|
| 65 | |
---|
| 66 | delta_x(nlat,1:nlon)=delta_x(nlat-1,1:nlon); |
---|
| 67 | delta_y(nlat,1:nlon)=delta_y(nlat-1,1:nlon); |
---|
| 68 | |
---|
| 69 | |
---|
| 70 | %% Reshape wind data |
---|
| 71 | %it, ilon,ilat,u,v |
---|
| 72 | ntime=length(unique(uvb(:,1))); |
---|
| 73 | Ub=reshape(uvb(:,end-1),nlat,nlon,ntime); |
---|
| 74 | Vb=reshape(uvb(:,end),nlat,nlon,ntime); |
---|
| 75 | |
---|
| 76 | Ubms = Ub .* repmat(delta_x,1,1,ntime); |
---|
| 77 | Vbms = Vb .* repmat(delta_x,1,1,ntime); |
---|
| 78 | |
---|
| 79 | %it, ilon,ilat,u,v |
---|
| 80 | ntime=length(unique(uvr(:,1))); |
---|
| 81 | Ur=reshape(uvr(:,end-1),nlat,nlon,ntime); |
---|
| 82 | Vr=reshape(uvr(:,end),nlat,nlon,ntime); |
---|
| 83 | Ur(repmat(Mask,1,1,ntime)==0)=0; |
---|
| 84 | Vr(repmat(Mask,1,1,ntime)==0)=0; |
---|
| 85 | |
---|
| 86 | %Ur(Ub<1e-15)=0; |
---|
| 87 | %Vr(Vb<1e-15)=0; |
---|
| 88 | Urms = Ur .* repmat(delta_x,1,1,ntime); |
---|
| 89 | Vrms = Vr .* repmat(delta_x,1,1,ntime); |
---|
| 90 | |
---|
| 91 | %For the quiver legend |
---|
| 92 | [~,jleg] = min(abs(lon-35.3)); |
---|
| 93 | [~,ileg] = min(abs(lat-32.8)); |
---|
| 94 | Urms(ileg,jleg,:) = 1; |
---|
| 95 | %Ubms(ileg-1,jleg,:) = 1; |
---|
| 96 | |
---|
| 97 | %% Plots |
---|
| 98 | figure(1) |
---|
| 99 | clf |
---|
| 100 | Hr=quiver(Lon(1:jech:end,1:iech:end),Lat(1:jech:end,1:iech:end),mean(Urms(1:jech:end,1:iech:end,1:end-1),3), ... |
---|
| 101 | mean(Vrms(1:jech:end,1:iech:end,1:end-1),3),'r'); |
---|
| 102 | hold on |
---|
| 103 | Hb=quiver(Lon(1:jech:end,1:iech:end),Lat(1:jech:end,1:iech:end),mean(Ubms(1:jech:end,1:iech:end,1:end-1),3),mean(Vbms(1:jech:end,1:iech:end,1:end-1),3),'b'); |
---|
| 104 | |
---|
| 105 | |
---|
| 106 | |
---|
| 107 | |
---|
| 108 | for j=1:nfloat |
---|
| 109 | |
---|
| 110 | rfloat=rsim(rsim(:,2)==j-1&rsim(:,end)>0,[end-1 end]); |
---|
| 111 | rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1); |
---|
| 112 | rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1); |
---|
| 113 | %Hsim=plot(rfloat_lon,rfloat_lat,'.-r'); |
---|
| 114 | rfloat=robs(robs(:,2)==j-1&robs(:,end)>0,[end-1 end]); |
---|
| 115 | rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1); |
---|
| 116 | rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1); |
---|
| 117 | Hobs=plot(rfloat_lon,rfloat_lat,'m','Color',[0.3 0.3 0.3]); |
---|
| 118 | |
---|
| 119 | end |
---|
| 120 | |
---|
| 121 | |
---|
| 122 | |
---|
| 123 | |
---|
| 124 | |
---|
| 125 | |
---|
| 126 | %Ajout des trajectoires réelles |
---|
| 127 | axis(LIM); |
---|
[222] | 128 | mapshow(S,'Color',[0 0 0]); |
---|
| 129 | xlabel('Degrees East') |
---|
| 130 | ylabel('Degrees North') |
---|
[170] | 131 | |
---|
[191] | 132 | legend([Hb Hr Hobs],'Averaged over 9 days : background','Corrected','Drifter trajectory','Location','SouthWest'); |
---|
[170] | 133 | text(35.33,32.75,'1 m/s') |
---|
| 134 | %print -dpng ../exp_forw_dan/simu_drifter_dan.png |
---|
| 135 | |
---|
| 136 | if osave |
---|
| 137 | print('-dpng','-r300',[outdir 'Eddy_velocity.png']); |
---|
| 138 | savefig([outdir 'Eddy_velocity.fig']); |
---|
| 139 | end |
---|
| 140 | |
---|
| 141 | |
---|
| 142 | %% eddy |
---|
| 143 | cax = [-10 10].*1e-10; |
---|
| 144 | cmap = scol_dif_lin(100,cax); |
---|
| 145 | |
---|
| 146 | figure(2); |
---|
| 147 | clf |
---|
| 148 | [sX,sY] = meshgrid(34.3, 33:34.5); |
---|
| 149 | %streamline(Lon,Lat,mean(Ut(:,:,end-1),3),mean(Vt(:,:,end-1),3),sX,sY) |
---|
| 150 | [Uxb,Uyb]=gradient(mean(Ub(:,:,1:end-1),3),1,1); |
---|
| 151 | [Vxb,Vyb]=gradient(mean(Vb(:,:,1:end-1),3),1,1); |
---|
| 152 | |
---|
| 153 | %Uxt = Uxt .* (pi/180) .* R .* cos(Lat * pi/180); |
---|
| 154 | %Vxt = Vxt .* (pi/180) .* R .* cos(Lat * pi/180); |
---|
| 155 | %Uyt = Uyt .* (pi/180) .* R ; |
---|
| 156 | %Vyt = Vyt .* (pi/180) .* R ; |
---|
| 157 | |
---|
| 158 | |
---|
| 159 | owb = okubo_weiss(Uxb,Uyb,Vxb,Vyb); |
---|
| 160 | |
---|
| 161 | imagesc(lon',lat,owb(1:iech:end,1:iech:end)); |
---|
| 162 | caxis(cax); |
---|
| 163 | axis xy |
---|
| 164 | shading interp |
---|
| 165 | hold on |
---|
| 166 | mapshow(S,'Color',[0 0 0]); |
---|
| 167 | colorbar |
---|
| 168 | axis(LIM); |
---|
| 169 | colormap(cmap); |
---|
| 170 | [Uxr,Uyr]=gradient(mean(Ur(:,:,1:end-1),3),1,1); |
---|
| 171 | [Vxr,Vyr]=gradient(mean(Vr(:,:,1:end-1),3),1,1); |
---|
| 172 | |
---|
| 173 | %Uxr = Uxr .* (pi/180) .* R .* cos(Lat * pi/180); |
---|
| 174 | %Vxr = Vxr .* (pi/180) .* R .* cos(Lat * pi/180); |
---|
| 175 | %Uyr = Uyr .* (pi/180) .* R ; |
---|
| 176 | %Vyr = Vyr .* (pi/180) .* R ; |
---|
| 177 | |
---|
| 178 | owr = okubo_weiss(Uxr,Uyr,Vxr,Vyr); |
---|
| 179 | title('Okubo-Weiss parameter of the Background'); |
---|
[222] | 180 | xlabel('Degrees East') |
---|
| 181 | ylabel('Degrees North') |
---|
| 182 | |
---|
[170] | 183 | if osave |
---|
| 184 | print('-dpng','-r300',[outdir 'okubo_weiss_aviso.png']); |
---|
| 185 | savefig([outdir 'okubo_weiss_aviso.fig']); |
---|
| 186 | end |
---|
| 187 | |
---|
| 188 | figure(3) |
---|
| 189 | clf |
---|
| 190 | imagesc(lon',lat,owr(1:iech:end,1:iech:end)); |
---|
| 191 | caxis(cax); |
---|
| 192 | shading interp |
---|
| 193 | axis xy |
---|
| 194 | hold on |
---|
| 195 | mapshow(S,'Color',[0 0 0]) |
---|
| 196 | colorbar |
---|
| 197 | axis(LIM); |
---|
| 198 | colormap(cmap); |
---|
[191] | 199 | title('Okubo-Weiss parameter of the corrected field'); |
---|
[222] | 200 | xlabel('Degrees East') |
---|
| 201 | ylabel('Degrees North') |
---|
[170] | 202 | |
---|
| 203 | if osave |
---|
| 204 | print('-dpng','-r300',[outdir 'okubo_weiss_analyse.png']); |
---|
| 205 | savefig([outdir 'okubo_weiss_analyse.fig']); |
---|
| 206 | end |
---|
| 207 | |
---|