Changeset 144

Timestamp:

04/01/15 18:20:38 (9 years ago)

Author:

jbrlod

Message:

adding feature for loading obs and filter design + dan_experiment

Location:

altifloat

Files:

• matlab_toolbox/make_filter.m (added)
• matlab_toolbox/script_perso.m (modified) (2 diffs)
• matlab_toolbox/script_prepare_exp_modele.m (modified) (1 diff)
• scripts/floater_model_assim.i (added)
• scripts/floater_model_forward.i (modified) (1 diff)
• src/floater.h (modified) (3 diffs)
• src/floater_delta.d (modified) (2 diffs)

altifloat/matlab_toolbox/script_perso.m

@@ -7 +7 @@
 %vt=load([ indir 'vzero.dat']);
 robs=load([ indir 'rfloat_total.dat']);
+rsim=load([ indir 'rfloat_total_jum0.dat']);
 uvb=load([ indir 'uv_bck_dan.dat']);
-uv=load([indir 'uv_total.dat']);
+uv=load([indir 'uv_total_jum0.dat']);
 meshg=load([indir 'meshgrid_dan.dat']);
 drifterdir='../altifloat/';
@@ -57 +58 @@
 
 for j=1:nfloat     
+  data= load([drifterdir fname{j}]);
+  
   rfloat=robs(robs(:,2)==j-1&robs(:,end)>0,[end-1 end]);
   rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1);
   rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1);
+  Hobs=plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'.-b');
+  
+  rfloat=rsim(rsim(:,2)==j-1&rsim(:,end)>0,[end-1 end]);
+  rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1);
+  rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1);
+  Hsim=plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'.-r');
+  
+
   %   plot(rfloat(:,2)+1,rfloat(:,1)+1,'.-m');
-  data= load([drifterdir fname{j}]);
-  Hobs=plot(data.lon(1:4),data.lat(1:4),'.-k')
+
+  Hreal=plot(data.lon(1:4),data.lat(1:4),'.-k')
   
-  Hsim=plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'.-m');
 end
 

altifloat/matlab_toolbox/script_prepare_exp_modele.m

@@ -39 +39 @@
 end % for j
 dlmwrite(['../exp_forw_dan/' initfile],MM,'\t');
+
+%STEP OBS : modify obs file
+

altifloat/scripts/floater_model_forward.i

@@ -16 +16 @@
 
 #Conditions initiales :
-#loadstate r_float 0 i 0 A 1 ../obs_float/init.dat D
+
 setstate r_float_d 0
 
-loadstate cfil 0 ij 0 A 1 ../obs_float/CC.dat D
-loadstate dfil 0 ij 0 A 1 ../obs_float/DD.dat D
-loadstate mask 0 ij 0 A 1 ../obs_float/Mask.dat D
-loadstate u_norm 0 ij 0 A 1 ../obs_float/NU.dat D
-loadstate v_norm 0 ij 0 A 1 ../obs_float/NV.dat D
+loadstate cfil 0 ij 0 A 1 ../exp_forw_dan/CC.dat D
+loadstate dfil 0 ij 0 A 1 ../exp_forw_dan/DD.dat D
+loadstate u_norm 0 ij 0 A 1 ../exp_forw_dan/NU.dat D
+loadstate v_norm 0 ij 0 A 1 ../exp_forw_dan/NV.dat D
 
 #Observations :

altifloat/src/floater.h

@@ -39 +39 @@
 //YREAL c2,c3,c4,c5;
 //YREAL dtfil=26042000;
-YREAL dtfil=2.604166666666667e+07;
-
+//YREAL dtfil=2.604166666666667e+07;
+YREAL dtfil=200000;
 #endif
 
@@ -213 +213 @@
         status=fscanf(fid,"%f",&val);
         mask[ilon][ilat]=(int)val;
+        YS1_mask(ilon,ilat)=mask[ilon][ilat];
 if (status!=1) {
         fprintf(stderr,"(read_mask) unable to read mask value (lon=%d,lat=%d)\n",ilon,ilat);
@@ -267 +268 @@
  
 }
- 
+
+void read_obs_2(int argv, char *argc[]) {
+  /* Read observation in a ascii file compatible with save_output_rfloat
+     format 
+     Column 1 (%d) : time step (from 0 to jptfl-1)
+     Column 2 (%d) : idfloat (from 0 to jpnfl-1)
+     Column 3 (%f) : grid point piret (lon)
+     Column 4 (%f) : grid point pjret (lat)
+  */
+  if (argv!=2) {
+    printf("(read_obs) Error : read_obs with %d (1 nedded)\n",argv-1);
+    printf("No obs read\n");
+    return;
+  }
+ FILE *fid;
+  int count=0;
+  fid=fopen(argc[1],"r");
+  if (fid==NULL) {
+    printf("(read_obs) Error : unable to open %s\n",argc[1]);
+    return;
+  }
+ int it,id;
+ YREAL posi,posj;
+  int status=4;
+  while(status==4) {
+    status=fscanf(fid,"%d %d %lf %lf",&it,&id,&posi,&posj);
+    if (status!=4) break;
+    count++;
+    piobs[it][id]=posi;
+    pjobs[it][id]=posj;
+    pmask[it][id]=1;
+  }//while (status==4)
+  fclose(fid);
+  fprintf(stdout,"(read_obs) %d data were loaded as observations\n",count);
+}
+
 void read_obs(int argv, char *argc[]) {
   /* Read observation in a ascii file

altifloat/src/floater_delta.d

@@ -18 +18 @@
 #ifdef FILTER
 defval K_FILTER //create a define K_FILTER in the source code
-defval OFIL 6  //order of the filter
+defval OFIL 5  //order of the filter
 #else
 defval OFIL 0 //Offtime of the main trajectory (0 since there is no filter)
@@ -237 +237 @@
 insert_fct arg run_inc
 insert_fct arg read_obs
+insert_fct arg read_obs_2
 insert_fct multirun
 insert_fct arg save_output_uv