source: altifloat/src/floater.h

Last change on this file was 217, checked in by jbrlod, 8 years ago

correcting bug when specifigyer time0 in load_wind
File size: 21.7 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 = (int *)malloc(sizeof(int));
242      *time0 = atoi(argc[3]);
243    }
244  }
245  init_wind(fname,iter0,time0);
246}
247
248void init_wind(char fname[],int iter0,int *ptime0) {
249  int ncid,retval;
250  int dlatid,dlonid,dtimeid;
251  int latid,lonid,timeid,u10id,v10id;
252  size_t Nlon,Nlat,Ntime;
253  if ((retval = nc_open(fname,NC_NOWRITE, &ncid)))
254    ERR(retval);
255  if ((retval = nc_inq_dimid(ncid,LONDIM,&dlonid)))
256    ERR(retval);
257  if ((retval = nc_inq_dimid(ncid,LATDIM,&dlatid)))
258    ERR(retval);
259  if ((retval = nc_inq_dimid(ncid,"time",&dtimeid)))
260    ERR(retval);
261  if ((retval = nc_inq_dimlen(ncid,dlonid,&Nlon)))
262    ERR(retval);
263  if ((retval = nc_inq_dimlen(ncid,dlatid,&Nlat)))
264    ERR(retval);
265  if ((retval = nc_inq_dimlen(ncid,dtimeid,&Ntime)))
266    ERR(retval);
267 
268  if ((Nlon != nlon) | (Nlat != nlat)) {
269   fprintf(stderr,
270           "(init_wind) Wrong dimension (nlat=%d,nlon=%d) for the wind_field %s\n",
271           (int)Nlat, (int)Nlon, fname);
272  }
273  if ((retval = nc_inq_varid(ncid,LATDIM,&latid)))
274    ERR(retval);
275  if ((retval = nc_inq_varid(ncid,LONDIM,&lonid)))
276    ERR(retval);
277  if ((retval = nc_inq_varid(ncid,"time",&timeid)))
278    ERR(retval);
279  if ((retval = nc_inq_varid(ncid,UVAR,&u10id)))
280    ERR(retval);
281  if ((retval = nc_inq_varid(ncid,VVAR,&v10id)))
282    ERR(retval);
283 
284  int time[Ntime];
285  float (*u10)[nlat][nlon] = new float[Ntime][nlat][nlon];
286  float (*v10)[nlat][nlon] = new float[Ntime][nlat][nlon];
287
288  //  float u10[Ntime][Nlat][Nlon];
289  //float v10[Ntime][Nlat][Nlon];
290 
291  if ((retval = nc_get_var_float(ncid,u10id,&u10[0][0][0])))
292    ERR(retval);
293 
294  if ((retval = nc_get_var_float(ncid,v10id,&v10[0][0][0])))
295    ERR(retval);
296
297 
298  if ((retval = nc_get_var_int(ncid,timeid,&time[0])))
299    ERR(retval);
300 
301
302  //Load time dimensions
303  if (ptime0 == NULL)
304    {
305      ptime0 = (int *)malloc(sizeof(int));
306      *ptime0 = time[0];
307    }
308
309  //Copy arrays
310  int nt=0;
311  int realtime=*ptime0; //time value from time (in hour)
312  int timeindex=iter0; //time index for uwind and vwind
313  int it;
314  YREAL alpha; //interp coefficient
315  for (it=0;it<(int)Ntime-1;it++) {
316   
317// ionterp whil relatime is between 2 step time
318    while (time[it]<=realtime && time[it+1]>realtime) {
319     
320
321      if (timeindex >= jptfl) {
322        fprintf(stderr,"(init_wind) Warning: time values in %s don't fill the jptfl\n",fname);
323        return;
324      }
325      alpha = ((double)realtime - time[it])/(time[it+1]-time[it]);
326      for (int ilon=0;ilon<nlon;ilon++)
327        for (int ilat=0;ilat<nlat;ilat++) {
328          uwind[ilon][ilat][timeindex] = (1-alpha)*u10[it][ilat][ilon] + alpha*u10[it+1][ilat][ilon];
329          vwind[ilon][ilat][timeindex] = (1-alpha)*v10[it][ilat][ilon] + alpha*v10[it+1][ilat][ilon];
330         
331        }
332      //To avoid cumulative rounding mistake, realtime is calculate with :
333      nt++;
334      realtime = *ptime0 + nt*(double)rdtflo/3600;
335      timeindex ++;
336
337    } //while time[it]<=realtime
338  } //for it
339  if (time[it]==realtime && timeindex < jptfl  ) //last step time
340    {
341      for (int ilon=0;ilon<nlon;ilon++)
342        for (int ilat=0;ilat<nlat;ilat++) {
343          uwind[ilon][ilat][timeindex] = u10[it][ilat][ilon] ;
344          vwind[ilon][ilat][timeindex] = v10[it][ilat][ilon] ;
345         
346        }
347    }
348 
349}
350
351
352
353void read_mask(int argv, char *argc[]) {
354char fname[50];
355if (argv==1) {
356        sprintf(fname,"../obs_float/mask.dat");
357}
358if (argv==2) {
359        sprintf(fname,argc[1]);
360}
361if (argv>2)
362 {
363        fprintf(stderr,"(read_mask) Error %d : bad number of arguments\n",argv);
364        return;
365}
366FILE *fid;
367fid=fopen(fname,"r");
368if (fid==NULL) {
369   fprintf(stderr,"(read_mask) Error : unable to open %s\n",fname);
370   return;
371 }
372 int ilon,ilat,status;
373float val;
374for (ilat=0;ilat<nlat;ilat++)
375        for (ilon=0;ilon<nlon;ilon++) {
376        status=fscanf(fid,"%f",&val);
377        mask[ilon][ilat]=(int)val;
378        YS1_mask(ilon,ilat)=mask[ilon][ilat];
379if (status!=1) {
380        fprintf(stderr,"(read_mask) unable to read mask value (lon=%d,lat=%d)\n",ilon,ilat);
381        fclose(fid);   
382        return;
383}
384}
385fclose(fid);
386
387}
388
389void read_bck (int argv, char *argc[]) {
390/* Read background in a ascii file
391   format
392   Column 1 (%d) : time step (from 1 to jptfl)
393   Column 2 (%d) : lon (Yi) (from 1 to nlon)
394   Column 3 (%d) : lat (Yj) (from 1 to nlat)
395   Column 4 (%f) : u value
396   Column 5 (%f) : v value
397*/
398
399if (argv!=2) {
400  fprintf(stderr,"(read_bck) Error : read_bck with %d (1 nedded)\n",argv-1);
401  fprintf(stderr,"No background read\n");
402  return;
403 }
404 FILE *fid;
405 fid=fopen(argc[1],"r");
406 if (fid==NULL) {
407   fprintf(stderr,"(read_bck) Error : unable to open %s\n",argc[1]);
408   return;
409 }
410 
411 int it,ilon,ilat,count;
412 count=0;
413 YREAL uval,vval;
414 int status=5;
415 while (status==5) {
416   status=fscanf(fid,"%d %d %d %lf %lf",&it,&ilon,&ilat,&uval,&vval);
417   if (status!=5) break;
418   if (ilon>nlon || ilon<=0 || ilat>nlat || ilat<=0 || it<=0 || it>jptfl) {
419     fprintf(stderr,"(read_bck) Error : bad coordinate values in %s (ilon=%d, ilat=%d, it=%d)\n",argc[1],ilon,ilat,it);
420     fclose(fid);
421     return;
422   }
423   ubck[ilon-1][ilat-1][it-1]=uval;
424   vbck[ilon-1][ilat-1][it-1]=vval;
425   count++;
426 }
427 if (count%(nlon*nlat)!=0) 
428   fprintf(stderr,"(read_bck) Warning : only %d data was loaded(nlon*nlat=%d)\n",count,nlon*nlat);
429 fclose(fid);
430 fprintf(stdout,"(read_bck) %d was loaded as background\n",count);
431 
432}
433
434void read_obs_2(int argv, char *argc[]) {
435  /* Read observation in a ascii file compatible with save_output_rfloat
436     format
437     Column 1 (%d) : time step (from 0 to jptfl-1)
438     Column 2 (%d) : idfloat (from 0 to jpnfl-1)
439     Column 3 (%f) : grid point piret (lat)
440     Column 4 (%f) : grid point pjret (lon)
441  */
442  if (argv!=2) {
443    printf("(read_obs) Error : read_obs with %d (1 nedded)\n",argv-1);
444    printf("No obs read\n");
445    return;
446  }
447 FILE *fid;
448  int count=0;
449  fid=fopen(argc[1],"r");
450  if (fid==NULL) {
451    printf("(read_obs) Error : unable to open %s\n",argc[1]);
452    return;
453  }
454 int it,id;
455 YREAL posi,posj;
456  int status=4;
457  while(status==4) {
458    status=fscanf(fid,"%d %d %lf %lf",&it,&id,&posi,&posj);
459    if (status!=4) break;
460    count++;
461    piobs[it][id]=posi;
462    pjobs[it][id]=posj;
463    pmask[it][id]=1;
464  }//while (status==4)
465  fclose(fid);
466  fprintf(stdout,"(read_obs) %d data were loaded as observations\n",count);
467}
468
469void read_obs(int argv, char *argc[]) {
470  /* Read observation in a ascii file
471     format
472     Column 1 (%d) : dimension (1 : lon(Yi), 2 lat(Yj))
473     Column 2 (%d) : time step (from 1 to jptfl)
474     Column 3 (%d) : idfloat
475     Column 4 (%f) : grid point of the floater
476  */
477
478
479  if (argv!=2) {
480    printf("(read_obs) Error : read_obs with %d (1 nedded)\n",argv-1);
481    printf("No obs read\n");
482    return;
483  }
484  FILE *fid;
485  int count=0;
486  fid=fopen(argc[1],"r");
487  if (fid==NULL) {
488    printf("(read_obs) Error : unable to open %s\n",argc[1]);
489    return;
490  }
491  int i,it,id;
492  YREAL pos;
493  int status=4;
494  while(status==4) {
495    status=fscanf(fid,"%d %d %d %lf",&i,&it,&id,&pos);
496    if (status!=4) break;
497    pjobs[it-1][id-1]=pos;
498    status=fscanf(fid,"%d %d %d %lf",&i,&it,&id,&pos);
499    if (status!=4 && i!=2) {
500      printf("(read_obs) Error : unable to read obs in %s\n",argc[1]);
501      return;
502    }
503    count++;
504    piobs[it-1][id-1]=pos;
505    pmask[it-1][id-1]=1;
506  }
507  fclose(fid);
508  fprintf(stdout,"(read_obs) %d data were loaded as observations\n",count);
509}
510
511void update_uv(){
512  int j,k,t;
513  for (t=0;t<jtlag;t++)
514    for (j=0;j<nlon;j++)
515      for (k=0;k<nlat;k++) {
516        YS1_u(j,k,t+1)+=YS1_u_d(j,k);
517        YS1_v(j,k,t+1)+=YS1_v_d(j,k);
518        //printf("u=%f delu=%f \n",YS1_u(j,k),YS1_u_d(j,k));   
519
520}
521}
522
523short is_activ(int jpfl, int it) {
524  /* return 1 if float jpfl is active after it (stricly) */
525
526  //If forward mode, dont test the flaot ater it
527  if (tassim==FORW)
528    return 1;
529 
530  int jptend=it+jtlag+1; 
531  it++;
532 
533  while(it<jptend && pmask[it][jpfl]==0)
534    it++;
535  return(it<jptend && pmask[it][jpfl]==1);
536}
537
538void load_eb(int iti, int itn) {
539  /*Load background from ubck and vbck from time iti to itn*/
540  int j,k,it;
541    for (j=0;j<nlon;j++)
542      for (k=0;k<nlat;k++) {
543        for (it=iti;it<itn;it++)
544          {
545            //+1 to tak uptime into account
546            YS1_u(j,k,it-iti+1)=ubck[j][k][it];
547            YS1_v(j,k,it-iti+1)=vbck[j][k][it];
548            YS1_uwind(j,k,it-iti+1)=uwind[j][k][it];
549            YS1_vwind(j,k,it-iti+1)=vwind[j][k][it];
550           
551          }
552#ifdef K_FILTER
553        YS1_uc_d(j,k)=0;
554        YS1_vc_d(j,k)=0;
555        // Yobs_insert_data (char *nmmod, int sortie, int iaxe, int jaxe, int kaxe,int pdt, YREAL val)
556       
557        Yobs_insert_data("uc_d",0,j,k,0,-1,0);
558        Yobs_insert_data("vc_d",0,j,k,0,-1,0);
559#else
560        Yobs_insert_data("u_d",0,j,k,0,-1,0);
561        Yobs_insert_data("v_d",0,j,k,0,-1,0);
562       
563#endif
564       
565      }
566}
567
568void save_uv(int iti, int itn) {
569  /* save u,v fields in umod,vmod
570   and save piret,pjret*/
571  int j,k,it;
572  int Yifloat=0;
573  for (it=iti;it<itn;it++) {
574    Yifloat=0;
575    for (j=0;j<nlon;j++)
576      for (k=0;k<nlat;k++) {
577        umod[j][k][it]=YS1_u(j,k,it-iti+1);
578        vmod[j][k][it]=YS1_v(j,k,it-iti+1);
579        utotal[j][k][it]=YS1_utot(j,k,it-iti+1);
580        vtotal[j][k][it]=YS1_vtot(j,k,it-iti+1);
581#ifdef UPDATE_BCK
582        ubck[j][k][it]=YS1_u(j,k,it-iti+1);
583        vbck[j][k][it]=YS1_v(j,k,it-iti+1);
584#endif
585      }
586    for (j=0;j<jpnfl;j++) 
587      if (pmask[iti][j]==1 && is_activ(j,iti)>0 && it>iti) {
588        pjret[it][j]=YS1_r_float(Yifloat,it-iti);
589        piret[it][j]=YS2_r_float(Yifloat,it-iti);
590        Yifloat++;
591      }
592  }
593}
594
595void save_output_rfloat (int argc, char *argv[]) {
596  FILE *fid;
597  fid=fopen(argv[1],"w");
598  if (fid==NULL) {
599    printf("\nfailed to open %s",argv[1]);
600    exit(3);
601  }
602  int j,it;
603  for (j=0;j<jpnfl;j++)
604    for (it=0;it<jptfl;it++)
605     
606      fprintf(fid,"%d %d %f %f\n",it,j,piret[it][j],pjret[it][j]);
607       
608                                             
609  fclose(fid);
610
611}
612
613
614void save_output_wind (int argc, char *argv[]) {
615  //Function used only for test
616FILE *fid;
617  fid=fopen(argv[1],"w");
618  if (fid==NULL) {
619    printf("\nfailed to open %s",argv[1]);
620    exit(3);
621  }
622  int j,k,it;
623    for (it=0;it<jptfl;it++)
624      for (j=0;j<nlon;j++)
625        for (k=0;k<nlat;k++) {
626          fprintf(fid,"%d %d %d %f %f\n",it,j,k,uwind[j][k][it],vwind[j][k][it]);
627        }
628                                             
629  fclose(fid);
630
631}
632
633void save_output_uv (int argc, char *argv[]) {
634  FILE *fid;
635  fid=fopen(argv[1],"w");
636  if (fid==NULL) {
637    printf("\nfailed to open %s",argv[1]);
638    exit(3);
639  }
640  int j,k,it;
641    for (it=0;it<jptfl;it++)
642      for (j=0;j<nlon;j++)
643        for (k=0;k<nlat;k++) {
644          fprintf(fid,"%d %d %d %f %f\n",it,j,k,umod[j][k][it],vmod[j][k][it]);
645        }
646                                             
647  fclose(fid);
648
649}
650
651void save_output_uvtot (int argc, char *argv[]) {
652  FILE *fid;
653  fid=fopen(argv[1],"w");
654  if (fid==NULL) {
655    printf("\nfailed to open %s",argv[1]);
656    exit(3);
657  }
658  int j,k,it;
659    for (it=0;it<jptfl;it++)
660      for (j=0;j<nlon;j++)
661        for (k=0;k<nlat;k++) {
662          fprintf(fid,"%d %d %d %f %f\n",it,j,k,utotal[j][k][it],vtotal[j][k][it]);
663        }
664                                             
665  fclose(fid);
666
667}
668
669void save_output_uinter (int argc, char *argv[]) {
670  FILE *fid;
671  fid=fopen(argv[1],"w");
672  if (fid==NULL) {
673    printf("\nfailed to open %s",argv[1]);
674    exit(3);
675  }
676  int j,it;
677  for (j=0;j<jpnfl;j++)
678    for (it=0;it<jptfl;it++)
679     
680      fprintf(fid,"%d %d %f %f\n",it,j,piret[it][j],pjret[it][j]);
681       
682                                             
683  fclose(fid);
684
685}
686
687void div_init() {
688  //Make pseudo-obs of divergence to constrain
689  int i,j,it;
690  for (it=0;it<jtlag;it++)
691    for (j=0;j<nlon;j++)
692      for (i=0;i<nlat;i++) {
693        Yobs_insert_data("cout_div",0,j,i,0,it,0);
694      }
695
696}
697
698void load_init(int it) {
699  /*Load floater position at time it as init.
700  Only floats that have another position after it
701  Warning : nbre de flotteurs maxiums : nfloat */
702  int jfl;
703  int Yifloat=0;
704  int nfobs=0;
705  int ninit=0;
706  int jptend=it+1+jtlag;
707  for (jfl=0;jfl<jpnfl;jfl++) {
708    nfobs=0;
709    ninit=0;
710    if (pmask[it][jfl]==1 && is_activ(jfl,it)>0) {
711      //Chargemente de l'init
712      YS1_r_float(Yifloat,0)=pjobs[it][jfl];
713      YS2_r_float(Yifloat,0)=piobs[it][jfl];
714      ninit++;
715      //Chargement des obs
716      for (int itfl=it+1;itfl<jptend;itfl++) {
717        if (pmask[itfl][jfl]==1) {
718          // Yobs_insert_data (char *nmmod, int sortie, int iaxe, int jaxe, int kaxe,int pdt, YREAL val)
719          Yobs_insert_data("r_cout_d",0,Yifloat,0,0,itfl-it,pjobs[itfl][jfl]);
720          Yobs_insert_data("r_cout_d",1,Yifloat,0,0,itfl-it,piobs[itfl][jfl]);
721        nfobs++;
722        } //if pmask==1
723        } //for itfl
724     
725      Yifloat++; 
726      printf("--Float nr %d : %d initial position / %d obs loaded\n",Yifloat,ninit,nfobs); 
727    } //if pmask==1
728
729    } //for jfl
730 
731 
732  }
733
734void myforward() {
735  int niti=0; //Temps init
736  int nitn=jptfl-1; //Temps fin
737  char sactiv[20]="activ";
738  char sMD[20]="M";
739  char sonly[20]="only";
740  char *liste[3];
741  liste[0]=sactiv;
742  liste[1]=sMD;
743  liste[2]=sonly;
744 
745  //Mode forward
746  tassim=FORW;
747
748  //Effacement des obs potentielles
749  erase_lobs();
750 
751  // Chargement de toutes les init et obs dans les bon pas de temps
752  load_init(niti);
753 
754  //obs pour divergence
755  div_init();
756 
757  //Réglage de l'ébauche
758   load_eb(niti,nitn);
759 
760   //Forward
761   strcpy(liste[1],"M");
762   Yactraj(3, liste);
763   Yset_modeltime(0);
764   before_it(1);
765
766   Yforward(-1, 0);
767
768   //Save outputs
769   save_uv(niti,nitn);
770
771  //On se remet au mode ASSIM par défaut
772  tassim=ASSIM;
773}
774
775void multirun() {
776  int niti=0; //Temps init
777  int nitn=jtlag;//Temps fin
778    while(nitn<jptfl) {
779     
780      //Effacement de toutes les obs précédentes
781      erase_lobs();
782
783      // Chargement de toutes les init et obs dans les bon pas de temps
784      load_init(niti);
785
786 
787      //obs pour divergence
788      div_init();
789
790      // Réglage de l'ébauche
791      load_eb(niti,nitn);
792
793      // Run_inc
794      run_inc(0,NULL);
795
796      // Sauvegarde du résultat dans les tableaux
797      save_uv(niti,nitn);
798
799
800      niti=nitn;
801      nitn=niti+jtlag;
802    } //while nitn<jptfl
803}
804
805void erase_udvd() {
806  int j,k;
807  for (j=0;j<nlon;j++)
808        for (k=0;k<nlat;k++) {
809#ifdef K_FILTER
810          YS1_uc_d(j,k)=0;
811          YS1_vc_d(j,k)=0;
812#else     
813          YS1_u_d(j,k)=0;
814          YS1_v_d(j,k)=0;
815#endif
816}
817}
818//void load_aviso(int argc, char *argv[]){
819       
820//      }
821
822void run_inc(int argv, char *argc[]) {
823//Run the incremental optimisation
824//1 forward of the complete model
825//2 initialize du et dv to zero
826//3 minimize the incremental cost function
827//4 update u and v
828
829/*
830Options
831nb_extiter : number of extern loop (default=6)
832inc_save : save option
833        0  : no save
834        1 (default): save after each extern loop r_float,u and v
835dirsave : directory in case inc_save>0 (default="../obs_float/")
836savename : basename of save fils (default "optim")
837*/
838char filesave[100];
839int i;
840
841char sactiv[20]="activ";
842char sMD[20]="M";
843char sonly[20]="only";
844char *liste[3];
845liste[0]=sactiv;
846liste[1]=sMD;
847liste[2]=sonly;
848/*Check M1QN3 config*/
849 if (YNbItRun<=0)
850                {        printf("runm(2): number of run iteration not seted; use set_nbiter please\n");
851                         //return(0);
852                }
853        if (Y3ddf1<=0)
854                {       printf("runm(2): expected positive fcost decrease missed; use setm_ddf1\n");
855                         //return(0);
856                }
857        if (YioInsertObsCtr<0)
858                   printf("runm(2): warning : oh oh, run with no obs !!! \n");
859               
860        YTypeAdjust = ADJUST_M1QN3; //d'office avec M1QN3
861       
862//Savestate config
863 YioModulot = OFF; //On sauvegarde tous les pas de temps
864    YioWrite = ON;
865    YioRead= OFF;
866   // Yio_savestate(cdes[1], cdes[3], 0, cdes[7]);
867    YioState=0 ;//Save all states
868    YioBin=OFF; //ascii output
869    YioAscii=ON; //ascii output
870    YioAxes=ON; //save axe numbers
871    printf("Run with %d ext_iter with %d obs\n",nb_extiter,YioInsertObsCtr);
872for (i=0;i<nb_extiter;i++) {
873printf("Ext loop %d/%d\n",i+1,nb_extiter);
874     strcpy(liste[1],"M");
875         Yactraj(3, liste);
876         Yset_modeltime(0);
877         before_it(1);
878         //printf("---forward(i=%d)---\n",i);
879         Yforward(-1, 0);
880         
881         if (inc_save==1 && i>0) {
882          //save rfloat
883          YioTime=ON;
884          sprintf(filesave,"%sr_float_%s_%d.dat",dirsave,savename,i);
885          Yio_savestate("r_float","i",0,filesave);
886         }
887             strcpy(liste[1],"D");
888         Yactraj(3,liste);
889         Yset_modeltime(0);
890         
891       
892        Y3run ('0');
893#ifdef K_FILTER
894        //compute values for u_d and v_d
895        strcpy(liste[1],"Tfil");
896        Yactraj(3,liste);
897        Yset_modeltime(0);
898        before_it(1);
899        Yforward(-1, 0);
900        strcpy(liste[1],"T_euler_d");
901        Yactraj(3,liste);
902        Yset_modeltime(0);
903        before_it(1);
904        Yforward(-1, 0);
905#endif
906
907
908    update_uv();
909    erase_udvd();
910
911    if (inc_save==1) {
912   
913   
914    //save u
915    YioTime=ON;
916    sprintf(filesave,"%su%s%d.dat",dirsave,savename,i+1);
917    Yio_savestate("u","ij",0,filesave);
918   
919    //save v
920    sprintf(filesave,"%sv%s%d.dat",dirsave,savename,i+1);
921    Yio_savestate("v","ij",0,filesave);
922   
923    }
924    }
925    /*End of the optimization, we run another time the model to compute the final value of r_float
926    */
927     strcpy(liste[1],"M");
928    Yactraj(3, liste);
929         Yset_modeltime(0);
930         before_it(1);
931         Yforward(-1, 0);
932 if (inc_save==1) {
933          //save rfloat
934          YioTime=ON;
935          sprintf(filesave,"%sr_float_%s_%d.dat",dirsave,savename,nb_extiter);
936          Yio_savestate("r_float","i",0,filesave);
937         }
938}
939
940
Note: See TracBrowser for help on using the repository browser.