Index: /branches/branche-mb/scripts/sinobad.i
===================================================================
--- /branches/branche-mb/scripts/sinobad.i	(revision 67)
+++ /branches/branche-mb/scripts/sinobad.i	(revision 68)
@@ -18,5 +18,10 @@
 xwriteout 0 true ../data_out/exp_T.nc
 #setstate nu 1
-#init_nu
+init_nu
+xset_fdt_fco th 100. 
+xset_fdt_fco tz 100.
+xset_fdt_fco kh 5.e+02
+xset_fdt_fco kz 1.2e-05
+
 print_time OFF
 set_modeltime 0
@@ -36,5 +41,5 @@
 savestate ub  1   ijk   50    A       1       ../data_out/ub_obs_48.dat   
 savestate vb  1   ijk   50    A       1       ../data_out/vb_obs_48.dat   
-#GOTO TESTAD 
+GOTO TESTAD 
 
 loadobs sshb  1   ij   50    A       1       ../data_out/sshb_obs_48.dat D   
Index: /branches/branche-mb/scripts/sinobad_2D.i
===================================================================
--- /branches/branche-mb/scripts/sinobad_2D.i	(revision 67)
+++ /branches/branche-mb/scripts/sinobad_2D.i	(revision 68)
@@ -2,4 +2,5 @@
 # 		fichier de D'INSTRUCTIONS de rÃ©fÃ©rence
 #Julien Brajard
+# modified by mohamed Berrada (filtre August 2011)
 #==============================================================================
 #script faisant l'assimilation de rÃ©fÃ©rence
@@ -14,4 +15,10 @@
 load_eb ta_c   // pour tenir en compte ta_c sinon sera 0 cause filtre
 true_target_in_tab ta_c
+
+xset_fdt_fco th 1000. 
+xset_fdt_fco tz 10.
+xset_fdt_fco kh 5.e+05
+xset_fdt_fco kz 1.2e-03
+
 xinitnc ../data_out/exp_T_2D.nc
 
@@ -40,5 +47,5 @@
 setstate nu 0
 outoebx nu 1 0
-set_bcoef nu 1.e-05
+set_bcoef nu 1.e-8
 set_modeltime 0
 print_time OFF
@@ -69,10 +76,18 @@
 goto nomulti
 load_eb ta_c  // charger l'ebauche '
-set_bcoef nu 1.e-06
+set_bcoef nu 5.e-05
 runm
 
 load_eb ta_c  // charger l'ebauche '
-set_bcoef nu 1.e-07
+set_bcoef nu 1.e-04
 runm
+
+##SAUVEGARDE
+xwriteout 0 fin ../data_out/exp_T_2D.nc
+xwriteout 50 forwfin_48 ../data_out/exp_T_2D.nc
+
+xrst_save ../data_in/file_rest/GYRE_00000400_restart_CONTROL_T_2D.nc
+
+
 nomulti
 
Index: /branches/branche-mb/src/cost.dat
===================================================================
--- /branches/branche-mb/src/cost.dat	(revision 67)
+++ /branches/branche-mb/src/cost.dat	(revision 68)
@@ -2,2 +2,40 @@
  0.0000000000000000e+00
  0.0000000000000000e+00
