forward (YREAL xu_d, YREAL uij, YREAL uip1j, YREAL uijp1, YREAL uijm1, YREAL uim1j, YREAL vij, YREAL vip1jm1, YREAL vijm1, YREAL vip1j, YREAL c1, YREAL c2,YREAL c3,YREAL c4,YREAL c5,YREAL c6,YREAL c7,YREAL c8,YREAL c9,YREAL mask)
{
 
 if (Yt==1) { //first step time
  //if (Yt<1000000) { //first step time  
    YS1=xu_d;
    //printf("ucd0=%f  \n",YS1);	
  }
  else
    {
      /* Boundary condition ? 
	 if(Yi==0 || Yi==nlon-1 || Yj==0 || Yj==nlat-1)
      */
      //Example: a weighted mean
      YS1=uij + dtfil*mask*(c1*uij + c2*uip1j + c3*uijp1 + c4*uijm1 + c5*uim1j + c6*vij + c7*vip1jm1 + c8*vijm1 + c9*vip1j);
     // printf("ucd=%f  \n",YS1);	
   }
}


backward (YREAL xu_d, YREAL uij, YREAL uip1j, YREAL uijp1, YREAL uijm1, YREAL uim1j, YREAL vij, YREAL vip1jm1, YREAL vijm1, YREAL vip1j, YREAL c1, YREAL c2,YREAL c3,YREAL c4,YREAL c5,YREAL c6,YREAL c7,YREAL c8,YREAL c9,YREAL mask){
  /*
  if (Yi==9 && Yj==16)
    {
      printf("YG=%g\n",YG1_u_d(Yi,Yj));
    }
  */
  
 

if (Yt==1) { //first step time
 // if (Yt<1000000) { //first step time  
    YJ1I1=1;
  } 
  else
    {
      YJ1I2=1+c1*mask*dtfil;
      YJ1I3=c2*mask*dtfil;
      YJ1I4=c3*mask*dtfil;
      YJ1I5=c4*mask*dtfil;
      YJ1I6=c5*mask*dtfil;
      YJ1I7=c6*mask*dtfil;
      YJ1I8=c7*mask*dtfil;
      YJ1I9=c8*mask*dtfil;
      YJ1I10=c9*mask*dtfil;
    }

}