source: trunk/src/sinobad.h @ 64

///////////////////////////////////////////////////////////////////////////////
// define, globales, fonctions perso (obligatoire et autres) ...
// mohamed.berrada@upmc.fr
///////////////////////////////////////////////////////////////////////////////
#include"../include/ncutil.h"//netcdf
#define key_zco
#include"../include/phy_cst_file.h"
#include"../include/namelist.h"
//#define PATH_NCFILES "/usr/home/mblod/neuro/nemo/nemo_opa/TRY/modipsl/config/GYRE/EXP00"
#define PATH_NCFILES "../data_in"
#include"../include/meshmask.h"

// For PPCA 
/*double shfs_ta[NZ*NY*NX][NPCA];
double moy_ta[NZ*NY*NX];
double stdev_ta[NPCA];
#define shfs_ta( i, j, k ,n )  (shfs_ta[(i)*(NY*NZ)+(j)*(NZ)+(k)][(n)]) // shape functions
#define moy_ta( i, j, k )  (moy_ta[(i)*(NY*NZ)+(j)*(NZ)+(k)]) // mean
void load_shape_func (int argc, char *argv[]);  // load PCA axes from file with mean
void load_stdev_pac (int argc, char *argv[]);  // load stdev from file
void load_mean (int argc, char *argv[]);  // load mean 
*/
////#define bmask( i, j )  (tmask( i, j ,0))  // ! elliptic equation is written at t-point
int tniter[TA+TU]; // nbre d'iterations pour backward de solsor_dynspg_flt.h
int ndiag_rst=0;  // diagnostic restart if =1

//Arrays for neuler=1 (neuler=0::Euler 1st time step)
#define tb_neuler1( i, j, k) (tb_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)]) 
#define sb_neuler1( i, j, k) (sb_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)]) 
#define ub_neuler1( i, j, k) (ub_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)]) 
#define vb_neuler1( i, j, k) (vb_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)]) 
#define sshb_neuler1( i, j) (sshb_neuler1[(j)*(NX)+(i)]) 
#define rotn_at_TU( i, j, k) (rotn_at_TU[(k)*(NY*NX)+(j)*(NX)+(i)]) 
#define hdivn_at_TU( i, j, k) (hdivn_at_TU[(k)*(NY*NX)+(j)*(NX)+(i)]) 
#define gcx_at_TU( i, j) (gcx_at_TU[(j)*(NX)+(i)]) 
double  tb_neuler1[NZ*NY*NX]; //
double  sb_neuler1[NZ*NY*NX]; //
double  ub_neuler1[NZ*NY*NX]; //
double  vb_neuler1[NZ*NY*NX]; //
double  sshb_neuler1[NY*NX]; //
double  rotn_at_TU[NZ*NY*NX]; //
double  hdivn_at_TU[NZ*NY*NX]; //
double  gcx_at_TU[NY*NX]; //

void true_target_in_tab(int argc, char *argv[]); // charge YS_*a_c(0,i,j,k) dans true_*a_c(i,j,k)
void normsup_alla_c(); // diagnostic of minimization
double gdsr[NZ];
int nksr;
void xtraqsr_init();// init of solar radiation penetration
double Ytold; //pour voir se pas de temps
double ktinit; //premier pas ocean
#   define  rdttra( k ) rdt
double r2dt;

double   atfp1       =    1. - 2. * atfp;// !: asselin time filter coeff. (atfp1= 1-2*atfp)
//------->zdfmxl
double      avt_c = 5.e-4; //  ! Kz criterion for the turbocline depth
double      rho_c = 0.01;//      ! density criterion for mixed layer depth

double ppgphi0 = 29.; // latitude of first raw column T-point
//latitude for the Coriolis or Beta parameter

//mode de sauvegarde (==0 3D, ==1 2D surface)
short savemode=0;


int n_zdfexp=3;
int jphgr_msh=5;//pour ff :beta-plane with regular grid-spacing and rotated domain (GYRE config)
// a faire
/*nyear   =   ndastp / 10000 : current year
  nmonth    : current month of the year nyear
  nday_year : current day of the year nyea
  ndastp    : = nyear*10000 + nmonth*100 + nday*/

void cal_ff();
void xdisplay() ; 
void xistate_gyre();
void xistate_yao_file(char * pth_file_yao); // a faire
void xdom_init();
void xflt_rst();
void xsolmat_init();
void xistate_init(int argc, char *argv[]);
void xrst_save(int argc, char *argv[]);