+ 2.1474823164582097e+02
+ 1.2274234090104164e+01
+ 2.3881198434552111e+01
+ 8.8771059162513382e+00
+ 5.5557975737505458e+00
+ 3.4803779742760317e+00
+ 1.7128978455386892e+00
+ 1.1973455752757605e+00
+ 9.8335187344394392e-01
+ 7.6117334075192200e-01
+ 5.9846985043808543e-01
+ 6.7366624706570466e-01
+ 5.5874401856066103e-01
+ 7.3595639204467600e-01
+ 5.0953722678532765e-01
+ 1.1965283691171660e+00
+ 6.5986339544255534e-01
+ 6.6013045954860117e-01
+ 6.5392286280927792e-01
+ 5.6394373996837344e-01
+ 5.7610237016601717e-01
+ 5.8032116526810040e-01
+ 5.8309360765246787e-01
+ 5.8519910883568982e-01
+ 5.8664512715255190e-01
+ 5.8761359807265934e-01
+ 5.8831541023256795e-01
+ 5.8880781610867361e-01
+ 5.8915693207370023e-01
+ 5.0893021741293576e-01
+ 5.0885469878063883e-01
+ 5.8920763747258753e-01
+ 5.0875296736771736e-01
+ 5.0874187732734499e-01
+ 5.8921509431560293e-01
+ 5.8922030000369108e-01
+ 5.8922395197080846e-01
+ 5.8922651344473387e-01
Index: /branches/branche-mb/src/dta_lfexp.h
===================================================================
--- /branches/branche-mb/src/dta_lfexp.h	(revision 67)
+++ /branches/branche-mb/src/dta_lfexp.h	(revision 68)
@@ -37,18 +37,18 @@
 	tim1= 0.;
 #ifdef K_NEU
-	tim1= - bim1 * x2;
+	tim1= bim1 * x2;
 #endif
       }
       else{
-	tim1= -bim1 * x3;
+	tim1= bim1 * x3;
       }
       if (Yi==NX-2){
 	tip1= 0.;
 #ifdef K_NEU
-	tip1= - bi * x2;
+	tip1= bi * x2;
 #endif
       }
       else{
-	tip1= -bi * x4;
+	tip1= bi * x4;
       }
       //cdt au bord j
@@ -56,21 +56,21 @@
 	tjm1= 0.;
 #ifdef K_NEU
-	tjm1= - bjm1 * x2;
+	tjm1= bjm1 * x2;
 #endif
       }
       else{
-	tjm1= -bjm1 * x5;
+	tjm1= bjm1 * x5;
       }
       if (Yj==NY-2){
 	tjp1= 0.;
 #ifdef K_NEU
-	tjp1= - bj * x2;
+	tjp1= bj * x2;
 #endif
       }
       else{
-	tjp1= -bj * x6;
+	tjp1= bj * x6;
       }
 
-      tij= (1./zbtr+bi+bim1+bj+bjm1) * x2;
+      tij= (1./zbtr-bi-bim1-bj-bjm1) * x2;
       YS1 = zbtr*(tip1 + tim1 + tjp1 + tjm1 + tij);
     }
@@ -105,18 +105,18 @@
 	dtim1_x2= 0.;
 #ifdef K_NEU
-	dtim1_x2= - bim1;
+	dtim1_x2= bim1;
 #endif
       }
       else{
-	dtim1_x3= -bim1;
+	dtim1_x3= bim1;
       }
       if (Yi==NX-2){
 	dtip1_x2= 0.;
 #ifdef K_NEU
-	dtip1_x2= - bi;
+	dtip1_x2= bi;
 #endif
       }
       else{
-	dtip1_x4= -bi;
+	dtip1_x4= bi;
       }
       //cdt au bord j
@@ -124,20 +124,20 @@
 	dtjm1_x2= 0.;
 #ifdef K_NEU
-	dtjm1_x2= - bjm1;
+	dtjm1_x2= bjm1;
 #endif
       }
       else{
-	dtjm1_x5= -bjm1;
+	dtjm1_x5= bjm1;
       }
       if (Yj==NY-2){
 	dtjp1_x2= 0.;
 #ifdef K_NEU
-	dtjp1_x2= - bj;
+	dtjp1_x2= bj;
 #endif
       }
       else{
-	dtjp1_x6= -bj;
+	dtjp1_x6= bj;
       }
-      dtij_x2= (1./zbtr+bi+bim1+bj+bjm1);
+      dtij_x2= (1./zbtr-bi-bim1-bj-bjm1);
 
       YJ1I1 = 0.;
