source: trunk/src/filtre.h @ 186

/***************************************************************************
                 module classe filtre_lap.h  -  description
***************************************************************************/
// Mohamed Berrada
// locean-ipsl.upmc, Paris, July 19, 2011
//===========================================================================
//                          methode forward
#define dta_c( i, j, k) (dta_c[(k)*(NY*NX)+(j)*(NX)+(i)]) // l'increment de t calculer par le fitre
#define ztu( i, j, k) (ztu[(k)*(NY*NX)+(j)*(NX)+(i)]) 
#define ztv( i, j, k) (ztv[(k)*(NY*NX)+(j)*(NX)+(i)]) 
#define zwy( i, j, k) (zwy[(k)*(NY*NX)+(j)*(NX)+(i)]) 

double  dta_c[NZ*NY*NX]; // l'increment de t calculer par le fitre
double  ztu[NZ*NY*NX]; //
double  ztv[NZ*NY*NX]; //
double  zwy[NZ*NY*NX]; //

forward ()
{
  int n,ji,jj,jk;
  double zabe1,zabe2,zbtr,zta;
  double zave3r;
  double fsahtu_fil; // coef diff hori
  double fsahtv_fil; // coef diff hori
  double avt_fil; // coef diffu verti
  double pdtv_fil, pdtz_fil; // pas du temps verti et horiz
  //boundarie conditions 
  for(jk=0;jk<NZ;jk++)
    for(jj=0;jj<NY;jj++){
      YS1_nu(0,jj,jk)=0.;      YS1_nu(NX-1,jj,jk)=0.;
  }
  for(jk=0;jk<NZ;jk++)
    for(ji=0;ji<NX;ji++){
      YS1_nu(ji,0,jk)=0.;      YS1_nu(ji,NY-1,jk)=0.;
  }
  //end boundarie conditions

  // initialisation de l'increment et des tableaux temp
  for(jk=0;jk<NZ;jk++)
    for(jj=0;jj<NY;jj++)
      for(ji=0;ji<NX;ji++){
        dta_c(ji,jj,jk)=YS1_nu(ji,jj,jk);
        ztu(ji,jj,jk)=0.;
        ztv(ji,jj,jk)=0.;
        zwy(ji,jj,jk)=0.;
      }


  // sqrt W^(-1/2) a faire ici

  // sqrt horiz diffu Lh^1/2
  for(int n=0;n<MHs2;n++){  
    for(jk = 0; jk<NZ-1,jk++){                                  
      //! 1. First derivative (gradient)
      for( jj = 1; jj<NY-1;jj++){
        for(ji = 1; ji<NX-1;ji++){
          zabe1 = fsahtu_fil * umask(ji,jj,jk) * ze1ur(ji,jj);
          zabe2 = fsahtv_fil * vmask(ji,jj,jk) * ze2vr(ji,jj);
          ztu(ji,jj,jk) = zabe1 * ( dta_c(ji+1,jj  ,jk) - dta_c(ji,jj,jk) );
          ztv(ji,jj,jk) = zabe2 * ( dta_c(ji  ,jj+1,jk) - dta_c(ji,jj,jk) );
        }
      }  
      //! 2. Second derivative (divergence)
      //! --------------------
      for(jj = 1; jj<NY-1; jj++){
        for(ji = 1; ji<NX-1;ji++){
          zbtr = zbtr2(ji,jj);
          //! horizontal diffusive trends
          zta = pdth_fil*zbtr*(ztu(ji,jj,jk)- ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk));
          //! add it to the general tracer trends
          dta_c(ji,jj,jk) = dta_c(ji,jj,jk) + zta;
        }
      }
    }
  }

  // sqrt verti diffu Lv^1/2
  for(n=0;n<MZs2;n++){   
    for(jk=1;jk<NZ;jk++){ //! first vertical derivative (gradient)
      for(jj = 1;jj<NY-1;jj++){
        for(ji = 1;ji<NX-1;ji++){ 
          zave3r = 1.e0 / fse3w(ji,jj,jk); //! 1 sur le pas veritcal au point w (Arakawa mesh)
          zwy(ji,jj,jk) =   avt_fil * ( dta_c(ji,jj,jk-1) - dta_c(ji,jj,jk) ) * zave3r;
        }
      }
    }
    
    for(jk=0;jk<NZ-1;jk++){ //! second vertical derivative (divergence)
        for(jj=1;jj<NY-1;jj++){ 
          for(ji=1;ji<NX-1;ji++){ 
            ze3tr=1.e0 / fse3t(ji,jj,jk);// ! 1 sur le pas veritcal au point t 
            zta = pdtz_fil* ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) * ze3tr;
            dta_c(ji,jj,jk) = dta_c(ji,jj,jk) + zta;
          }
        }
    }
  }
  // 
}

//===========================================================================
//                          methode backward and linward
backward ()
{

  if(Ycurward==BACKWARD){
  }
  else if(Ycurward==LINWARD){
  }
}