//////////////////////////////////////////////////////////////
// Les fonctions OBLIGATOIRES
//_____________________________________________________________________________
void appli_start (int argc, char *argv[])
{       // permet si besoin de prendre la main des le debut de l'appl;

  printf("////////////////////////////////////////////////////////////////////////\n");
  printf("//                       NEMO/YAO PROJECT                             //\n");
  printf("//                     M. Berrada 02-2009                            //\n");
  printf("//                  LOCEAN-IPSL.UPMC (Paris 6)                        //\n");
  printf("//====================================================================//\n");
  printf("\n");
#ifdef K_SOLSORYAO
  printf("option SOLSORYAO enable\n");
#else
  printf("option SOLSORYAO disable\n");
#endif
#ifdef K_OPTIMORDER
  printf("option OPTIMORDER enable\n");
#else
  printf("option OPTIMORDER disable\n");
#endif

  //Vérification du type de réel utilisé
  //A modifier - pas génial
#ifdef YFLOAT
#ifndef YYFLOAT
  printf("Incoherent real type between O_REAL option and ncutil.h\n");
  exit(1);
#endif
#endif
#ifdef YDOUBLE
#ifndef YYDOUBLE
  printf("Incoherent real type between O_REAL option and ncutil.h\n");
  exit(1);
#endif
#endif

  
  for (int i=0;i<TA+TU;i++) tniter[i]=0;//nmax;// pour backward de solsor_dynspg_flt.h
  phy_cst();
  xinit_mesh_mask_nc();// lire le fichier mesh_mask.nc
  cal_ff();
  xdom_init();
  xsolmat_init();
  xflt_rst();
  // xistate_init(); basculer dans .i
  xdisplay();
  xtraqsr_init();
  ktinit=0; //A mettre dans la namelist ????
}
//____________________________________________________________________________
int test=0;
void before_it (int nit) // permet d'intervenir si besoin avant une iteration
{
  // static int test=0;
  FILE *p;
  char* c="a";
  if(!test) c="w";
  test=1;
  p=fopen("cost.dat",c);
  fprintf(p,"%23.16e\n",YTotalCost);
  fclose(p);
  Ytold=Yt;
  kt=ktinit;
}
//_____________________________________________________________________________
void cost_function (int pdt)
{                       // fonction de cout (si les standards ne conviennent pas)
}
//_____________________________________________________________________________
void adjust_target ()
{                       // fonction d'ajustement (si la standard ne convient pas)
}
//_____________________________________________________________________________
void after_it (int nit)
{
  xdisplay();
}
//_____________________________________________________________________________
void forward_before (int ctrp)
{       // permet d'intervenir si besoin avant le forward
 

  kt=ktinit + (Yt-TU);
  printf("kt=%lf\n",kt);
    if(Yt==TU+1 && neuler==0)
      r2dt=rdt;
    else
      r2dt=2.*rdt;  
  }
//_____________________________________________________________________________
void forward_after (int ctrp)
{                       // permet d'intervenir si besoin aprÚs le forward
  //printf("neuler=%d\n",neuler);
}
//_____________________________________________________________________________
void backward_before (int ctrp)
{                       // permet d'intervenir si besoin avant le backward

    if( ndiag_rst==1){
      printf("backward_before:: ndiag_rst will be =0\n");
      exit(99);
    }
    if(Yt==TU+1 && neuler==0)
      r2dt=rdt;
    else
      r2dt=2.*rdt;  
}
//_____________________________________________________________________________
void backward_after (int ctrp)
{                       // permet d'intervenir si besoin apres le backward
}
//_____________________________________________________________________________
short select_io(int indic, char *nmmod, int sortie, int iaxe, int jaxe, int kaxe, int pdt, YREAL *val)
{                       // Pour faire des selections sur les fonctions d'entrees sorties de Yao ou autre en generale.
                        // Doit retourner 1 si l'element dont les caracteristiques sont passes en parametre doit
                        // etre retenu, pris, selectionne ; et 0 sinon
                        // indic permet de savoir de quelle instance provient l'appel :
                        // =YIO_LOADSTATE        => appel via loadstate
                        // =YIO_SAVESTATE        => appel via savestate
                        // =YIO_LOADOBS    => appel via loadobs
                        // =YIO_OUTOOBS    => appel via outoobs

   if(indic==YIO_SAVESTATE && !strcmp(nmmod,"tb") && savemode){
    if(kaxe==0 )
      return(1);
    else 
      return(0);
      }
    
  /*  if(indic==YIO_OUTOOBS){
    if(iaxe==9 && jaxe==9 )
      return(1);
    else 
      return(0);
      }*/
  return(1);
}
//======================== Fin des foncts obligatoires =========================
//==============================================================================

//Les autres fonctions

//_____________________________________________________________________________
void    xdisplay()
{
  normsup_alla_c();
}

//_____________________________________________________________________________
#define hu( i, j )  (hu[(j)*(NX)+(i)])
#define hv( i, j )  (hv[(j)*(NX)+(i)])
#define hur( i, j )  (hur[(j)*(NX)+(i)])
#define hvr( i, j )  (hvr[(j)*(NX)+(i)])

double hu[NY*NX]; // xdom_init
double hv[NY*NX]; // xdom_init
double hur[NY*NX]; // xdom_init
double hvr[NY*NX]; // xdom_init

void xdom_init(){
  for(int j=0;j<NY;j++)
    for(int i=0;i<NX;i++){
      hu(i,j) =0.;
      hv(i,j) =0.;
      for(int k=0;k<NZ;k++){
        hu(i,j) = hu(i,j) + fse3u(i,j,k) * umask(i,j,k);
        hv(i,j) = hv(i,j) + fse3v(i,j,k) * vmask(i,j,k);
      }
      
      //      ! Inverse of the local depth
      hur(i,j) = fse3u(i,j,0);             //! Lower bound : thickness of the first model level
      hvr(i,j) = fse3v(i,j,0);
      for(int k=1;k<NZ;k++){ // ! Sum of the vertical scale factors
        hur(i,j) = hur(i,j) + fse3u(i,j,k) * umask(i,j,k);
        hvr(i,j) = hvr(i,j) + fse3v(i,j,k) * vmask(i,j,k);
      }
      // ! Cte and mask the inverse of the local depth
      hur(i,j) = 1. / hur(i,j) * umask(i,j,0);
      hvr(i,j) = 1. / hvr(i,j) * vmask(i,j,0);
      
    }
}