Index: /branches/branche-mb/src/dta_zfimp.h
===================================================================
--- /branches/branche-mb/src/dta_zfimp.h	(revision 68)
+++ /branches/branche-mb/src/dta_zfimp.h	(revision 68)
@@ -0,0 +1,62 @@
+
+/***************************************************************************
+                 module classe dta_zfimp.h  -  description
+***************************************************************************/
+// Mohamed Berrada
+// locean-ipsl.upmc, Paris, August 1, 2011
+//===========================================================================
+//                          methode forward
+forward (YREAL x1,YREAL x2,YREAL x3)
+{ /*	1	from	nu       		1  i 	j    k    
+	2	from	dta_zfimp    		1  i 	j    k+1   t
+	3	from	dta_zfimp_t    		1  i 	j    k     t */
+
+  if(Yt==1){
+    double wm12=1./sqrt(fse3t(Yi,Yj,Yk)); //1/sqrt(dz)
+    YS1=wm12*x1;
+    if(Yi==0 || Yi==NX-1 || Yj==0 || Yj==NY-1 || Yk==NZ-1){
+      YS1=0.;
+    }
+  }
+  else{
+    if(Yi==0 || Yi==NX-1 || Yj==0 || Yj==NY-1 || Yk==NZ-1){
+      YS1=0.;
+    }
+    else{
+      if(Yk==NZ-2)
+	YS1 = x3;
+      else
+	YS1 = x3 - tab_luzimp( Yi, Yj, Yk+1 , 0 ) * x2;
+    }
+  }
+  //
+}
+      
+//===========================================================================
+//                         methode  backward
+
+backward (YREAL x1,YREAL x2,YREAL x3)
+{
+  YJ1I1=0.;     YJ1I2=0.;     YJ1I3=0.; 
+  if(Yt==1){
+    double wm12=1./sqrt(fse3t(Yi,Yj,Yk)); //1/sqrt(dz)
+    YJ1I1=wm12;
+    if(Yi==0 || Yi==NX-1 || Yj==0 || Yj==NY-1 || Yk==NZ-1){
+      YJ1I1=0.;
+    }
+  }
+  else{
+    if(Yi==0 || Yi==NX-1 || Yj==0 || Yj==NY-1 || Yk==NZ-1){
+      YJ1I1=0.;
+    }
+    else{
+      if(Yk==NZ-2)
+	YJ1I3 = 1.;
+      else{
+	YJ1I2 = -tab_luzimp( Yi, Yj, Yk+1 , 0 );
+	YJ1I3 = 1.;
+      } 
+    }
+  }
+  //
+}
Index: /branches/branche-mb/src/dta_zfimp_t.h
===================================================================
--- /branches/branche-mb/src/dta_zfimp_t.h	(revision 68)
+++ /branches/branche-mb/src/dta_zfimp_t.h	(revision 68)
@@ -0,0 +1,54 @@
+
+/***************************************************************************
+                 module classe dta_zfimp_t.h  -  description
+***************************************************************************/
+// Mohamed Berrada
+// locean-ipsl.upmc, Paris, August 1, 2011
+//===========================================================================
+//                          methode forward
+forward (YREAL x1,YREAL x2)
+{ /*	1	from	dta_zfimp_t    		1  i 	j    k-1   t  
+	2	from	dta_zfimp      		1  i 	j    k     t-1*/
+
+  if(Yt==1){
+    YS1=0;
+  }
+  else{
+    if(Yi==0 || Yi==NX-1 || Yj==0 || Yj==NY-1 || Yk==NZ-1){
+      YS1=0.;
+    }
+    else{
+      double a0 =avt_fil*pdtz_fil/fse3w(Yi,Yj,Yk); 
+      if(Yk==0)
+	YS1 = x2/tab_luzimp( Yi, Yj, Yk , 1 ); 
+      else
+	YS1 = (x2 - a0 /fse3t(Yi,Yj,Yk)* x1)/tab_luzimp( Yi, Yj, Yk , 1 );
+    }
+  }
+  //
+}
+      
+//===========================================================================
+//                         methode  backward
+
+backward (YREAL x1,YREAL x2)
+{
+  YJ1I1=0.;     YJ1I2=0.; 
+  if(Yt==1){
+    YJ1I1=0.;     YJ1I2=0.; 
+  }
+  else{
+    if(Yi==0 || Yi==NX-1 || Yj==0 || Yj==NY-1 || Yk==NZ-1){
+      YJ1I1=0.;     YJ1I2=0.; 
+    }
+    else{
+      double a0 =avt_fil*pdtz_fil/fse3w(Yi,Yj,Yk); 
+      if(Yk==0)
+	YJ1I2 = 1./tab_luzimp( Yi, Yj, Yk , 1 ); 
+      else{
+	YJ1I1= -a0/fse3t(Yi,Yj,Yk)/tab_luzimp( Yi, Yj, Yk , 1 );
+	YJ1I2=1./tab_luzimp( Yi, Yj, Yk , 1 ); 
+      }
+    }
+  }
+}
Index: /branches/branche-mb/src/sinobad.d
===================================================================
--- /branches/branche-mb/src/sinobad.d	(revision 67)
+++ /branches/branche-mb/src/sinobad.d	(revision 68)
@@ -1,7 +1,7 @@
-defval K_NEU
+//defval K_NEU
 #define FZIMP
 #define FILTER
