Changeset 160 for altifloat/matlab_toolbox

Timestamp:

07/20/15 11:51:06 (9 years ago)

Author:

jbrlod

Message:

finalize wind effect

Location:

altifloat/matlab_toolbox

Files:

• interp_wind.m (modified) (3 diffs)
• script_perso.m (modified) (4 diffs)

altifloat/matlab_toolbox/interp_wind.m

@@ -22 +22 @@
 GR_long=reshape(mg(:,1),Nlong,Nlat); 
 GR_lat=reshape(mg(:,2),Nlong,Nlat); 
+
+%Grid point size
+R_earth=6371229; %in meters
+delta{1}=zeros(Nlong,Nlat); %deltax
+delta{2}=zeros(Nlong,Nlat); %deltay (should correspond to fields order)
+
+for ilon=1:Nlong-1;
+    for jlat=1:Nlat-1;
+        
+    delta{1}(ilon,jlat)=R_earth*(2*pi/360)*(GR_long(ilon+1,jlat)-GR_long(ilon,jlat))...
+        .*cos(2*pi*GR_lat(ilon,jlat)/360); 
+    delta{2}(ilon,jlat)=R_earth*(2*pi/360)*(GR_lat(ilon,jlat+1)-GR_lat(ilon,jlat));  
+    
+    
+end 
+end
+
+
+
+delta{1}(Nlong,1:Nlat)=delta{1}(Nlong-1,1:Nlat);
+delta{2}(Nlong,1:Nlat)=delta{2}(Nlong-1,1:Nlat);
+
+
+delta{1}(1:Nlong,Nlat)=delta{1}(1:Nlong,Nlat-1);
+delta{2}(1:Nlong,Nlat)=delta{2}(1:Nlong,Nlat-1);
 
 
@@ -73 +98 @@
   data_fi = nan*ones(Nlong,Nlat,size(data,3));
   for it=1:size(data,3)
-  data_fi(:,:,it) = interp2(LAT,LON,data_f(:,:,it),GR_lat,GR_long,'nearest');
+  data_fi(:,:,it) = interp2(LAT,LON,data_f(:,:,it),GR_lat,GR_long,'nearest')./delta{k};
   end
   
@@ -80 +105 @@
   dimids = [londimid latdimid timeid];
   newid = netcdf.defVar(ncid,[fields{k} '_f'],'float',dimids);
-  netcdf.copyAtt(ncid,varid,'units',ncid,newid);
+  %netcdf.copyAtt(ncid,varid,'units',ncid,newid);
+  netcdf.putAtt(ncid,newid,'units','gridpoint per second');
   netcdf.copyAtt(ncid,varid,'long_name',ncid,newid);
   netcdf.putAtt(ncid,newid,'Description',['Convert to floating point and interpolated using' ...

altifloat/matlab_toolbox/script_perso.m

@@ -1 +1 @@
 close all
 clear all
-indir='../exp_dan_lr/';
+indir='../exp_hr/';
 %indir='../exp_twin_aviso/';
 %indir='../../code_leila/altifloat/obs_float/';
@@ -7 +7 @@
 %ut=load([ indir 'uzero.dat']);
 %vt=load([ indir 'vzero.dat']);
-robs=load([ indir 'obs_jum.dat']);
+robs=load([ indir 'obs_jum_cut_18.dat']);
 rsim=load([ indir 'rfloat_total_jum1.dat']);
-uvb=load([ indir 'uv_bck.dat']);
-uvt=load([ indir 'vel_true_all168_dan_sept1_sept7.dat']);
-uvr=load([indir 'uv_total_jum1.dat']);
-meshg=load([indir 'meshgrid_aviso.dat']);
+uvb=load([ indir 'uv_back_aviso_0901_0902_sp_t2.dat']);
+uvt=load([ indir 'uv_truthdan_0901_0903_lin_t2.dat']);
+uvr=load([indir 'uv_tot_jum1.dat']);
+uvgeo=load([indir 'uv_geos_jum1.dat']);
+meshg=load([indir 'meshgrid_aviso_t2.dat']);
 Div_ana_file=[indir 'Div_ana.dat'];
 Div_fg_file=[indir 'Div_fg.dat'];
@@ -61 +62 @@
 Vr=reshape(uvr(:,end),nlat,nlon,ntime);
 
+if exist('uvgeo')
+    ntime=length(unique(uvgeo(:,1)));
+Ugeo=reshape(uvgeo(:,end-1),nlat,nlon,ntime);
+Vgeo=reshape(uvgeo(:,end),nlat,nlon,ntime);
+end
+
 
 Hr=quiver(Lon(1:jech:end,1:iech:end),Lat(1:jech:end,1:iech:end),mean(Ur(1:jech:end,1:iech:end,1:end-1),3), ...
@@ -92 +99 @@
 Ht=quiver(Lon(1:jech:end,1:iech:end),Lat(1:jech:end,1:iech:end),mean(Ut(1:jech:end,1:iech:end,1:end-1),3),mean(Vt(1:jech:end,1:iech:end,1:end-1),3),'m');
 axis(LIM);
+