//_____________________________________________________________________________
#define gcp( i, j, k ,t )  (gcp[(j)*(NX)+(i)][(k)][(t)]) // a t0
#define gcdmat( i, j ,t )  (gcdmat[(j)*(NX)+(i)][(t)])
#define gcdprc( i, j ,t )  (gcdprc[(j)*(NX)+(i)][(t)])
#define gccd( i, j ,t )  (gccd[(j)*(NX)+(i)][(t)])

#define gcx( i, j )  (gcx[(i)][(j)])
double gcx[NX][NY];   // xsolsor
#define G_gcx( i , j )  (G_gcx[(j)*(NX)+(i)])
double G_gcx[NY*NX];   // xsolsor
#define G_gcx0( i , j )  (G_gcx0[(j)*(NX)+(i)])
double G_gcx0[NY*NX];   // xsolsor
#define Ytb_gcx( i , j )  (Ytb_gcx[(j)*(NX)+(i)])
double Ytb_gcx[NY*NX];   // xsolsor
#define gcr( i, j )  (gcr[(j)*(NX)+(i)])
double gcr[NY*NX];   // xsolsor

double gcp[NY*NX][4][2]; // xsolmat_init
double gcdmat[NY*NX][2]; // xsolmat_init
double gcdprc[NY*NX][2]; // xsolmat_init
double gccd[NY*NX][2];   // xsolmat_init

int nn_rstssh=0;
void xflt_rst(){
  for(int j=0;j<NY;j++)
    for(int i=0;i<NX;i++){
      YS_gcx_dynspg_flt(0,i,j,TU-1)=0.;   YS_gcx_dynspg_flt(0,i,j,TU)=0.;
      YS_gcx2(0,i,j,TU-1)=0.;   YS_gcx2(0,i,j,TU)=0.;
      gcx(i,j)=0.;
    }
   if( nn_rstssh == 1 ) {
    for(int j=0;j<NY;j++)
      for(int i=0;i<NX;i++){
        YS_sshb(0,i,j,TU)=0.;   YS_sshn(0,i,j,TU)=0.;
      }
  }
}

//_____________________________________________________________________________
void xsolmat_init(){
  //      ! initialize to zero
  double  z2dt_t0,z2dt_t1;//zcoef_t0 = 0.e0,zcoef_t1 = 0.e0;
  for(int j=0;j<NY;j++)
    for(int i=0;i<NX;i++){// pour neuler=0
      gcp(i,j,0,0) = 0.e0;  gcp(i,j,0,1) = 0.e0;
      gcp(i,j,1,0) = 0.e0;  gcp(i,j,1,1) = 0.e0;
      gcp(i,j,2,0) = 0.e0;  gcp(i,j,2,1) = 0.e0;
      gcp(i,j,3,0) = 0.e0;  gcp(i,j,3,1) = 0.e0;

      gcdprc(i,j,0) = 0.e0;  gcdprc(i,j,1) = 0.e0;
      gcdmat(i,j,0) = 0.e0;  gcdmat(i,j,1) = 0.e0;
    }
  z2dt_t0 = rdt;
  z2dt_t1 = 2.*rdt;
  double zcoef_t0,zcoefs_t0,zcoefw_t0,zcoefe_t0,zcoefn_t0;
  double zcoef_t1,zcoefs_t1,zcoefw_t1,zcoefe_t1,zcoefn_t1;
  //      ! defined the coefficients for free surface elliptic system
  for(int j=1;j<NY-1;j++)
    for(int i=1;i<NX-1;i++){
      zcoef_t0 = z2dt_t0 * z2dt_t0 * grav * rnu * bmask(i,j);
      zcoefs_t0 = -zcoef_t0 * hv(i  ,j-1) * e1v(i  ,j-1) / e2v(i  ,j-1);//    ! south coefficient
      zcoefw_t0 = -zcoef_t0 * hu(i-1,j  ) * e2u(i-1,j  ) / e1u(i-1,j  );//    ! west coefficient
      zcoefe_t0 = -zcoef_t0 * hu(i  ,j  ) * e2u(i  ,j  ) / e1u(i  ,j  );//    ! east coefficient
      zcoefn_t0 = -zcoef_t0 * hv(i  ,j  ) * e1v(i  ,j  ) / e2v(i  ,j  );//    ! north coefficient
      zcoef_t1 = z2dt_t1 * z2dt_t1 * grav * rnu * bmask(i,j);
      zcoefs_t1 = -zcoef_t1 * hv(i  ,j-1) * e1v(i  ,j-1) / e2v(i  ,j-1);//    ! south coefficient
      zcoefw_t1 = -zcoef_t1 * hu(i-1,j  ) * e2u(i-1,j  ) / e1u(i-1,j  );//    ! west coefficient
      zcoefe_t1 = -zcoef_t1 * hu(i  ,j  ) * e2u(i  ,j  ) / e1u(i  ,j  );//    ! east coefficient
      zcoefn_t1 = -zcoef_t1 * hv(i  ,j  ) * e1v(i  ,j  ) / e2v(i  ,j  );//    ! north coefficient

      gcp(i,j,0,0) = zcoefs_t0;  gcp(i,j,0,1) = zcoefs_t1;
      gcp(i,j,1,0) = zcoefw_t0;  gcp(i,j,1,1) = zcoefw_t1;
      gcp(i,j,2,0) = zcoefe_t0;  gcp(i,j,2,1) = zcoefe_t1;
      gcp(i,j,3,0) = zcoefn_t0;  gcp(i,j,3,1) = zcoefn_t1;
      gcdmat(i,j,0) = e1t(i,j) * e2t(i,j) * bmask(i,j)    //         ! diagonal coefficient
        - zcoefs_t0 -zcoefw_t0 -zcoefe_t0 -zcoefn_t0;
      gcdmat(i,j,1) = e1t(i,j) * e2t(i,j) * bmask(i,j)    //         ! diagonal coefficient
        - zcoefs_t1 -zcoefw_t1 -zcoefe_t1 -zcoefn_t1;
 
   }
   //  ! SOR and PCG solvers
  for(int j=0;j<NY;j++)
    for(int i=0;i<NX;i++){
      if( bmask(i,j) != 0 ){
        gcdprc(i,j,0) = 1.e0 / gcdmat(i,j,0);
        gcdprc(i,j,1) = 1.e0 / gcdmat(i,j,1);
      }
    }
         
  for(int j=0;j<NY;j++)
    for(int i=0;i<NX;i++){
      gcp(i,j,0,0) = gcp(i,j,0,0) * gcdprc(i,j,0); gcp(i,j,0,1) = gcp(i,j,0,1) * gcdprc(i,j,1);
      gcp(i,j,1,0) = gcp(i,j,1,0) * gcdprc(i,j,0); gcp(i,j,1,1) = gcp(i,j,1,1) * gcdprc(i,j,1);
      gcp(i,j,2,0) = gcp(i,j,2,0) * gcdprc(i,j,0); gcp(i,j,2,1) = gcp(i,j,2,1) * gcdprc(i,j,1);
      gcp(i,j,3,0) = gcp(i,j,3,0) * gcdprc(i,j,0); gcp(i,j,3,1) = gcp(i,j,3,1) * gcdprc(i,j,1);
      gccd(i,j,0) = sor * gcp(i,j,1,0); gccd(i,j,1) = sor * gcp(i,j,1,1);
    }
}