-//#define FILZNOL
-#define FILLNOZ
+#define FILZNOL
+//#define FILLNOZ
 //#define SOLSORYAO
 //#define OPTIMORDER
@@ -22,7 +22,7 @@
 defval TU       2 | start time step
 #ifdef FILTER
-defval OFTL      10
-defval OFTZ      0
-defval OFT       10
+defval OFTL      0
+defval OFTZ      4
+defval OFT       4
 #endif
 #ifndef FILTER
@@ -329,4 +329,5 @@
 //#---------->dta_zfil
 #ifndef FILLNOZ
+
 #ifndef FILZNOL
 #ifndef FZIMP
@@ -367,5 +368,7 @@
 ctin	 ta_c   	1	from	dta_zfexp      		1  i 	j    k   t
 #endif
-#endif
+
+#endif
+
 #ifdef FILLNOZ
 ctin	 ta_c   	1	from	dta_lfexp      		1  i 	j    k   t
@@ -1209,5 +1212,6 @@
 #ifdef FILTER
 insert_fct  arg load_eb
-insert_fct  init_nu
+insert_fct  arg init_nu
+insert_fct  arg xset_fdt_fco
 #endif
 insert_fct arg xwriteout
Index: /branches/branche-mb/src/sinobad.h
===================================================================
--- /branches/branche-mb/src/sinobad.h	(revision 67)
+++ /branches/branche-mb/src/sinobad.h	(revision 68)
@@ -87,9 +87,10 @@
 double pdth_fil=100.;//7200.;
 double pdtz_fil=100.;//100.;//7200.;
-double fsaht_fil= 5.e+2;//1000.;
+double fsaht_fil= 5.e+02;//1000.;
 double avt_fil=1.2e-05;//1.2e-5;
 //double lscale= 212000.; // length-scale
 double norm_fil; //coef de normalisation
 void normfil();// calcul du coef de normalisation norm_fil
+void xset_fdt_fco();// set fdt fcoef
 #endif
 
@@ -259,10 +260,11 @@
 //_____________________________________________________________________________
 #ifdef K_FILTER
