source: trunk/src/sinobad.h @ 61

Last change on this file since 61 was 61, checked in by jbrlod, 13 years ago

correction du bug du LT avec solsor
  • Property svn:eol-style set to native
File size: 28.9 KB
Line 
1
2
3///////////////////////////////////////////////////////////////////////////////
4// define, globales, fonctions perso (obligatoire et autres) ...
5// mohamed.berrada@upmc.fr
6///////////////////////////////////////////////////////////////////////////////
7#include"../include/ncutil.h"//netcdf
8#define key_zco
9#include"../include/phy_cst_file.h"
10#include"../include/namelist.h"
11//#define PATH_NCFILES "/usr/home/mblod/neuro/nemo/nemo_opa/TRY/modipsl/config/GYRE/EXP00"
12#define PATH_NCFILES "../data_in"
13#include"../include/meshmask.h"
14
15// For PPCA
16/*double shfs_ta[NZ*NY*NX][NPCA];
17double moy_ta[NZ*NY*NX];
18double stdev_ta[NPCA];
19#define shfs_ta( i, j, k ,n )  (shfs_ta[(i)*(NY*NZ)+(j)*(NZ)+(k)][(n)]) // shape functions
20#define moy_ta( i, j, k )  (moy_ta[(i)*(NY*NZ)+(j)*(NZ)+(k)]) // mean
21void load_shape_func (int argc, char *argv[]);  // load PCA axes from file with mean
22void load_stdev_pac (int argc, char *argv[]);  // load stdev from file
23void load_mean (int argc, char *argv[]);  // load mean
24*/
25////#define bmask( i, j )  (tmask( i, j ,0))  // ! elliptic equation is written at t-point
26int tniter[TA+TU]; // nbre d'iterations pour backward de solsor_dynspg_flt.h
27int ndiag_rst=0;  // diagnostic restart if =1
28
29//Arrays for neuler=1 (neuler=0::Euler 1st time step)
30#define tb_neuler1( i, j, k) (tb_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)])
31#define sb_neuler1( i, j, k) (sb_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)])
32#define ub_neuler1( i, j, k) (ub_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)])
33#define vb_neuler1( i, j, k) (vb_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)])
34#define sshb_neuler1( i, j) (sshb_neuler1[(j)*(NX)+(i)])
35#define rotn_at_TU( i, j, k) (rotn_at_TU[(k)*(NY*NX)+(j)*(NX)+(i)])
36#define hdivn_at_TU( i, j, k) (hdivn_at_TU[(k)*(NY*NX)+(j)*(NX)+(i)])
37#define gcx_at_TU( i, j) (gcx_at_TU[(j)*(NX)+(i)])
38double  tb_neuler1[NZ*NY*NX]; //
39double  sb_neuler1[NZ*NY*NX]; //
40double  ub_neuler1[NZ*NY*NX]; //
41double  vb_neuler1[NZ*NY*NX]; //
42double  sshb_neuler1[NY*NX]; //
43double  rotn_at_TU[NZ*NY*NX]; //
44double  hdivn_at_TU[NZ*NY*NX]; //
45double  gcx_at_TU[NY*NX]; //
46
47void true_target_in_tab(int argc, char *argv[]); // charge YS_*a_c(0,i,j,k) dans true_*a_c(i,j,k)
48void normsup_alla_c(); // diagnostic of minimization
49double gdsr[NZ];
50int nksr;
51void xtraqsr_init();// init of solar radiation penetration
52
53#   define  rdttra( k ) rdt
54double r2dt;
55
56double   atfp1       =    1. - 2. * atfp;// !: asselin time filter coeff. (atfp1= 1-2*atfp)
57//------->zdfmxl
58double      avt_c = 5.e-4; //  ! Kz criterion for the turbocline depth
59double      rho_c = 0.01;//      ! density criterion for mixed layer depth
60
61double ppgphi0 = 29.; // latitude of first raw column T-point
62//latitude for the Coriolis or Beta parameter
63
64//mode de sauvegarde (==0 3D, ==1 2D surface)
65short savemode=0;
66
67
68int n_zdfexp=3;
69int jphgr_msh=5;//pour ff :beta-plane with regular grid-spacing and rotated domain (GYRE config)
70// a faire
71/*nyear   =   ndastp / 10000 : current year
72  nmonth    : current month of the year nyear
73  nday_year : current day of the year nyea
74  ndastp    : = nyear*10000 + nmonth*100 + nday*/
75
76void cal_ff();
77void xdisplay() ; 
78void xistate_gyre();
79void xistate_yao_file(char * pth_file_yao); // a faire
80void xdom_init();
81void xflt_rst();
82void xsolmat_init();
83void xistate_init(int argc, char *argv[]);
84void xrst_save(int argc, char *argv[]);
85
86//////////////////////////////////////////////////////////////
87// Les fonctions OBLIGATOIRES
88//_____________________________________________________________________________
89void appli_start (int argc, char *argv[])
90{       // permet si besoin de prendre la main des le debut de l'appl;
91
92  printf("////////////////////////////////////////////////////////////////////////\n");
93  printf("//                       NEMO/YAO PROJECT                             //\n");
94  printf("//                     M. Berrada 02-2009                            //\n");
95  printf("//                  LOCEAN-IPSL.UPMC (Paris 6)                        //\n");
96  printf("//====================================================================//\n");
97  printf("\n");
98#ifdef K_SOLSORYAO
99  printf("option SOLSORYAO enable\n");
100#else
101  printf("option SOLSORYAO disable\n");
102#endif
103#ifdef K_OPTIMORDER
104  printf("option OPTIMORDER enable\n");
105#else
106  printf("option OPTIMORDER disable\n");
107#endif
108
109  //Vérification du type de réel utilisé
110  //A modifier - pas génial
111#ifdef YFLOAT
112#ifndef YYFLOAT
113  printf("Incoherent real type between O_REAL option and ncutil.h\n");
114  exit(1);
115#endif
116#endif
117#ifdef YDOUBLE
118#ifndef YYDOUBLE
119  printf("Incoherent real type between O_REAL option and ncutil.h\n");
120  exit(1);
121#endif
122#endif
123
124 
125  for (int i=0;i<TA+TU;i++) tniter[i]=0;//nmax;// pour backward de solsor_dynspg_flt.h
126  phy_cst();
127  xinit_mesh_mask_nc();// lire le fichier mesh_mask.nc
128  cal_ff();
129  xdom_init();
130  xsolmat_init();
131  xflt_rst();
132  // xistate_init(); basculer dans .i
133  xdisplay();
134  xtraqsr_init();
135 
136}
137//____________________________________________________________________________
138int test=0;
139void before_it (int nit) // permet d'intervenir si besoin avant une iteration
140{
141  // static int test=0;
142  FILE *p;
143  char* c="a";
144  if(!test) c="w";
145  test=1;
146  p=fopen("cost.dat",c);
147  fprintf(p,"%23.16e\n",YTotalCost);
148  fclose(p);
149}
150//_____________________________________________________________________________
151void cost_function (int pdt)
152{                       // fonction de cout (si les standards ne conviennent pas)
153}
154//_____________________________________________________________________________
155void adjust_target ()
156{                       // fonction d'ajustement (si la standard ne convient pas)
157}
158//_____________________________________________________________________________
159void after_it (int nit)
160{
161  xdisplay();
162}
163//_____________________________________________________________________________
164void forward_before (int ctrp)
165{       // permet d'intervenir si besoin avant le forward
166  kt=Yt-TU; //oceanic time step
167 
168    if(Yt==TU+1 && neuler==0)
169      r2dt=rdt;
170    else
171      r2dt=2.*rdt; 
172  }
173//_____________________________________________________________________________
174void forward_after (int ctrp)
175{                       // permet d'intervenir si besoin aprÚs le forward
176  printf("neuler=%d\n",neuler);
177}
178//_____________________________________________________________________________
179void backward_before (int ctrp)
180{                       // permet d'intervenir si besoin avant le backward
181
182    if( ndiag_rst==1){
183      printf("backward_before:: ndiag_rst will be =0\n");
184      exit(99);
185    }
186    if(Yt==TU+1 && neuler==0)
187      r2dt=rdt;
188    else
189      r2dt=2.*rdt; 
190}
191//_____________________________________________________________________________
192void backward_after (int ctrp)
193{                       // permet d'intervenir si besoin apres le backward
194}
195//_____________________________________________________________________________
196short select_io(int indic, char *nmmod, int sortie, int iaxe, int jaxe, int kaxe, int pdt, YREAL *val)
197{                       // Pour faire des selections sur les fonctions d'entrees sorties de Yao ou autre en generale.
198                        // Doit retourner 1 si l'element dont les caracteristiques sont passes en parametre doit
199                        // etre retenu, pris, selectionne ; et 0 sinon
200                        // indic permet de savoir de quelle instance provient l'appel :
201                        // =YIO_LOADSTATE        => appel via loadstate
202                        // =YIO_SAVESTATE        => appel via savestate
203                        // =YIO_LOADOBS    => appel via loadobs
204                        // =YIO_OUTOOBS    => appel via outoobs
205
206   if(indic==YIO_SAVESTATE && !strcmp(nmmod,"tb") && savemode){
207    if(kaxe==0 )
208      return(1);
209    else 
210      return(0);
211      }
212   
213  /*  if(indic==YIO_OUTOOBS){
214    if(iaxe==9 && jaxe==9 )
215      return(1);
216    else
217      return(0);
218      }*/
219  return(1);
220}
221//======================== Fin des foncts obligatoires =========================
222//==============================================================================
223
224//Les autres fonctions
225
226//_____________________________________________________________________________
227void    xdisplay()
228{
229  normsup_alla_c();
230}
231
232//_____________________________________________________________________________
233#define hu( i, j )  (hu[(j)*(NX)+(i)])
234#define hv( i, j )  (hv[(j)*(NX)+(i)])
235#define hur( i, j )  (hur[(j)*(NX)+(i)])
236#define hvr( i, j )  (hvr[(j)*(NX)+(i)])
237
238double hu[NY*NX]; // xdom_init
239double hv[NY*NX]; // xdom_init
240double hur[NY*NX]; // xdom_init
241double hvr[NY*NX]; // xdom_init
242
243void xdom_init(){
244  for(int j=0;j<NY;j++)
245    for(int i=0;i<NX;i++){
246      hu(i,j) =0.;
247      hv(i,j) =0.;
248      for(int k=0;k<NZ;k++){
249        hu(i,j) = hu(i,j) + fse3u(i,j,k) * umask(i,j,k);
250        hv(i,j) = hv(i,j) + fse3v(i,j,k) * vmask(i,j,k);
251      }
252     
253      //      ! Inverse of the local depth
254      hur(i,j) = fse3u(i,j,0);             //! Lower bound : thickness of the first model level
255      hvr(i,j) = fse3v(i,j,0);
256      for(int k=1;k<NZ;k++){ // ! Sum of the vertical scale factors
257        hur(i,j) = hur(i,j) + fse3u(i,j,k) * umask(i,j,k);
258        hvr(i,j) = hvr(i,j) + fse3v(i,j,k) * vmask(i,j,k);
259      }
260      // ! Cte and mask the inverse of the local depth
261      hur(i,j) = 1. / hur(i,j) * umask(i,j,0);
262      hvr(i,j) = 1. / hvr(i,j) * vmask(i,j,0);
263     
264    }
265}
266
267//_____________________________________________________________________________
268#define gcp( i, j, k ,t )  (gcp[(j)*(NX)+(i)][(k)][(t)]) // a t0
269#define gcdmat( i, j ,t )  (gcdmat[(j)*(NX)+(i)][(t)])
270#define gcdprc( i, j ,t )  (gcdprc[(j)*(NX)+(i)][(t)])
271#define gccd( i, j ,t )  (gccd[(j)*(NX)+(i)][(t)])
272
273#define gcx( i, j )  (gcx[(i)][(j)])
274double gcx[NX][NY];   // xsolsor
275#define G_gcx( i , j )  (G_gcx[(j)*(NX)+(i)])
276double G_gcx[NY*NX];   // xsolsor
277#define G_gcx0( i , j )  (G_gcx0[(j)*(NX)+(i)])
278double G_gcx0[NY*NX];   // xsolsor
279#define Ytb_gcx( i , j )  (Ytb_gcx[(j)*(NX)+(i)])
280double Ytb_gcx[NY*NX];   // xsolsor
281#define gcr( i, j )  (gcr[(j)*(NX)+(i)])
282double gcr[NY*NX];   // xsolsor
283
284double gcp[NY*NX][4][2]; // xsolmat_init
285double gcdmat[NY*NX][2]; // xsolmat_init
286double gcdprc[NY*NX][2]; // xsolmat_init
287double gccd[NY*NX][2];   // xsolmat_init
288
289int nn_rstssh=0;
290void xflt_rst(){
291  for(int j=0;j<NY;j++)
292    for(int i=0;i<NX;i++){
293      YS_gcx_dynspg_flt(0,i,j,TU-1)=0.;   YS_gcx_dynspg_flt(0,i,j,TU)=0.;
294      YS_gcx2(0,i,j,TU-1)=0.;   YS_gcx2(0,i,j,TU)=0.;
295      gcx(i,j)=0.;
296    }
297   if( nn_rstssh == 1 ) {
298    for(int j=0;j<NY;j++)
299      for(int i=0;i<NX;i++){
300        YS_sshb(0,i,j,TU)=0.;   YS_sshn(0,i,j,TU)=0.;
301      }
302  }
303}
304
305//_____________________________________________________________________________
306void xsolmat_init(){
307  //      ! initialize to zero
308  double  z2dt_t0,z2dt_t1;//zcoef_t0 = 0.e0,zcoef_t1 = 0.e0;
309  for(int j=0;j<NY;j++)
310    for(int i=0;i<NX;i++){// pour neuler=0
311      gcp(i,j,0,0) = 0.e0;  gcp(i,j,0,1) = 0.e0;
312      gcp(i,j,1,0) = 0.e0;  gcp(i,j,1,1) = 0.e0;
313      gcp(i,j,2,0) = 0.e0;  gcp(i,j,2,1) = 0.e0;
314      gcp(i,j,3,0) = 0.e0;  gcp(i,j,3,1) = 0.e0;
315
316      gcdprc(i,j,0) = 0.e0;  gcdprc(i,j,1) = 0.e0;
317      gcdmat(i,j,0) = 0.e0;  gcdmat(i,j,1) = 0.e0;
318    }
319  z2dt_t0 = rdt;
320  z2dt_t1 = 2.*rdt;
321  double zcoef_t0,zcoefs_t0,zcoefw_t0,zcoefe_t0,zcoefn_t0;
322  double zcoef_t1,zcoefs_t1,zcoefw_t1,zcoefe_t1,zcoefn_t1;
323  //      ! defined the coefficients for free surface elliptic system
324  for(int j=1;j<NY-1;j++)
325    for(int i=1;i<NX-1;i++){
326      zcoef_t0 = z2dt_t0 * z2dt_t0 * grav * rnu * bmask(i,j);
327      zcoefs_t0 = -zcoef_t0 * hv(i  ,j-1) * e1v(i  ,j-1) / e2v(i  ,j-1);//    ! south coefficient
328      zcoefw_t0 = -zcoef_t0 * hu(i-1,j  ) * e2u(i-1,j  ) / e1u(i-1,j  );//    ! west coefficient
329      zcoefe_t0 = -zcoef_t0 * hu(i  ,j  ) * e2u(i  ,j  ) / e1u(i  ,j  );//    ! east coefficient
330      zcoefn_t0 = -zcoef_t0 * hv(i  ,j  ) * e1v(i  ,j  ) / e2v(i  ,j  );//    ! north coefficient
331      zcoef_t1 = z2dt_t1 * z2dt_t1 * grav * rnu * bmask(i,j);
332      zcoefs_t1 = -zcoef_t1 * hv(i  ,j-1) * e1v(i  ,j-1) / e2v(i  ,j-1);//    ! south coefficient
333      zcoefw_t1 = -zcoef_t1 * hu(i-1,j  ) * e2u(i-1,j  ) / e1u(i-1,j  );//    ! west coefficient
334      zcoefe_t1 = -zcoef_t1 * hu(i  ,j  ) * e2u(i  ,j  ) / e1u(i  ,j  );//    ! east coefficient
335      zcoefn_t1 = -zcoef_t1 * hv(i  ,j  ) * e1v(i  ,j  ) / e2v(i  ,j  );//    ! north coefficient
336
337      gcp(i,j,0,0) = zcoefs_t0;  gcp(i,j,0,1) = zcoefs_t1;
338      gcp(i,j,1,0) = zcoefw_t0;  gcp(i,j,1,1) = zcoefw_t1;
339      gcp(i,j,2,0) = zcoefe_t0;  gcp(i,j,2,1) = zcoefe_t1;
340      gcp(i,j,3,0) = zcoefn_t0;  gcp(i,j,3,1) = zcoefn_t1;
341      gcdmat(i,j,0) = e1t(i,j) * e2t(i,j) * bmask(i,j)    //         ! diagonal coefficient
342        - zcoefs_t0 -zcoefw_t0 -zcoefe_t0 -zcoefn_t0;
343      gcdmat(i,j,1) = e1t(i,j) * e2t(i,j) * bmask(i,j)    //         ! diagonal coefficient
344        - zcoefs_t1 -zcoefw_t1 -zcoefe_t1 -zcoefn_t1;
345 
346   }
347   //  ! SOR and PCG solvers
348  for(int j=0;j<NY;j++)
349    for(int i=0;i<NX;i++){
350      if( bmask(i,j) != 0 ){
351        gcdprc(i,j,0) = 1.e0 / gcdmat(i,j,0);
352        gcdprc(i,j,1) = 1.e0 / gcdmat(i,j,1);
353      }
354    }
355         
356  for(int j=0;j<NY;j++)
357    for(int i=0;i<NX;i++){
358      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);
359      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);
360      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);
361      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);
362      gccd(i,j,0) = sor * gcp(i,j,1,0); gccd(i,j,1) = sor * gcp(i,j,1,1);
363    }
364}
365
366//_____________________________________________________________________________
367void xistate_init(int argc, char *argv[]){
368  int itest=atoi(argv[1]);
369  if(itest==1){// from restart file
370    printf("start from rest_file\n ");
371    neuler=0;//=1;// Restart from a file // Set time-step indicator at nit000 (leap-frog)
372    // ifdef YE_*a_c alors *a=*b=*a_c
373    char rest_file[200];
374    if (argc!=3)
375      {
376        printf("\nncf file for restart not specified in command xistate_init\n");
377        exit(1);
378      }
379    else
380      {
381        strcpy(rest_file,argv[2]);
382      }
383    //char const    *rest_file    = PATH_NCFILES"/rst3Yao.nc";//GYRE_00000003_restart.nc";
384    int rest_file_id;
385    rest_file_id=Ouvre_nc(rest_file);
386    xistate_rest_file(rest_file_id);
387    nc_close( rest_file_id);
388    //printf("kt=%f, ndastp=%f, adatrj=%f",kt,ndastp,adatrj);
389
390    for(int i=0;i<NX;i++)
391      for(int j=0;j<NY;j++){
392        for(int k=0;k<NZ;k++){
393          YS_ua_c(0,i,j,k)=YS_ua(0,i,j,k,TU);
394          YS_va_c(0,i,j,k)=YS_va(0,i,j,k,TU);
395          YS_ta_c(0,i,j,k)=YS_ta(0,i,j,k,TU);
396          YS_sa_c(0,i,j,k)=YS_sa(0,i,j,k,TU);
397
398          tb_neuler1(i,j,k)=YS_tb(0,i,j,k,TU);
399          sb_neuler1(i,j,k)=YS_sb(0,i,j,k,TU);
400          ub_neuler1(i,j,k)=YS_ub(0,i,j,k,TU);
401          vb_neuler1(i,j,k)=YS_vb(0,i,j,k,TU);
402          rotn_at_TU( i, j, k)=YS_rotn(0,i,j,k,TU);
403          hdivn_at_TU( i, j, k)=YS_hdivn(0,i,j,k,TU);
404        }
405        YS_sshn_c(0,i,j)=YS_sshn(0,i,j,TU);     
406        sshb_neuler1(i,j)=YS_sshb(0,i,j,TU);
407        gcx_at_TU( i, j)=YS_gcx2(0,i,j,TU);
408     }
409  }
410  else if (itest==2){
411    neuler=0;//
412    xistate_gyre(); // GYRE  configuration : start from pre-defined temperature
413    //            and salinity fields
414  }
415  else if (itest==3){
416    neuler=0;//
417    xistate_yao_file("data_in/file_yao/"); // GYRE  configuration : start from pre-defined temperature
418    //            and salinity fields
419  }
420}
421
422//_____________________________________________________________________________
423void xeos_init(){
424}
425
426//_____________________________________________________________________________
427void xistate_gyre(){
428  printf("start from gyre analytic initialization\n");
429  for(int i=0;i<NX;i++){
430    for(int j=0;j<NY;j++){
431      YS_sshn(0,i,j,TU)=0.;
432      YS_sshb(0,i,j,TU)=0.;
433      YS_gcx2(0,i,j,TU)=0.;
434      YS_gcx2(0,i,j,TU-1)=0.;
435      for(int k=0;k<NZ;k++){
436        YS_ua_c(0,i,j,k)=0.;
437        YS_ua(0,i,j,k,TU)=0.;
438        YS_ub(0,i,j,k,TU)=0.;
439        YS_va_c(0,i,j,k)=0.;
440        YS_va(0,i,j,k,TU)=0.;
441        YS_vb(0,i,j,k,TU)=0.;
442        YS_rotn(0,i,j,k,TU)=0.;
443        YS_rotn(0,i,j,k,TU-1)=0.;
444        YS_hdivn(0,i,j,k,TU)=0.;
445        YS_hdivn(0,i,j,k,TU-1)=0.;
446        YS_ta_c(0,i,j,k) = (16.-12.*tanh((fsdept(i,j,k)-400.)/700.))
447          * (-tanh((500.-fsdept(i,j,k))/150.)+1.)/2.
448          + (15.*(1.-tanh((fsdept(i,j,k)-50.)/1500.))-1.4*tanh((fsdept(i,j,k)-100.)/100.) 
449             + 7.*(1500.-fsdept(i,j,k))/1500.) 
450          *(-tanh((fsdept(i,j,k)-500.)/150.)+1.)/2.;
451        YS_ta_c(0,i,j,k) = YS_ta_c(0,i,j,k)* tmask(i,j,k);
452        YS_ta(0,i,j,k,TU) = YS_ta_c(0,i,j,k);
453        YS_tb(0,i,j,k,TU) = YS_ta_c(0,i,j,k);
454       
455        YS_sa_c(0,i,j,k) = (36.25-1.13*tanh((fsdept(i,j,k)-305.)/460.))
456          *(-tanh((500.-fsdept(i,j,k))/150.)+1.)/2.
457          +(35.55+1.25*(5000.-fsdept(i,j,k))/5000.
458            -1.62*tanh((fsdept(i,j,k)-60.)/650.)
459            +0.2*tanh((fsdept(i,j,k)-35.)/100.)
460            +0.2*tanh((fsdept(i,j,k)-1000.)/5000.))
461          *(-tanh((fsdept(i,j,k)-500.)/150.)+1.)/2.;
462        YS_sa_c(0,i,j,k) = YS_sa_c(0,i,j,k)*tmask(i,j,k);
463        YS_sa(0,i,j,k,TU) = YS_sa_c(0,i,j,k);
464        YS_sb(0,i,j,k,TU) = YS_sa_c(0,i,j,k);
465
466        rotn_at_TU( i, j, k)=YS_rotn(0,i,j,k,TU);
467        hdivn_at_TU( i, j, k)=YS_hdivn(0,i,j,k,TU);
468
469      }   
470      YS_sshn_c(0,i,j)=YS_sshn(0,i,j,TU);       
471      gcx_at_TU( i, j)=YS_gcx2(0,i,j,TU);
472    }
473  }
474}
475
476void xchangesavemode(int argc, char *argv[]){
477  short newmode=atoi(argv[1]);
478  savemode=newmode;
479
480}
481
482//_____________________________________________________________________________
483void xistate_yao_file(char * pth_file_yao){ // a faire
484  /*  printf("start from yao file\n");
485  FILE *pf;
486  char *file_yao;
487  sprintf(file_yao,"%stb.dat",pth_file_yao);
488  pf=fopen(file_yao,'r');
489  fscanf(p,"%lf",);
490  fclose(p);
491 for(int j=0;j<NY;j++)
492    for(int i=0;i<NX;i++){
493      for(int k=0;k<NZ;k++){
494        YS_ua_c(0,i,j,k)=0.;
495        YS_va_c(0,i,j,k)=0.;
496        YS_ta_c(0,i,j,k) =;
497        YS_sa_c(0,i,j,k) =;
498      }   
499      YS_sshn_c(0,i,j)=;
500    }*/
501}
502//_____________________________________________________________________________
503void xtraqsr_init(){// init of solar radiation penetration
504  //! z-coordinate with or without partial step : same w-level everywhere inside the ocean
505  for(int k=0;k<NZ;k++) gdsr[k] = 0.e0;
506  for(int k=0;k<NZ;k++) {
507    double zdp1 = -gdepw_0[k]; //#define fsdepw( i , j , k )  gdepw_0[k]
508    gdsr[k] = ro0cpr*(rabs*exp(zdp1/xsi1)+(1.-rabs)*exp(zdp1/xsi2));
509    if(gdsr[k]<= 1.e-10) break;
510  }
511  int      indic = 0;
512  for(int k=0;k<NZ;k++) {
513    if( gdsr[k]<=1.e-15 && indic == 0){
514      gdsr[k] = 0.e0;
515      nksr = k+1;
516      indic = 1;
517    }
518    if(nksr>NZ-1) nksr=NZ-1;
519  }
520}                                                         
521
522//_____________________________________________________________________________
523void    cal_ff() //domhgr.F90     Coriolis parameter
524{
525  if(jphgr_msh==5){// beta-plane and rotated domain (gyre configuration)
526    double zbeta = 2.*omega*cos(rad*ppgphi0)/ra ;// beta at latitude ppgphi0
527    double zphi0 = 15.;//latitude of the first row F-points
528    double zf0   = 2.*omega*sin(rad*zphi0);// compute f0 1st point south
529    for(int i=0;i<NX;i++)
530      for(int j=0;j<NY;j++)
531        ff(i,j) = (zf0+zbeta*fabs(gphif(i,j)-zphi0)*rad*ra);// f = f0 +beta* y ( y=0 at south)
532  }
533}
534//_____________________________________________________________________________
535/*void load_shape_func (int argc, char *argv[])  // load PCA axes from file
536{
537  FILE *ffp;
538  if ((ffp=fopen( argv[2], "r")) == NULL){
539    printf( "error: could not find file of shape functions '%s'\n", argv[2]);
540    exit(99);// abort program
541  }
542  if (!strcmp(argv[1],"ta_c")){
543    for(int i=0;i<NX;i++)
544      for(int j=0;j<NY;j++)
545        for(int k=0;k<NZ;k++){
546          for(int n=0;n<NPCA;n++) shfs_ta( i, j, k ,n )=0.0;
547        }
548   
549    int ctrl=-1;
550    while (ctrl<NZ*NY*NX-1){//!feof(ffp)) {
551      ctrl=ctrl+1;
552      for(int n=0;n<NPCA;n++) fscanf(ffp,"%lf",&shfs_ta[ctrl][n]);
553    }
554  }
555  else{
556    printf( "load_shape_func:: pas encore fait\n");
557    exit(99);// abort program
558  }
559 
560  fclose(ffp);
561 
562 
563  ////////////////
564  }*/
565
566//______________________________________________________________________________
567 /*void load_stdev_pca (int argc, char *argv[])  // load stdev from file
568{
569  FILE *ffp;
570  if ((ffp=fopen( argv[2], "r")) == NULL){
571    printf( "error: could not find file of stdev '%s'\n", argv[1]);
572    exit(99);// abort program
573  }
574  if (!strcmp(argv[1],"ta_c")){
575    for(int n=0;n<NPCA;n++) fscanf(ffp,"%lf",&stdev_ta[n]);
576  }
577  else{
578    printf( "load_shape_func:: pas encore fait\n");
579    exit(99);// abort program
580  }
581  fclose(ffp);
582 
583 
584  ////////////////
585  }*/
586//______________________________________________________________________________
587  /*void load_mean (int argc, char *argv[]){  // load mean to generate obs
588  if (!strcmp(argv[2],"ta_c")){
589    if(atoi(argv[1])==0){
590      for(int i=0;i<NX;i++)
591        for(int j=0;j<NY;j++)
592          for(int k=0;k<NZ;k++){
593            moy_ta( i, j, k )=0.0;
594          }
595      for(int i=0;i<NX;i++)
596        for(int j=0;j<NY;j++)
597          for(int k=0;k<NZ;k++){
598            moy_ta( i, j, k )=YS_ta_c(0,i,j,k);
599          }
600    }
601    else if(atoi(argv[1])==1){
602      FILE *ffp;
603      if ((ffp=fopen( argv[3], "r")) == NULL){
604        printf( "error: could not find file of shape functions '%s'\n", argv[3]);
605        exit(99);// abort program
606      }
607      for(int i=0;i<NX;i++)
608        for(int j=0;j<NY;j++)
609          for(int k=0;k<NZ;k++){
610            moy_ta( i, j, k )=0.0;
611          }
612   
613      int ctrl=-1;
614      while (ctrl<NZ*NY*NX-1){//!feof(ffp)) {
615        ctrl=ctrl+1;
616        fscanf(ffp,"%lf",&moy_ta[ctrl]);
617      }
618    }
619  }
620  else{
621    printf( "load_shape_func:: pas encore fait\n");
622    exit(99);// abort program
623  }
624  }*/
625//______________________________________________________________________________
626
627//----------------------------------------------
628#if defined (YE_ta_c) || defined (YE_pca_ta)
629#define true_tac( i, j, k )   (true_tac[(k)*(NY*NX)+(j)*(NX)+(i)])
630double true_tac[NZ*NY*NX];
631#endif
632
633#ifdef YE_sa_c
634#define true_sac( i, j, k )   (true_sac[(k)*(NY*NX)+(j)*(NX)+(i)])
635double true_sac[NZ*NY*NX];
636#endif
637
638#ifdef YE_ua_c
639#define true_uac( i, j, k )   (true_uac[(k)*(NY*NX)+(j)*(NX)+(i)])
640double true_uac[NZ*NY*NX];
641#endif
642
643#ifdef YE_va_c
644#define true_vac( i, j, k )   (true_vac[(k)*(NY*NX)+(j)*(NX)+(i)])
645double true_vac[NZ*NY*NX];
646#endif
647//----------------------------------------------
648void true_target_in_tab(int argc, char *argv[]){
649  int i,j,k;
650  if(!strcmp(argv[1],"ta_c")){
651#if defined (YE_ta_c) || defined (YE_pca_ta)
652    for(i=0;i<NX;i++)
653      for(j=0;j<NY;j++)
654        for(k=0;k<NZ;k++)
655          true_tac( i, j , k ) =YS_ta_c(0,i,j,k);
656#endif
657  } 
658
659  if(!strcmp(argv[1],"sa_c")){
660#ifdef YE_sa_c
661    for(i=0;i<NX;i++)
662      for(j=0;j<NY;j++)
663        for(k=0;k<NZ;k++)
664          true_sac( i, j , k ) =YS_sa_c(0,i,j,k);
665#endif
666  } 
667
668  if(!strcmp(argv[1],"ua_c")){
669#ifdef YE_ua_c
670    for(i=0;i<NX;i++)
671      for(j=0;j<NY;j++)
672        for(k=0;k<NZ;k++)
673          true_uac( i, j , k ) =YS_ua_c(0,i,j,k);
674#endif
675  } 
676
677  if(!strcmp(argv[1],"va_c")){
678#ifdef YE_va_c
679    for(i=0;i<NX;i++)
680      for(j=0;j<NY;j++)
681        for(k=0;k<NZ;k++)
682          true_vac( i, j , k ) =YS_va_c(0,i,j,k);
683#endif
684  } 
685 
686}
687int iii=28,jjj=12;
688//---------------------------
689void normsup_alla_c(){
690  double sum=0.,sumn=0.;
691  int i,j,k;
692 
693#if defined (YE_ta_c) || defined (YE_pca_ta)
694  sum=0.;sumn=0.;
695  for(i=0;i<NX;i++)
696    for(j=0;j<NY;j++)
697      // if(i==iii && j==jjj){
698      for(k=0;k<NZ;k++){
699        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));
700        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));
701      }//}
702  sumn=sqrt(sumn);
703  printf("norm_sup_diff_ta_c : %23.16e\n",sum);
704  printf("norm_2_diff_ta_c : %23.16e\n",sumn);
705#endif
706 
707#ifdef YE_sa_c
708  sum=0.;sumn=0.;
709  for(i=0;i<NX;i++)
710    for(j=0;j<NY;j++)
711      //if(i==iii && j==jjj){
712      for(k=0;k<NZ;k++){
713        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));
714      }//}
715  printf("norm_sup_diff_sa_c : %23.16e\n",sum);
716#endif
717 
718#ifdef YE_ua_c
719  sum=0.;sumn=0.;
720  for(i=0;i<NX;i++)
721    for(j=0;j<NY;j++)
722      for(k=0;k<NZ;k++){
723        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));
724      }
725  printf("norm_sup_diff_ua_c : %23.16e\n",sum);
726#endif
727 
728#ifdef YE_va_c
729  sum=0.;sumn=0.;
730  for(i=0;i<NX;i++)
731    for(j=0;j<NY;j++)
732      for(k=0;k<NZ;k++){
733        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));
734      }
735  printf("norm_sup_diff_va_c : %23.16e\n",sum);
736#endif
737 
738}
739
740//__________________________________________________________________
741void xcomparYF();
742void xcomparYF(){
743  //comparer Fortran vs Yao
744  int i,j,k,t,ks=1,ite=0;
745  double v,sum;
746  FILE *pf;
747  char* fua="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_pua.dat";
748  if((pf=fopen(fua,"r"))==NULL){
749    printf("file in function xcomparYF\n");
750    exit(99);
751  }
752  sum=0.;
753  while(!feof(pf)){
754    fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v);
755    if(t==TA-3 && k==ks){
756      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));
757      ite=1;
758    }
759  }
760  if(ite==0)
761    printf("xcomparYF: problem dans ua\n");
762  else
763    printf("norm_sup(uaF-uaY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum);
764  fclose(pf);
765 
766  char* fva="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_pva.dat";
767  if((pf=fopen(fva,"r"))==NULL){
768    printf("file in function xcomparYF\n");
769    exit(99);
770  }
771 
772  ite=0;
773  sum=0.;
774  while(feof(pf)==0){
775    fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v);
776    if(t==TA-3 && k==ks){
777      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));
778      ite=1;
779    }
780  }
781  if(ite==0)
782    printf("xcomparYF: problem dans va\n");
783  else
784    printf("norm_sup(vaF-vaY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum);
785  fclose(pf);
786 
787  char* fwa="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_pwn.dat";
788  if((pf=fopen(fwa,"r"))==NULL){
789    printf("file in function xcomparYF\n");
790    exit(99);
791  }
792 
793  ite=0;
794  sum=0.;
795  while(!feof(pf)){
796    fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v);
797    if(t==TA-3 && k==ks){
798      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));
799      ite=1;
800    }
801  }
802  if(ite==0)
803    printf("xcomparYF: problem dans wa\n");
804  else
805    printf("norm_sup(waF-waY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum);
806  fclose(pf);
807 
808  char* fta="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_pta.dat";
809  if((pf=fopen(fta,"r"))==NULL){
810    printf("file in function xcomparYF\n");
811    exit(99);
812  }
813 
814  ite=0;
815  sum=0.;
816  while(!feof(pf)){
817    fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v);
818    if(t==TA-3 && k==ks){
819      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));
820      ite=1;
821    }
822  }
823  if(ite==0)
824    printf("xcomparYF: problem dans ta\n");
825  else
826    printf("norm_sup(taF-taY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum);
827  fclose(pf);
828 
829  char* fsa="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_psa.dat";
830  if((pf=fopen(fsa,"r"))==NULL){
831    printf("file in function xcomparYF\n");
832    exit(99);
833  }
834 
835  ite=0;
836  sum=0.;
837  while(!feof(pf)){
838    fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v);
839    if(t==TA-3 && k==ks){
840      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));
841      ite=1;
842    }
843  }
844  if(ite==0)
845    printf("xcomparYF: problem dans sa\n");
846  else
847    printf("norm_sup(saF-saY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum);
848  fclose(pf);
849
850
851  char* fsshn="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_psshn.dat";
852  if((pf=fopen(fsshn,"r"))==NULL){
853    printf("file in function xcomparYF\n");
854    exit(99);
855  }
856 
857  ite=0;
858  sum=0.;
859  while(!feof(pf)){
860    fscanf(pf," %i %i %i %lf",&t,&i,&j,&v);
861    if(t==TA-3){
862      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));
863      ite=1;
864    }
865  }
866  if(ite==0)
867    printf("xcomparYF: problem dans sshn\n");
868  else
869    printf("norm_sup(sshnF-sshnY)(t=%i) : %23.16e\n",TA,sum);
870  fclose(pf);
871 
872}
873
874//Gestion des fichiers netcdf
875
876void xrst_save(int argc, char *argv[]) {
877  //Sauvegarde un restart
878  if (argc!=2)
879    {
880      printf("\ncf file for saving restart not specified in command xrst_save\n");
881      exit(1);
882    }
883  else
884    {
885      char rest_file[200];
886      int rest_file_id;
887      int dimid[4];
888      strcpy(rest_file,argv[1]);
889      rest_file_id=Ouvre_nc_write(rest_file); //Ouvre le fichier netcds pour l'écriture
890      write_rst_global_att(rest_file_id); //Ecrit les global attributes
891      define_dim(rest_file_id,dimid); //Définit les dimensions
892      write_rst_var(TU,rest_file_id,dimid); //Ecrit les variables
893      nc_close( rest_file_id); //ferme le fichier
894    }
895}
896
897void xwriteout (int argc, char *argv[]) {
898  //xwriteout t suff ncfile
899  //Ecrit dans le fichier ncfile les variables u,v,t,s,ssh avec le suffixe suff au pas de temps t
900  if (argc!=4)
901   { 
902     printf("\n not enough arguments in writeout\n");
903     exit(1);
904   }
905
906  int t;
907  char suff[200];
908  char ncfile[200];
909  int ncfile_id;
910  t=atoi(argv[1]);
911  strcpy(suff,argv[2]);
912  strcpy(ncfile,argv[3]);
913  ncfile_id=Ouvre_nc_add(ncfile);
914  write_out_var(t,ncfile_id,suff);
915  nc_close(ncfile_id);
916
917}
918void xwritegrad (int argc, char *argv[]) {
919  //xwritegrad t suff ncfile
920  //Ecrit dans le fichier ncfile les grad des variables u,v,t,s,ssh avec le suffixe suff au pas de temps t
921  if (argc!=4)
922   { 
923     printf("\n not enough arguments in writeout\n");
924     exit(1);
925   }
926
927  int t;
928  char suff[200];
929  char ncfile[200];
930  int ncfile_id;
931  t=atoi(argv[1]);
932  strcpy(suff,argv[2]);
933  strcpy(ncfile,argv[3]);
934  ncfile_id=Ouvre_nc_add(ncfile);
935  write_grad_var(t,ncfile_id,suff);
936  nc_close(ncfile_id);
937
938}
939
940void xinitnc(int argc,char *argv[]) {
941  //xinit ncfile
942  //init le fichier netcdf pour sauvegarde
943  if (argc!=2)
944   { 
945     printf("\ncf file for saving output not specified in command xinitnc\n");
946     exit(1);
947   }
948  else
949    {
950      char ncfile[200];
951      int nc_file_id;
952      int dimid[4];
953      strcpy(ncfile,argv[1]);
954      nc_file_id=Ouvre_nc_write(ncfile);
955      define_dim(nc_file_id,dimid);
956      write_out_var_init(nc_file_id,dimid);
957      nc_close(nc_file_id);
958    }
959 
960}
Note: See TracBrowser for help on using the repository browser.