//_____________________________________________________________________________
void xistate_init(int argc, char *argv[]){
  int itest=atoi(argv[1]);
  if(itest==1){// from restart file
    printf("start from rest_file\n ");
    neuler=0;//=1;// Restart from a file // Set time-step indicator at nit000 (leap-frog)
    // ifdef YE_*a_c alors *a=*b=*a_c
    char rest_file[200];
    if (argc!=3)
      {
        printf("\nncf file for restart not specified in command xistate_init\n");
        exit(1);
      }
    else
      {
        strcpy(rest_file,argv[2]);
      }
    //char const    *rest_file    = PATH_NCFILES"/rst3Yao.nc";//GYRE_00000003_restart.nc";
    int rest_file_id;
    rest_file_id=Ouvre_nc(rest_file);
    xistate_rest_file(rest_file_id);
    nc_close( rest_file_id);
    //printf("kt=%f, ndastp=%f, adatrj=%f",kt,ndastp,adatrj);

    for(int i=0;i<NX;i++)
      for(int j=0;j<NY;j++){
        for(int k=0;k<NZ;k++){
          YS_ua_c(0,i,j,k)=YS_ua(0,i,j,k,TU);
          YS_va_c(0,i,j,k)=YS_va(0,i,j,k,TU);
          YS_ta_c(0,i,j,k)=YS_ta(0,i,j,k,TU);
          YS_sa_c(0,i,j,k)=YS_sa(0,i,j,k,TU);

          tb_neuler1(i,j,k)=YS_tb(0,i,j,k,TU);
          sb_neuler1(i,j,k)=YS_sb(0,i,j,k,TU);
          ub_neuler1(i,j,k)=YS_ub(0,i,j,k,TU);
          vb_neuler1(i,j,k)=YS_vb(0,i,j,k,TU);
          rotn_at_TU( i, j, k)=YS_rotn(0,i,j,k,TU);
          hdivn_at_TU( i, j, k)=YS_hdivn(0,i,j,k,TU);
        }
        YS_sshn_c(0,i,j)=YS_sshn(0,i,j,TU);     
        sshb_neuler1(i,j)=YS_sshb(0,i,j,TU);
        gcx_at_TU( i, j)=YS_gcx2(0,i,j,TU);
     }
  }
  else if (itest==2){
    neuler=0;// 
    xistate_gyre(); // GYRE  configuration : start from pre-defined temperature
    //            and salinity fields 
  }
  else if (itest==3){
    neuler=0;// 
    xistate_yao_file("data_in/file_yao/"); // GYRE  configuration : start from pre-defined temperature
    //            and salinity fields 
  }
}

//_____________________________________________________________________________
void xeos_init(){
}