+//_____________________________________________________________________________
 void normfil(){
   double norm_zfil=1.;
   double norm_lfil=1.;
-#if defined (K_FILZNOL) /* && !defined (K_FILLNOZ)*/ 
-  norm_zfil= 2.*sqrt(PI*avt_fil*pdtz_fil*2.*OFTZ);
-#ifdef K_FZIMP
+  norm_lfil= 2.*sqrt(PI*fsaht_fil*pdth_fil*2.*(OFTL-1.));// schema explicite
+  norm_zfil= 2.*sqrt(PI*avt_fil*pdtz_fil*2.*(OFTZ-1.));// schema explicite 
+#ifdef K_FZIMP  // schema implicite 
   double fact=1.; 
   int M=2*(OFTZ-1);
@@ -271,21 +273,67 @@
   norm_zfil=2.*fact*sqrt(avt_fil*pdtz_fil);
 #endif
-#elif  defined (K_FILLNOZ) /*&& !defined (K_FILZNOL)*/  
-  norm_lfil= 2.*sqrt(PI*fsaht_fil*pdth_fil*2.*OFTL);
-#endif
+
+#if defined (K_FILZNOL)
+  norm_lfil=1.;
+#endif
+#if  defined (K_FILLNOZ)
+  norm_zfil=1.;
+#endif
+
   norm_fil= norm_zfil*norm_lfil;
   printf("norm_fil = %23.16e\n",norm_fil);
 }
 //_____________________________________________________________________________
-void init_nu(){
-
-  for(int tt=0;tt<24;tt++)
+void xset_fdt_fco(int argc, char *argv[]){// set fdt fcoef
+  if (argc!=3){
+    printf("\n inappropriate arguments choice in command xset_fdt_fco\n");
+    exit(99);
+  }
+  if(!strcmp(argv[1],"th")){
+    pdth_fil=atof(argv[2]);
+  }
+  else if(!strcmp(argv[1],"tz")){
+    pdtz_fil=atof(argv[2]);
+  }
+  else if(!strcmp(argv[1],"kh")){
+    fsaht_fil=atof(argv[2]);
+  }
+  else if(!strcmp(argv[1],"kz")){
+    avt_fil=atof(argv[2]);
+  }
+  else{
+    printf("\n inappropriate first argument choice in command xset_fdt_fco\n");
+    exit(99);
+  }
+}
+//_____________________________________________________________________________
+void init_nu(int argc, char *argv[]){
+  if (argc==1){
+    for(int tt=0;tt<24;tt++)
+      for(int k=0;k<NZ;k++)
+	for(int j=0;j<NY;j++)
+	  for(int i=0;i<NX;i++){
+	    YS_nu(0,i,j,k)=YS_nu(0,i,j,k)+0.1*exp(-((i+1.-(3.*NX)/4.)*(i+1.-(3.*NX)/4.)+ 
+						    (j+1.-(3.*NY)/4.)*(j+1.-(3.*NY)/4.))/4.)
+	      *(1.-exp(-double(k+1.-NZ)/NZ))*(1.-exp(-double(tt)/24.));
+	  }
+  }
+  else if  (argc==5){
     for(int k=0;k<NZ;k++)
       for(int j=0;j<NY;j++)
 	for(int i=0;i<NX;i++){
-	  YS_nu(0,i,j,k)=YS_nu(0,i,j,k)+0.1*exp(-((i+1.-(3.*NX)/4.)*(i+1.-(3.*NX)/4.)+ 
-					      (j+1.-(3.*NY)/4.)*(j+1.-(3.*NY)/4.))/4.)
-	    *(1.-exp(-double(k+1.-NZ)/NZ))*(1.-exp(-double(tt)/24.));
-	    }
+	  YS_nu(0,i,j,k)=0.;
+	  }
+
+    float v=atof(argv[1]);
+    int ji=atoi(argv[2])-1;
+    int jj=atoi(argv[3])-1;
+    int jk=atoi(argv[4])-1;
+    YS_nu(0,ji,jj,jk)=v;
+  }
+  else{
+    printf("\n inappropriate arguments choice in command init_nu\n");
+    exit(99);
+  }
 }
 //_____________________________________________________________________________
@@ -300,31 +348,40 @@
       }
   double a0,a1;
