source: altifloat/src/floater.h @ 160

Last change on this file since 160 was 160, checked in by jbrlod, 9 years ago

finalize wind effect
File size: 21.6 KB
Line 
1#define UPDATE_BCK
2
3#include <netcdf.h>
4
5//constantes
6//YREAL rdtflo=2*3600;
7YREAL rdtflo=3600;
8//YREAL rdtflo=1;
9YREAL ur_dx=0;
10YREAL ur_dy=0;
11YREAL vr_dx=0;
12YREAL vr_dy=0;
13
14
15//Mode Forward ou assim;
16typedef enum {ASSIM,FORW} assim_t;
17
18assim_t tassim = ASSIM;
19
20//Tableaux
21
22//Variable d'état
23YREAL umod[nlon][nlat][jptfl];
24YREAL vmod[nlon][nlat][jptfl];
25YREAL uwind[nlon][nlat][jptfl];
26YREAL vwind[nlon][nlat][jptfl];
27YREAL utotal[nlon][nlat][jptfl];
28YREAL vtotal[nlon][nlat][jptfl];
29
30//Background
31YREAL ubck[nlon][nlat][jptfl];
32YREAL vbck[nlon][nlat][jptfl];
33int mask[nlon][nlat];
34
35//Trajectoires des flotteurs (en points de grille)
36YREAL piobs[jptfl][jpnfl];
37YREAL pjobs[jptfl][jpnfl];
38YREAL piret[jptfl][jpnfl];
39YREAL pjret[jptfl][jpnfl];
40short pmask[jptfl][jpnfl];
41
42//FILTER COEFFICIENTS
43#ifdef K_FILTER
44//YREAL c1;
45//YREAL c2,c3,c4,c5;
46//YREAL dtfil=26042000;
47//YREAL dtfil=2.604166666666667e+07;
48
49#ifdef INC_FIL
50#include "fil_def.h"
51#else
52YREAL dtfil=200000;
53#endif
54#endif
55
56#define ERR(e) {fprintf(stderr,"Error: %s\n", nc_strerror(e)); exit(EXIT_FAILURE);}
57#define LONDIM "lon_interp"
58#define LATDIM "lat_interp"
59#define UVAR   "u10_f"
60#define VVAR   "v10_f"
61
62//Options du run incr士ental
63int nb_extiter=8;
64short inc_save=0;
65char dirsave[50]="../obs_float/";
66char savename[50]="optim25float-filtre15-inc";
67extern double Y3ddf1;
68//Declaration pour multirun
69void erase_lobs();
70short is_activ(int jpfl, int it);
71void read_obs(int argv, char *argc[]);
72extern int Yobs_insert_data (char *nmmod, int sortie, int iaxe, int jaxe, int kaxe,int pdt, YREAL val);
73void run_inc(int argv, char *argc[]);
74void div_init();
75void init_wind(char fname[],int iter0,int *ptime0);
76
77void appli_start(int argc, char *argv[]){
78#ifdef K_FILTER
79  //c1=0.5;
80  //c2=c3=c4=c5
81#endif
82
83  //Init des tableaux
84  for (int i=0;i<nlon;i++)
85    for (int j=0;j<nlat;j++) {
86        mask[i][j]=1;   
87      for (int k=0;k<jptfl;k++) {
88        umod[i][j][k]=0;
89        vmod[i][j][k]=0;
90        ubck[i][j][k]=0;
91        vbck[i][j][k]=0;
92        utotal[i][j][k]=0;
93        vtotal[i][j][k]=0;
94      }
95}
96 for (int i=0;i<jptfl;i++)
97    for (int j=0;j<jpnfl;j++) {
98      piobs[i][j]=0;
99      pjobs[i][j]=0;
100      piret[i][j]=0;
101      pjret[i][j]=0; 
102      pmask[i][j]=0;
103      }
104
105
106 /*
107char L1[50]="read_obs";
108char L2[50]="../obs_float/obs_float_test.dat";
109char *L[2];
110 L[0]=L1;
111 L[1]=L2;
112 read_obs(2,L);*/
113 //Test de isactiv
114 //printf("isactiv(0,0):%d\n",is_activ(0,0));
115 //printf("isactiv(1,0):%d\n",is_activ(1,0));
116 //printf("isactiv(0,3):%d\n",is_activ(0,3));
117 //printf("isactiv(1,3):%d\n",is_activ(1,3));
118 //printf("isactiv(0,6):%d\n",is_activ(0,6));
119 //printf("isactiv(1,6):%d\n",is_activ(1,6));
120
121 
122 // erase_lobs();
123}
124void before_it(int nit){}
125void cost_function(int pdt){}
126void adjust_target(){}
127void after_it(int nit){
128       
129        //printf("u_d(9,16)=%f %f %f %f\n",YS1_u_d(9,16),YS1_u_d(10,16),YS1_u_d(10,17),YS1_u_d(9,17));
130        }
131void forward_before(int ctrp){
132}
133void forward_after(int ctrp){
134
135}
136
137void backward_before(int ctrp){
138}
139
140void backward_after(int ctrp){
141        }
142
143short select_io(int indic, char *nmod, int sortie, int iaxe,
144                int jaxe, int kaxe, int pdt, YREAL *val){
145        //if(indic==YIO_LOADSTATE && !strcmp(nmod,"u")) {
146        //printf("s=%d t=%d Yi=%d Yj=%d val=%f\n",sortie,pdt,iaxe,jaxe,*val);   
147        //}
148       
149        return(1);
150}
151void clear_Yst_nodo(struct Yst_nodo *nobs, int lev, int max) {
152  /* lev est le niveau courant (racine=0)
153max est le niveau des feuilles soit
1544 si dimod=1
1555 si dimod=2
1566 si dimod=3
157dimod est défini au niveau 2 avec dimod=YTabMod[imod].dim
158
159
160  */
161 
162  if (lev==max) {
163    //on nettoie une feuille (ni frere ni fils)
164    free(nobs);
165    nobs = NULL ;
166  }
167 
168  else {
169
170    if (lev==2) {
171      //on peut déterminer le niveau max
172      int dimod=YTabMod[nobs->iind].dim;
173      max=dimod+3;
174    }
175
176    if (nobs->fils!=NULL && nobs->frere!=NULL) {
177      clear_Yst_nodo(nobs->frere,lev,max);
178      nobs->frere=NULL;
179    }   
180    if (nobs->fils!=NULL) {
181      clear_Yst_nodo(nobs->fils,lev+1,max);
182      nobs->fils=NULL;
183    }
184    free(nobs);
185    nobs=NULL;
186  } // else (pas feuille)
187 
188
189}
190
191//in case we want to see how it s loading the data we do this above
192void erase_lobs() {
193  //Efface toute l'arboresence des obs (y compris les ébauches)
194 
195  int itraj;
196  struct Yst_nodo *YRobs = NULL;
197  YioInsertObsCtr = 0 ;
198  for (itraj=0;itraj<YNBTRAJ;itraj++) {
199    //Effacement des obs liées au modul imod
200
201    //  if (YTabMod[itraj].is_target || !YTabMod[itraj].is_cout) {
202      YRobs=YTabTraj[itraj].YRobs;
203      clear_Yst_nodo(YRobs,0,10);
204
205      YTabTraj[itraj].YRobs=NULL;
206      //}
207
208  } //for itraj
209  for (int wi = 0; wi < YNBTRAJ; ++wi)
210    {
211      YTabTraj[wi].YRobs = (struct Yst_nodo*) malloc (sizeof (struct Yst_nodo)); /* init du Root des arborescences */
212      YTabTraj[wi].YRobs->frere = YTabTraj[wi].YRobs->fils = NULL;      /* d'observations pour chaque trajectoire */
213    }
214
215}
216
217
218void load_wind(int argv, char *argc[]) {
219  /* load_wind windname iter0 time0
220     load a wind netcdfile ( from ERAI)
221     time0 is the time of the first iteration considered in the file
222     (by default time0 is time(0) field in the netcdf file)
223     iter0 is the timestep corresponding to time0.
224     (by default iter0 is 0)
225   */
226  if ((argv<2) | (argv >4)) {
227    fprintf(stderr,"(load_wind): wrong number of arguments\n");
228    return;
229  }
230
231  //time0 is a pointer to pass the NULL address in case time0 is
232  //not specified in the call of the function
233  int *time0=NULL;
234
235  int iter0=0;
236  char fname[50];
237  sprintf(fname,argc[1]);
238  if (argv > 2) {
239    iter0 = atoi(argc[2]);
240    if (argv > 3) 
241      *time0 = atoi(argc[3]);
242  }
243  init_wind(fname,iter0,time0);
244}
245
246void init_wind(char fname[],int iter0,int *ptime0) {
247  int ncid,retval;
248  int dlatid,dlonid,dtimeid;
249  int latid,lonid,timeid,u10id,v10id;
250  size_t Nlon,Nlat,Ntime;
251  if ((retval = nc_open(fname,NC_NOWRITE, &ncid)))
252    ERR(retval);
253  if ((retval = nc_inq_dimid(ncid,LONDIM,&dlonid)))
254    ERR(retval);
255  if ((retval = nc_inq_dimid(ncid,LATDIM,&dlatid)))
256    ERR(retval);
257  if ((retval = nc_inq_dimid(ncid,"time",&dtimeid)))
258    ERR(retval);
259  if ((retval = nc_inq_dimlen(ncid,dlonid,&Nlon)))
260    ERR(retval);
261  if ((retval = nc_inq_dimlen(ncid,dlatid,&Nlat)))
262    ERR(retval);
263  if ((retval = nc_inq_dimlen(ncid,dtimeid,&Ntime)))
264    ERR(retval);
265 
266  if ((Nlon != nlon) | (Nlat != nlat)) {
267   fprintf(stderr,
268           "(init_wind) Wrong dimension (nlat=%d,nlon=%d) for the wind_field %s\n",
269           (int)Nlat, (int)Nlon, fname);
270  }
271  if ((retval = nc_inq_varid(ncid,LATDIM,&latid)))
272    ERR(retval);
273  if ((retval = nc_inq_varid(ncid,LONDIM,&lonid)))
274    ERR(retval);
275  if ((retval = nc_inq_varid(ncid,"time",&timeid)))
276    ERR(retval);
277  if ((retval = nc_inq_varid(ncid,UVAR,&u10id)))
278    ERR(retval);
279  if ((retval = nc_inq_varid(ncid,VVAR,&v10id)))
280    ERR(retval);
281 
282  int time[Ntime];
283  float (*u10)[nlat][nlon] = new float[Ntime][nlat][nlon];
284  float (*v10)[nlat][nlon] = new float[Ntime][nlat][nlon];
285
286  //  float u10[Ntime][Nlat][Nlon];
287  //float v10[Ntime][Nlat][Nlon];
288 
289  if ((retval = nc_get_var_float(ncid,u10id,&u10[0][0][0])))
290    ERR(retval);
291 
292  if ((retval = nc_get_var_float(ncid,v10id,&v10[0][0][0])))
293    ERR(retval);
294
295 
296  if ((retval = nc_get_var_int(ncid,timeid,&time[0])))
297    ERR(retval);
298 
299
300  //Load time dimensions
301  if (ptime0 == NULL)
302    {
303      ptime0 = (int *)malloc(sizeof(int));
304      *ptime0 = time[0];
305    }
306
307  //Copy arrays
308  int nt=0;
309  int realtime=*ptime0; //time value from time (in hour)
310  int timeindex=iter0; //time index for uwind and vwind
311  int it;
312  YREAL alpha; //interp coefficient
313  for (it=0;it<(int)Ntime-1;it++) {
314   
315// ionterp whil relatime is between 2 step time
316    while (time[it]<=realtime && time[it+1]>realtime) {
317     
318
319      if (timeindex >= jptfl) {
320        fprintf(stderr,"(init_wind) Warning: time values in %s don't fill the jptfl\n",fname);
321        return;
322      }
323      alpha = ((double)realtime - time[it])/(time[it+1]-time[it]);
324      for (int ilon=0;ilon<nlon;ilon++)
325        for (int ilat=0;ilat<nlat;ilat++) {
326          uwind[ilon][ilat][timeindex] = (1-alpha)*u10[it][ilat][ilon] + alpha*u10[it+1][ilat][ilon];
327          vwind[ilon][ilat][timeindex] = (1-alpha)*v10[it][ilat][ilon] + alpha*v10[it+1][ilat][ilon];
328         
329        }
330      //To avoid cumulative rounding mistake, realtime is calculate with :
331      nt++;
332      realtime = *ptime0 + nt*(double)rdtflo/3600;
333      timeindex ++;
334
335    } //while time[it]<=realtime
336  } //for it
337  if (time[it]==realtime && timeindex < jptfl  ) //last step time
338    {
339      for (int ilon=0;ilon<nlon;ilon++)
340        for (int ilat=0;ilat<nlat;ilat++) {
341          uwind[ilon][ilat][timeindex] = u10[it][ilat][ilon] ;
342          vwind[ilon][ilat][timeindex] = v10[it][ilat][ilon] ;
343         
344        }
345    }
346 
347}
348
349
350
351void read_mask(int argv, char *argc[]) {
352char fname[50];
353if (argv==1) {
354        sprintf(fname,"../obs_float/mask.dat");
355}
356if (argv==2) {
357        sprintf(fname,argc[1]);
358}
359if (argv>2)
360 {
361        fprintf(stderr,"(read_mask) Error %d : bad number of arguments\n",argv);
362        return;
363}
364FILE *fid;
365fid=fopen(fname,"r");
366if (fid==NULL) {
367   fprintf(stderr,"(read_mask) Error : unable to open %s\n",fname);
368   return;
369 }
370 int ilon,ilat,status;
371float val;
372for (ilat=0;ilat<nlat;ilat++)
373        for (ilon=0;ilon<nlon;ilon++) {
374        status=fscanf(fid,"%f",&val);
375        mask[ilon][ilat]=(int)val;
376        YS1_mask(ilon,ilat)=mask[ilon][ilat];
377if (status!=1) {
378        fprintf(stderr,"(read_mask) unable to read mask value (lon=%d,lat=%d)\n",ilon,ilat);
379        fclose(fid);   
380        return;
381}
382}
383fclose(fid);
384
385}
386
387void read_bck (int argv, char *argc[]) {
388/* Read background in a ascii file
389   format
390   Column 1 (%d) : time step (from 1 to jptfl)
391   Column 2 (%d) : lon (Yi) (from 1 to nlon)
392   Column 3 (%d) : lat (Yj) (from 1 to nlat)
393   Column 4 (%f) : u value
394   Column 5 (%f) : v value
395*/
396
397if (argv!=2) {
398  fprintf(stderr,"(read_bck) Error : read_bck with %d (1 nedded)\n",argv-1);
399  fprintf(stderr,"No background read\n");
400  return;
401 }
402 FILE *fid;
403 fid=fopen(argc[1],"r");
404 if (fid==NULL) {
405   fprintf(stderr,"(read_bck) Error : unable to open %s\n",argc[1]);
406   return;
407 }
408 
409 int it,ilon,ilat,count;
410 count=0;
411 YREAL uval,vval;
412 int status=5;
413 while (status==5) {
414   status=fscanf(fid,"%d %d %d %lf %lf",&it,&ilon,&ilat,&uval,&vval);
415   if (status!=5) break;
416   if (ilon>nlon || ilon<=0 || ilat>nlat || ilat<=0 || it<=0 || it>jptfl) {
417     fprintf(stderr,"(read_bck) Error : bad coordinate values in %s (ilon=%d, ilat=%d, it=%d)\n",argc[1],ilon,ilat,it);
418     fclose(fid);
419     return;
420   }
421   ubck[ilon-1][ilat-1][it-1]=uval;
422   vbck[ilon-1][ilat-1][it-1]=vval;
423   count++;
424 }
425 if (count%(nlon*nlat)!=0) 
426   fprintf(stderr,"(read_bck) Warning : only %d data was loaded(nlon*nlat=%d)\n",count,nlon*nlat);
427 fclose(fid);
428 fprintf(stdout,"(read_bck) %d was loaded as background\n",count);
429 
430}
431
432void read_obs_2(int argv, char *argc[]) {
433  /* Read observation in a ascii file compatible with save_output_rfloat
434     format
435     Column 1 (%d) : time step (from 0 to jptfl-1)
436     Column 2 (%d) : idfloat (from 0 to jpnfl-1)
437     Column 3 (%f) : grid point piret (lat)
438     Column 4 (%f) : grid point pjret (lon)
439  */
440  if (argv!=2) {
441    printf("(read_obs) Error : read_obs with %d (1 nedded)\n",argv-1);
442    printf("No obs read\n");
443    return;
444  }
445 FILE *fid;
446  int count=0;
447  fid=fopen(argc[1],"r");
448  if (fid==NULL) {
449    printf("(read_obs) Error : unable to open %s\n",argc[1]);
450    return;
451  }
452 int it,id;
453 YREAL posi,posj;
454  int status=4;
455  while(status==4) {
456    status=fscanf(fid,"%d %d %lf %lf",&it,&id,&posi,&posj);
457    if (status!=4) break;
458    count++;
459    piobs[it][id]=posi;
460    pjobs[it][id]=posj;
461    pmask[it][id]=1;
462  }//while (status==4)
463  fclose(fid);
464  fprintf(stdout,"(read_obs) %d data were loaded as observations\n",count);
465}
466
467void read_obs(int argv, char *argc[]) {
468  /* Read observation in a ascii file
469     format
470     Column 1 (%d) : dimension (1 : lon(Yi), 2 lat(Yj))
471     Column 2 (%d) : time step (from 1 to jptfl)
472     Column 3 (%d) : idfloat
473     Column 4 (%f) : grid point of the floater
474  */
475
476
477  if (argv!=2) {
478    printf("(read_obs) Error : read_obs with %d (1 nedded)\n",argv-1);
479    printf("No obs read\n");
480    return;
481  }
482  FILE *fid;
483  int count=0;
484  fid=fopen(argc[1],"r");
485  if (fid==NULL) {
486    printf("(read_obs) Error : unable to open %s\n",argc[1]);
487    return;
488  }
489  int i,it,id;
490  YREAL pos;
491  int status=4;
492  while(status==4) {
493    status=fscanf(fid,"%d %d %d %lf",&i,&it,&id,&pos);
494    if (status!=4) break;
495    pjobs[it-1][id-1]=pos;
496    status=fscanf(fid,"%d %d %d %lf",&i,&it,&id,&pos);
497    if (status!=4 && i!=2) {
498      printf("(read_obs) Error : unable to read obs in %s\n",argc[1]);
499      return;
500    }
501    count++;
502    piobs[it-1][id-1]=pos;
503    pmask[it-1][id-1]=1;
504  }
505  fclose(fid);
506  fprintf(stdout,"(read_obs) %d data were loaded as observations\n",count);
507}
508
509void update_uv(){
510  int j,k,t;
511  for (t=0;t<jtlag;t++)
512    for (j=0;j<nlon;j++)
513      for (k=0;k<nlat;k++) {
514        YS1_u(j,k,t+1)+=YS1_u_d(j,k);
515        YS1_v(j,k,t+1)+=YS1_v_d(j,k);
516        //printf("u=%f delu=%f \n",YS1_u(j,k),YS1_u_d(j,k));   
517
518}
519}
520
521short is_activ(int jpfl, int it) {
522  /* return 1 if float jpfl is active after it (stricly) */
523
524  //If forward mode, dont test the flaot ater it
525  if (tassim==FORW)
526    return 1;
527 
528  int jptend=it+jtlag+1; 
529  it++;
530 
531  while(it<jptend && pmask[it][jpfl]==0)
532    it++;
533  return(it<jptend && pmask[it][jpfl]==1);
534}
535
536void load_eb(int iti, int itn) {
537  /*Load background from ubck and vbck from time iti to itn*/
538  int j,k,it;
539    for (j=0;j<nlon;j++)
540      for (k=0;k<nlat;k++) {
541        for (it=iti;it<itn;it++)
542          {
543            //+1 to tak uptime into account
544            YS1_u(j,k,it-iti+1)=ubck[j][k][it];
545            YS1_v(j,k,it-iti+1)=vbck[j][k][it];
546            YS1_uwind(j,k,it-iti+1)=uwind[j][k][it];
547            YS1_vwind(j,k,it-iti+1)=vwind[j][k][it];
548           
549          }
550#ifdef K_FILTER
551        YS1_uc_d(j,k)=0;
552        YS1_vc_d(j,k)=0;
553        // Yobs_insert_data (char *nmmod, int sortie, int iaxe, int jaxe, int kaxe,int pdt, YREAL val)
554       
555        Yobs_insert_data("uc_d",0,j,k,0,-1,0);
556        Yobs_insert_data("vc_d",0,j,k,0,-1,0);
557#else
558        Yobs_insert_data("u_d",0,j,k,0,-1,0);
559        Yobs_insert_data("v_d",0,j,k,0,-1,0);
560       
561#endif
562       
563      }
564}
565
566void save_uv(int iti, int itn) {
567  /* save u,v fields in umod,vmod
568   and save piret,pjret*/
569  int j,k,it;
570  int Yifloat=0;
571  for (it=iti;it<itn;it++) {
572    Yifloat=0;
573    for (j=0;j<nlon;j++)
574      for (k=0;k<nlat;k++) {
575        umod[j][k][it]=YS1_u(j,k,it-iti+1);
576        vmod[j][k][it]=YS1_v(j,k,it-iti+1);
577        utotal[j][k][it]=YS1_utot(j,k,it-iti+1);
578        vtotal[j][k][it]=YS1_vtot(j,k,it-iti+1);
579#ifdef UPDATE_BCK
580        ubck[j][k][it]=YS1_u(j,k,it-iti+1);
581        vbck[j][k][it]=YS1_v(j,k,it-iti+1);
582#endif
583      }
584    for (j=0;j<jpnfl;j++) 
585      if (pmask[iti][j]==1 && is_activ(j,iti)>0 && it>iti) {
586        pjret[it][j]=YS1_r_float(Yifloat,it-iti);
587        piret[it][j]=YS2_r_float(Yifloat,it-iti);
588        Yifloat++;
589      }
590  }
591}
592
593void save_output_rfloat (int argc, char *argv[]) {
594  FILE *fid;
595  fid=fopen(argv[1],"w");
596  if (fid==NULL) {
597    printf("\nfailed to open %s",argv[1]);
598    exit(3);
599  }
600  int j,it;
601  for (j=0;j<jpnfl;j++)
602    for (it=0;it<jptfl;it++)
603     
604      fprintf(fid,"%d %d %f %f\n",it,j,piret[it][j],pjret[it][j]);
605       
606                                             
607  fclose(fid);
608
609}
610
611
612void save_output_wind (int argc, char *argv[]) {
613  //Function used only for test
614FILE *fid;
615  fid=fopen(argv[1],"w");
616  if (fid==NULL) {
617    printf("\nfailed to open %s",argv[1]);
618    exit(3);
619  }
620  int j,k,it;
621    for (it=0;it<jptfl;it++)
622      for (j=0;j<nlon;j++)
623        for (k=0;k<nlat;k++) {
624          fprintf(fid,"%d %d %d %f %f\n",it,j,k,uwind[j][k][it],vwind[j][k][it]);
625        }
626                                             
627  fclose(fid);
628
629}
630
631void save_output_uv (int argc, char *argv[]) {
632  FILE *fid;
633  fid=fopen(argv[1],"w");
634  if (fid==NULL) {
635    printf("\nfailed to open %s",argv[1]);
636    exit(3);
637  }
638  int j,k,it;
639    for (it=0;it<jptfl;it++)
640      for (j=0;j<nlon;j++)
641        for (k=0;k<nlat;k++) {
642          fprintf(fid,"%d %d %d %f %f\n",it,j,k,umod[j][k][it],vmod[j][k][it]);
643        }
644                                             
645  fclose(fid);
646
647}
648
649void save_output_uvtot (int argc, char *argv[]) {
650  FILE *fid;
651  fid=fopen(argv[1],"w");
652  if (fid==NULL) {
653    printf("\nfailed to open %s",argv[1]);
654    exit(3);
655  }
656  int j,k,it;
657    for (it=0;it<jptfl;it++)
658      for (j=0;j<nlon;j++)
659        for (k=0;k<nlat;k++) {
660          fprintf(fid,"%d %d %d %f %f\n",it,j,k,utotal[j][k][it],vtotal[j][k][it]);
661        }
662                                             
663  fclose(fid);
664
665}
666
667void save_output_uinter (int argc, char *argv[]) {
668  FILE *fid;
669  fid=fopen(argv[1],"w");
670  if (fid==NULL) {
671    printf("\nfailed to open %s",argv[1]);
672    exit(3);
673  }
674  int j,it;
675  for (j=0;j<jpnfl;j++)
676    for (it=0;it<jptfl;it++)
677     
678      fprintf(fid,"%d %d %f %f\n",it,j,piret[it][j],pjret[it][j]);
679       
680                                             
681  fclose(fid);
682
683}
684
685void div_init() {
686  //Make pseudo-obs of divergence to constrain
687  int i,j,it;
688  for (it=0;it<jtlag;it++)
689    for (j=0;j<nlon;j++)
690      for (i=0;i<nlat;i++) {
691        Yobs_insert_data("cout_div",0,j,i,0,it,0);
692      }
693
694}
695
696void load_init(int it) {
697  /*Load floater position at time it as init.
698  Only floats that have another position after it
699  Warning : nbre de flotteurs maxiums : nfloat */
700  int jfl;
701  int Yifloat=0;
702  int nfobs=0;
703  int ninit=0;
704  int jptend=it+1+jtlag;
705  for (jfl=0;jfl<jpnfl;jfl++) {
706    nfobs=0;
707    ninit=0;
708    if (pmask[it][jfl]==1 && is_activ(jfl,it)>0) {
709      //Chargemente de l'init
710      YS1_r_float(Yifloat,0)=pjobs[it][jfl];
711      YS2_r_float(Yifloat,0)=piobs[it][jfl];
712      ninit++;
713      //Chargement des obs
714      for (int itfl=it+1;itfl<jptend;itfl++) {
715        if (pmask[itfl][jfl]==1) {
716          // Yobs_insert_data (char *nmmod, int sortie, int iaxe, int jaxe, int kaxe,int pdt, YREAL val)
717          Yobs_insert_data("r_cout_d",0,Yifloat,0,0,itfl-it,pjobs[itfl][jfl]);
718          Yobs_insert_data("r_cout_d",1,Yifloat,0,0,itfl-it,piobs[itfl][jfl]);
719        nfobs++;
720        } //if pmask==1
721        } //for itfl
722     
723      Yifloat++; 
724      printf("--Float nr %d : %d initial position / %d obs loaded\n",Yifloat,ninit,nfobs); 
725    } //if pmask==1
726
727    } //for jfl
728 
729 
730  }
731
732void myforward() {
733  int niti=0; //Temps init
734  int nitn=jptfl-1; //Temps fin
735  char sactiv[20]="activ";
736  char sMD[20]="M";
737  char sonly[20]="only";
738  char *liste[3];
739  liste[0]=sactiv;
740  liste[1]=sMD;
741  liste[2]=sonly;
742 
743  //Mode forward
744  tassim=FORW;
745
746  //Effacement des obs potentielles
747  erase_lobs();
748 
749  // Chargement de toutes les init et obs dans les bon pas de temps
750  load_init(niti);
751 
752  //obs pour divergence
753  div_init();
754 
755  //Réglage de l'ébauche
756   load_eb(niti,nitn);
757 
758   //Forward
759   strcpy(liste[1],"M");
760   Yactraj(3, liste);
761   Yset_modeltime(0);
762   before_it(1);
763
764   Yforward(-1, 0);
765
766   //Save outputs
767   save_uv(niti,nitn);
768
769  //On se remet au mode ASSIM par défaut
770  tassim=ASSIM;
771}
772
773void multirun() {
774  int niti=0; //Temps init
775  int nitn=jtlag;//Temps fin
776    while(nitn<jptfl) {
777     
778      //Effacement de toutes les obs précédentes
779      erase_lobs();
780
781      // Chargement de toutes les init et obs dans les bon pas de temps
782      load_init(niti);
783
784 
785      //obs pour divergence
786      div_init();
787
788      // Réglage de l'ébauche
789      load_eb(niti,nitn);
790
791      // Run_inc
792      run_inc(0,NULL);
793
794      // Sauvegarde du résultat dans les tableaux
795      save_uv(niti,nitn);
796
797
798      niti=nitn;
799      nitn=niti+jtlag;
800    } //while nitn<jptfl
801}
802
803void erase_udvd() {
804  int j,k;
805  for (j=0;j<nlon;j++)
806        for (k=0;k<nlat;k++) {
807#ifdef K_FILTER
808          YS1_uc_d(j,k)=0;
809          YS1_vc_d(j,k)=0;
810#else     
811          YS1_u_d(j,k)=0;
812          YS1_v_d(j,k)=0;
813#endif
814}
815}
816//void load_aviso(int argc, char *argv[]){
817       
818//      }
819
820void run_inc(int argv, char *argc[]) {
821//Run the incremental optimisation
822//1 forward of the complete model
823//2 initialize du et dv to zero
824//3 minimize the incremental cost function
825//4 update u and v
826
827/*
828Options
829nb_extiter : number of extern loop (default=6)
830inc_save : save option
831        0  : no save
832        1 (default): save after each extern loop r_float,u and v
833dirsave : directory in case inc_save>0 (default="../obs_float/")
834savename : basename of save fils (default "optim")
835*/
836char filesave[100];
837int i;
838
839char sactiv[20]="activ";
840char sMD[20]="M";
841char sonly[20]="only";
842char *liste[3];
843liste[0]=sactiv;
844liste[1]=sMD;
845liste[2]=sonly;
846/*Check M1QN3 config*/
847 if (YNbItRun<=0)
848                {        printf("runm(2): number of run iteration not seted; use set_nbiter please\n");
849                         //return(0);
850                }
851        if (Y3ddf1<=0)
852                {       printf("runm(2): expected positive fcost decrease missed; use setm_ddf1\n");
853                         //return(0);
854                }
855        if (YioInsertObsCtr<0)
856                   printf("runm(2): warning : oh oh, run with no obs !!! \n");
857               
858        YTypeAdjust = ADJUST_M1QN3; //d'office avec M1QN3
859       
860//Savestate config
861 YioModulot = OFF; //On sauvegarde tous les pas de temps
862    YioWrite = ON;
863    YioRead= OFF;
864   // Yio_savestate(cdes[1], cdes[3], 0, cdes[7]);
865    YioState=0 ;//Save all states
866    YioBin=OFF; //ascii output
867    YioAscii=ON; //ascii output
868    YioAxes=ON; //save axe numbers
869    printf("Run with %d ext_iter with %d obs\n",nb_extiter,YioInsertObsCtr);
870for (i=0;i<nb_extiter;i++) {
871printf("Ext loop %d/%d\n",i+1,nb_extiter);
872     strcpy(liste[1],"M");
873         Yactraj(3, liste);
874         Yset_modeltime(0);
875         before_it(1);
876         //printf("---forward(i=%d)---\n",i);
877         Yforward(-1, 0);
878         
879         if (inc_save==1 && i>0) {
880          //save rfloat
881          YioTime=ON;
882          sprintf(filesave,"%sr_float_%s_%d.dat",dirsave,savename,i);
883          Yio_savestate("r_float","i",0,filesave);
884         }
885             strcpy(liste[1],"D");
886         Yactraj(3,liste);
887         Yset_modeltime(0);
888         
889       
890        Y3run ('0');
891#ifdef K_FILTER
892        //compute values for u_d and v_d
893        strcpy(liste[1],"Tfil");
894        Yactraj(3,liste);
895        Yset_modeltime(0);
896        before_it(1);
897        Yforward(-1, 0);
898        strcpy(liste[1],"T_euler_d");
899        Yactraj(3,liste);
900        Yset_modeltime(0);
901        before_it(1);
902        Yforward(-1, 0);
903#endif
904
905
906    update_uv();
907    erase_udvd();
908
909    if (inc_save==1) {
910   
911   
912    //save u
913    YioTime=ON;
914    sprintf(filesave,"%su%s%d.dat",dirsave,savename,i+1);
915    Yio_savestate("u","ij",0,filesave);
916   
917    //save v
918    sprintf(filesave,"%sv%s%d.dat",dirsave,savename,i+1);
919    Yio_savestate("v","ij",0,filesave);
920   
921    }
922    }
923    /*End of the optimization, we run another time the model to compute the final value of r_float
924    */
925     strcpy(liste[1],"M");
926    Yactraj(3, liste);
927         Yset_modeltime(0);
928         before_it(1);
929         Yforward(-1, 0);
930 if (inc_save==1) {
931          //save rfloat
932          YioTime=ON;
933          sprintf(filesave,"%sr_float_%s_%d.dat",dirsave,savename,nb_extiter);
934          Yio_savestate("r_float","i",0,filesave);
935         }
936}
937
938
Note: See TracBrowser for help on using the repository browser.