//_____________________________________________________________________________
void xistate_gyre(){
  printf("start from gyre analytic initialization\n");
  for(int i=0;i<NX;i++){
    for(int j=0;j<NY;j++){
      YS_sshn(0,i,j,TU)=0.;
      YS_sshb(0,i,j,TU)=0.;
      YS_gcx2(0,i,j,TU)=0.;
      YS_gcx2(0,i,j,TU-1)=0.;
      for(int k=0;k<NZ;k++){
        YS_ua_c(0,i,j,k)=0.;
        YS_ua(0,i,j,k,TU)=0.;
        YS_ub(0,i,j,k,TU)=0.;
        YS_va_c(0,i,j,k)=0.;
        YS_va(0,i,j,k,TU)=0.;
        YS_vb(0,i,j,k,TU)=0.;
        YS_rotn(0,i,j,k,TU)=0.;
        YS_rotn(0,i,j,k,TU-1)=0.;
        YS_hdivn(0,i,j,k,TU)=0.;
        YS_hdivn(0,i,j,k,TU-1)=0.;
        YS_ta_c(0,i,j,k) = (16.-12.*tanh((fsdept(i,j,k)-400.)/700.))
          * (-tanh((500.-fsdept(i,j,k))/150.)+1.)/2.
          + (15.*(1.-tanh((fsdept(i,j,k)-50.)/1500.))-1.4*tanh((fsdept(i,j,k)-100.)/100.)  
             + 7.*(1500.-fsdept(i,j,k))/1500.) 
          *(-tanh((fsdept(i,j,k)-500.)/150.)+1.)/2.;
        YS_ta_c(0,i,j,k) = YS_ta_c(0,i,j,k)* tmask(i,j,k);
        YS_ta(0,i,j,k,TU) = YS_ta_c(0,i,j,k);
        YS_tb(0,i,j,k,TU) = YS_ta_c(0,i,j,k);
        
        YS_sa_c(0,i,j,k) = (36.25-1.13*tanh((fsdept(i,j,k)-305.)/460.))
          *(-tanh((500.-fsdept(i,j,k))/150.)+1.)/2.
          +(35.55+1.25*(5000.-fsdept(i,j,k))/5000.
            -1.62*tanh((fsdept(i,j,k)-60.)/650.)
            +0.2*tanh((fsdept(i,j,k)-35.)/100.)
            +0.2*tanh((fsdept(i,j,k)-1000.)/5000.))
          *(-tanh((fsdept(i,j,k)-500.)/150.)+1.)/2.;
        YS_sa_c(0,i,j,k) = YS_sa_c(0,i,j,k)*tmask(i,j,k);
        YS_sa(0,i,j,k,TU) = YS_sa_c(0,i,j,k);
        YS_sb(0,i,j,k,TU) = YS_sa_c(0,i,j,k);

        rotn_at_TU( i, j, k)=YS_rotn(0,i,j,k,TU);
        hdivn_at_TU( i, j, k)=YS_hdivn(0,i,j,k,TU);

      }   
      YS_sshn_c(0,i,j)=YS_sshn(0,i,j,TU);       
      gcx_at_TU( i, j)=YS_gcx2(0,i,j,TU);
    }
  }
}

void xchangesavemode(int argc, char *argv[]){
  short newmode=atoi(argv[1]);
  savemode=newmode;

}

//_____________________________________________________________________________
void xistate_yao_file(char * pth_file_yao){ // a faire
  /*  printf("start from yao file\n");
  FILE *pf;
  char *file_yao;
  sprintf(file_yao,"%stb.dat",pth_file_yao);
  pf=fopen(file_yao,'r');
  fscanf(p,"%lf",);
  fclose(p);
 for(int j=0;j<NY;j++)
    for(int i=0;i<NX;i++){
      for(int k=0;k<NZ;k++){
        YS_ua_c(0,i,j,k)=0.;
        YS_va_c(0,i,j,k)=0.;
        YS_ta_c(0,i,j,k) =;
        YS_sa_c(0,i,j,k) =;
      }   
      YS_sshn_c(0,i,j)=;
    }*/
}
//_____________________________________________________________________________
void xtraqsr_init(){// init of solar radiation penetration
  //! z-coordinate with or without partial step : same w-level everywhere inside the ocean
  for(int k=0;k<NZ;k++) gdsr[k] = 0.e0;
  for(int k=0;k<NZ;k++) {
    double zdp1 = -gdepw_0[k]; //#define fsdepw( i , j , k )  gdepw_0[k]
    gdsr[k] = ro0cpr*(rabs*exp(zdp1/xsi1)+(1.-rabs)*exp(zdp1/xsi2));
    if(gdsr[k]<= 1.e-10) break;
  }
  int      indic = 0;
  for(int k=0;k<NZ;k++) {
    if( gdsr[k]<=1.e-15 && indic == 0){
      gdsr[k] = 0.e0;
      nksr = k+1;
      indic = 1;
    }
    if(nksr>NZ-1) nksr=NZ-1;
  }
}                                                          