+  double a2,b2,c2;
   for(int j=1;j<NY-1;j++)
     for(int i=1;i<NX-1;i++){
       a1=avt_fil*pdtz_fil/fse3w(i,j,1);
-      tab_luzimp( i, j, 0 , 0 ) = 0.;// 0-->gam
-      tab_luzimp( i, j, 0 , 1 ) = fse3t(i,j,0)+a1;// 1-->bet
+      b2= 1.-a1/fse3t(i,j,0);
 #ifndef K_NEUM
       a0=avt_fil*pdtz_fil/fse3w(i,j,0);
-      tab_luzimp( i, j, 0 , 1 ) += a0;// 1-->bet
-#endif
+      b2 -= (a0/fse3t(i,j,0));// 1-->bet
+#endif
+      tab_luzimp( i, j, 0 , 0 ) = 0.;// 0-->gam
+      tab_luzimp( i, j, 0 , 1 ) = b2;// 1-->bet
       for(int k=1;k<NZ-2;k++){
 	a0=avt_fil*pdtz_fil/fse3w(i,j,k);
 	a1=avt_fil*pdtz_fil/fse3w(i,j,k+1);
-	tab_luzimp( i, j, k , 0 ) = -a0/tab_luzimp( i, j, k-1 , 1);
-	tab_luzimp( i, j, k , 1 ) = fse3t(i,j,k)+a1+a0 + a0 * tab_luzimp( i, j, k , 0 );
+	a2 = a0/fse3t(i,j,k);
+	b2 =  1.-a1/fse3t(i,j,k)-a0/fse3t(i,j,k);
+	c2 = a0/fse3t(i,j,k-1);
+	tab_luzimp( i, j, k , 0 ) = c2/tab_luzimp( i, j, k-1 , 1);
+	tab_luzimp( i, j, k , 1 ) = b2 - a2 * tab_luzimp( i, j, k , 0 );
       }
       a0=avt_fil*pdtz_fil/fse3w(i,j,NZ-2);
-      tab_luzimp( i, j, NZ-2 , 0 ) = -a0/tab_luzimp( i, j, NZ-3 , 1);
-      tab_luzimp( i, j, NZ-2 , 1 ) = fse3t(i,j,NZ-2)+a0 + a0 * tab_luzimp( i, j, NZ-2 , 0 );
+      a2 = a0/fse3t(i,j,NZ-2);
+      b2 =  1.-a0/fse3t(i,j,NZ-2);
 #ifndef K_NEUM
       a1=avt_fil*pdtz_fil/fse3w(i,j,NZ-1);
-      tab_luzimp( i, j, NZ-2 , 1 ) += a1;
-#endif
-    }
-  /*  for(int k=0;k<NZ;k++){ // init a 0
-    printf("gam(%d)=%f \t",k+1,tab_luzimp( 5, 5, k , 0 ));
-    printf("bet(%d)=%f \n",k+1,tab_luzimp( 5, 5, k , 1 ));
-    }*/
+      b2 -= a1/fse3t(i,j,NZ-2);
+#endif
+      c2 = a0/fse3t(i,j,NZ-3);
+     
+      tab_luzimp( i, j, NZ-2 , 0 ) = c2/tab_luzimp( i, j, NZ-3 , 1);
+      tab_luzimp( i, j, NZ-2 , 1 ) = b2 - a2 * tab_luzimp( i, j, NZ-2 , 0 );
+    }
+  for(int k=0;k<NZ;k++){ // init a 0
+    printf("gam(%d)=%f \t",k+1,tab_luzimp( 23, 15, k , 0 ));
+    printf("bet(%d)=%f \n",k+1,tab_luzimp( 23, 15, k , 1 ));
+    }
 }
 #endif
Index: /branches/branche-mb/src/ta_c.h
===================================================================
--- /branches/branche-mb/src/ta_c.h	(revision 67)
+++ /branches/branche-mb/src/ta_c.h	(revision 68)
@@ -20,5 +20,5 @@
   ////
   //  printf("norm_fac= %24.16e\n",norm_fac);
-  if(fabs(norm_fil*x1)>5.) 
+  if(fabs(norm_fil*x1)>5) 
       printf("increm(%d,%d,%d) = %24.16e (%24.16e) \n",Yi+1,Yj+1,Yk+1,norm_fil*x1,YS1);
 }