//_____________________________________________________________________________
void    cal_ff() //domhgr.F90     Coriolis parameter
{
  if(jphgr_msh==5){// beta-plane and rotated domain (gyre configuration)
    double zbeta = 2.*omega*cos(rad*ppgphi0)/ra ;// beta at latitude ppgphi0
    double zphi0 = 15.;//latitude of the first row F-points
    double zf0   = 2.*omega*sin(rad*zphi0);// compute f0 1st point south
    for(int i=0;i<NX;i++)
      for(int j=0;j<NY;j++)
        ff(i,j) = (zf0+zbeta*fabs(gphif(i,j)-zphi0)*rad*ra);// f = f0 +beta* y ( y=0 at south)
  }
}
//_____________________________________________________________________________
/*void load_shape_func (int argc, char *argv[])  // load PCA axes from file
{
  FILE *ffp; 
  if ((ffp=fopen( argv[2], "r")) == NULL){
    printf( "error: could not find file of shape functions '%s'\n", argv[2]);
    exit(99);// abort program
  }
  if (!strcmp(argv[1],"ta_c")){
    for(int i=0;i<NX;i++)
      for(int j=0;j<NY;j++)
        for(int k=0;k<NZ;k++){
          for(int n=0;n<NPCA;n++) shfs_ta( i, j, k ,n )=0.0;
        }
    
    int ctrl=-1;
    while (ctrl<NZ*NY*NX-1){//!feof(ffp)) {
      ctrl=ctrl+1;
      for(int n=0;n<NPCA;n++) fscanf(ffp,"%lf",&shfs_ta[ctrl][n]);
    }
  }
  else{
    printf( "load_shape_func:: pas encore fait\n");
    exit(99);// abort program
  }
  
  fclose(ffp);
  
  
  ////////////////
  }*/

//______________________________________________________________________________
 /*void load_stdev_pca (int argc, char *argv[])  // load stdev from file
{
  FILE *ffp; 
  if ((ffp=fopen( argv[2], "r")) == NULL){
    printf( "error: could not find file of stdev '%s'\n", argv[1]);
    exit(99);// abort program
  }
  if (!strcmp(argv[1],"ta_c")){
    for(int n=0;n<NPCA;n++) fscanf(ffp,"%lf",&stdev_ta[n]);
  }
  else{
    printf( "load_shape_func:: pas encore fait\n");
    exit(99);// abort program
  }
  fclose(ffp);
  
  
  ////////////////
  }*/
//______________________________________________________________________________
  /*void load_mean (int argc, char *argv[]){  // load mean to generate obs
  if (!strcmp(argv[2],"ta_c")){
    if(atoi(argv[1])==0){
      for(int i=0;i<NX;i++)
        for(int j=0;j<NY;j++)
          for(int k=0;k<NZ;k++){
            moy_ta( i, j, k )=0.0;
          }
      for(int i=0;i<NX;i++)
        for(int j=0;j<NY;j++)
          for(int k=0;k<NZ;k++){
            moy_ta( i, j, k )=YS_ta_c(0,i,j,k);
          }
    }
    else if(atoi(argv[1])==1){
      FILE *ffp; 
      if ((ffp=fopen( argv[3], "r")) == NULL){
        printf( "error: could not find file of shape functions '%s'\n", argv[3]);
        exit(99);// abort program
      }
      for(int i=0;i<NX;i++)
        for(int j=0;j<NY;j++)
          for(int k=0;k<NZ;k++){
            moy_ta( i, j, k )=0.0;
          }
    
      int ctrl=-1;
      while (ctrl<NZ*NY*NX-1){//!feof(ffp)) {
        ctrl=ctrl+1;
        fscanf(ffp,"%lf",&moy_ta[ctrl]);
      }
    }
  }
  else{
    printf( "load_shape_func:: pas encore fait\n");
    exit(99);// abort program
  }
  }*/
//______________________________________________________________________________

//----------------------------------------------
#if defined (YE_ta_c) || defined (YE_pca_ta)
#define true_tac( i, j, k )   (true_tac[(k)*(NY*NX)+(j)*(NX)+(i)])
double true_tac[NZ*NY*NX];
#endif

#ifdef YE_sa_c
#define true_sac( i, j, k )   (true_sac[(k)*(NY*NX)+(j)*(NX)+(i)])
double true_sac[NZ*NY*NX];
#endif

#ifdef YE_ua_c
#define true_uac( i, j, k )   (true_uac[(k)*(NY*NX)+(j)*(NX)+(i)])
double true_uac[NZ*NY*NX];
#endif

#ifdef YE_va_c
#define true_vac( i, j, k )   (true_vac[(k)*(NY*NX)+(j)*(NX)+(i)])
double true_vac[NZ*NY*NX];
#endif
//----------------------------------------------
void true_target_in_tab(int argc, char *argv[]){
  int i,j,k;
  if(!strcmp(argv[1],"ta_c")){
#if defined (YE_ta_c) || defined (YE_pca_ta)
    for(i=0;i<NX;i++)
      for(j=0;j<NY;j++)
        for(k=0;k<NZ;k++)
          true_tac( i, j , k ) =YS_ta_c(0,i,j,k);
#endif
  } 

  if(!strcmp(argv[1],"sa_c")){
#ifdef YE_sa_c
    for(i=0;i<NX;i++)
      for(j=0;j<NY;j++)
        for(k=0;k<NZ;k++)
          true_sac( i, j , k ) =YS_sa_c(0,i,j,k);
#endif
  } 

  if(!strcmp(argv[1],"ua_c")){
#ifdef YE_ua_c
    for(i=0;i<NX;i++)
      for(j=0;j<NY;j++)
        for(k=0;k<NZ;k++)
          true_uac( i, j , k ) =YS_ua_c(0,i,j,k);
#endif
  } 

  if(!strcmp(argv[1],"va_c")){
#ifdef YE_va_c
    for(i=0;i<NX;i++)
      for(j=0;j<NY;j++)
        for(k=0;k<NZ;k++)
          true_vac( i, j , k ) =YS_va_c(0,i,j,k);
#endif
  } 
  
}
int iii=28,jjj=12;
//---------------------------
void normsup_alla_c(){
  double sum=0.,sumn=0.;
  int i,j,k;
  
#if defined (YE_ta_c) || defined (YE_pca_ta)
  sum=0.;sumn=0.;
  for(i=0;i<NX;i++)
    for(j=0;j<NY;j++)
      // if(i==iii && j==jjj){
      for(k=0;k<NZ;k++){
        if(sum<fabs(true_tac( i, j , k )-YS_ta_c(0,i,j,k))) sum=fabs(true_tac( i, j , k )-YS_ta_c(0,i,j,k));
        sumn+=(true_tac( i, j , k )-YS_ta_c(0,i,j,k))*(true_tac( i, j , k )-YS_ta_c(0,i,j,k));
      }//}
  sumn=sqrt(sumn);
  printf("norm_sup_diff_ta_c : %23.16e\n",sum);
  printf("norm_2_diff_ta_c : %23.16e\n",sumn);
#endif
  
#ifdef YE_sa_c
  sum=0.;sumn=0.;
  for(i=0;i<NX;i++)
    for(j=0;j<NY;j++)
      //if(i==iii && j==jjj){
      for(k=0;k<NZ;k++){
        if(sum<fabs(true_sac( i, j , k )-YS_sa_c(0,i,j,k))) sum=fabs(true_sac( i, j , k )-YS_sa_c(0,i,j,k));
      }//}
  printf("norm_sup_diff_sa_c : %23.16e\n",sum);
#endif
  
#ifdef YE_ua_c
  sum=0.;sumn=0.;
  for(i=0;i<NX;i++)
    for(j=0;j<NY;j++)
      for(k=0;k<NZ;k++){
        if(sum<fabs(true_uac( i, j , k )-YS_ua_c(0,i,j,k))) sum=fabs(true_uac( i, j , k )-YS_ua_c(0,i,j,k));
      }
  printf("norm_sup_diff_ua_c : %23.16e\n",sum);
#endif
  
#ifdef YE_va_c
  sum=0.;sumn=0.;
  for(i=0;i<NX;i++)
    for(j=0;j<NY;j++)
      for(k=0;k<NZ;k++){
        if(sum<fabs(true_vac( i, j , k )-YS_va_c(0,i,j,k))) sum=fabs(true_vac( i, j , k )-YS_va_c(0,i,j,k));
      }
  printf("norm_sup_diff_va_c : %23.16e\n",sum);
#endif
  
}

//__________________________________________________________________
void xcomparYF();
void xcomparYF(){
  //comparer Fortran vs Yao
  int i,j,k,t,ks=1,ite=0;
  double v,sum;
  FILE *pf;
  char* fua="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_pua.dat";
  if((pf=fopen(fua,"r"))==NULL){
    printf("file in function xcomparYF\n");
    exit(99);
  }
  sum=0.;
  while(!feof(pf)){
    fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v);
    if(t==TA-3 && k==ks){
      if(sum<fabs(v-YS_ua(0,i-1,j-1,k-1,TA-1))) sum=fabs(v-YS_ua(0,i-1,j-1,k-1,TA-1));
      ite=1;
    }
  }
  if(ite==0)
    printf("xcomparYF: problem dans ua\n");
  else
    printf("norm_sup(uaF-uaY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum);
  fclose(pf);
  
  char* fva="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_pva.dat";
  if((pf=fopen(fva,"r"))==NULL){
    printf("file in function xcomparYF\n");
    exit(99);
  }
  
  ite=0;
  sum=0.;
  while(feof(pf)==0){
    fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v);
    if(t==TA-3 && k==ks){
      if(sum<fabs(v-YS_va(0,i-1,j-1,k-1,TA-1))) sum=fabs(v-YS_va(0,i-1,j-1,k-1,TA-1));
      ite=1;
    }
  }
  if(ite==0)
    printf("xcomparYF: problem dans va\n");
  else
    printf("norm_sup(vaF-vaY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum);
  fclose(pf);
  
  char* fwa="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_pwn.dat";
  if((pf=fopen(fwa,"r"))==NULL){
    printf("file in function xcomparYF\n");
    exit(99);
  }
  
  ite=0;
  sum=0.;
  while(!feof(pf)){
    fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v);
    if(t==TA-3 && k==ks){
      if(sum<fabs(v-YS_wa(0,i-1,j-1,k-1,TA-1))) sum=fabs(v-YS_wa(0,i-1,j-1,k-1,TA-1));
      ite=1;
    }
  }
  if(ite==0)
    printf("xcomparYF: problem dans wa\n");
  else
    printf("norm_sup(waF-waY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum);
  fclose(pf);
  
  char* fta="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_pta.dat";
  if((pf=fopen(fta,"r"))==NULL){
    printf("file in function xcomparYF\n");
    exit(99);
  }
  
  ite=0;
  sum=0.;
  while(!feof(pf)){
    fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v);
    if(t==TA-3 && k==ks){
      if(sum<fabs(v-YS_ta(0,i-1,j-1,k-1,TA-1))) sum=fabs(v-YS_ta(0,i-1,j-1,k-1,TA-1));
      ite=1;
    }
  }
  if(ite==0)
    printf("xcomparYF: problem dans ta\n");
  else
    printf("norm_sup(taF-taY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum);
  fclose(pf);
  
  char* fsa="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_psa.dat";
  if((pf=fopen(fsa,"r"))==NULL){
    printf("file in function xcomparYF\n");
    exit(99);
  }
  
  ite=0;
  sum=0.;
  while(!feof(pf)){
    fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v);
    if(t==TA-3 && k==ks){
      if(sum<fabs(v-YS_sa(0,i-1,j-1,k-1,TA-1))) sum=fabs(v-YS_sa(0,i-1,j-1,k-1,TA-1));
      ite=1;
    }
  }
  if(ite==0)
    printf("xcomparYF: problem dans sa\n");
  else
    printf("norm_sup(saF-saY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum);
  fclose(pf);


  char* fsshn="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_psshn.dat";
  if((pf=fopen(fsshn,"r"))==NULL){
    printf("file in function xcomparYF\n");
    exit(99);
  }
  
  ite=0;
  sum=0.;
  while(!feof(pf)){
    fscanf(pf," %i %i %i %lf",&t,&i,&j,&v);
    if(t==TA-3){
      if(sum<fabs(v-YS_sshn(0,i-1,j-1,TA-1))) sum=fabs(v-YS_sshn(0,i-1,j-1,TA-1));
      ite=1;
    }
  }
  if(ite==0)
    printf("xcomparYF: problem dans sshn\n");
  else
    printf("norm_sup(sshnF-sshnY)(t=%i) : %23.16e\n",TA,sum);
  fclose(pf);
  
}

//Gestion des fichiers netcdf

void xrst_save(int argc, char *argv[]) {
  //Sauvegarde un restart
  if (argc!=2)
    {
      printf("\ncf file for saving restart not specified in command xrst_save\n");
      exit(1);
    }
  else
    {
      char rest_file[200];
      int rest_file_id;
      int dimid[4];
      strcpy(rest_file,argv[1]);
      rest_file_id=Ouvre_nc_write(rest_file); //Ouvre le fichier netcds pour l'écriture
      write_rst_global_att(rest_file_id); //Ecrit les global attributes
      define_dim(rest_file_id,dimid); //Définit les dimensions
      write_rst_var(TU,rest_file_id,dimid); //Ecrit les variables
      nc_close( rest_file_id); //ferme le fichier
    }
}

void xwriteout (int argc, char *argv[]) {
  //xwriteout t suff ncfile 
  //Ecrit dans le fichier ncfile les variables u,v,t,s,ssh avec le suffixe suff au pas de temps t
  if (argc!=4)
   { 
     printf("\n not enough arguments in writeout\n");
     exit(1);
   }

  int t;
  char suff[200];
  char ncfile[200];
  int ncfile_id;
  t=atoi(argv[1]);
  strcpy(suff,argv[2]);
  strcpy(ncfile,argv[3]);
  ncfile_id=Ouvre_nc_add(ncfile);
  write_out_var(t,ncfile_id,suff);
  nc_close(ncfile_id);

}
void xwritegrad (int argc, char *argv[]) {
  //xwritegrad t suff ncfile 
  //Ecrit dans le fichier ncfile les grad des variables u,v,t,s,ssh avec le suffixe suff au pas de temps t
  if (argc!=4)
   { 
     printf("\n not enough arguments in writeout\n");
     exit(1);
   }

  int t;
  char suff[200];
  char ncfile[200];
  int ncfile_id;
  t=atoi(argv[1]);
  strcpy(suff,argv[2]);
  strcpy(ncfile,argv[3]);
  ncfile_id=Ouvre_nc_add(ncfile);
  write_grad_var(t,ncfile_id,suff);
  nc_close(ncfile_id);

}

void xinitnc(int argc,char *argv[]) {
  //xinit ncfile
  //init le fichier netcdf pour sauvegarde
  if (argc!=2)
   { 
     printf("\ncf file for saving output not specified in command xinitnc\n");
     exit(1);
   }
  else
    {
      char ncfile[200];
      int nc_file_id;
      int dimid[4];
      strcpy(ncfile,argv[1]);
      nc_file_id=Ouvre_nc_write(ncfile);
      define_dim(nc_file_id,dimid);
      write_out_var_init(nc_file_id,dimid);
      nc_close(nc_file_id);
    }
  
}

void init_euler(int argc, char *argv[]){
  if (argc!=2)
    {
      printf("\nwrong number of argument in init_euler : specify 0 or 1\n");
      exit(1);
    }
  else
    {
      int n;
      n=atoi(argv[1]);
      if(n!=0 && n!=1)
        {
          printf("\nwrong argument in init_euler : specify 0 or 1 in command init_euler\n");
          exit(1);
        }
      neuler=n;
    }
}