source: trunk/src/sinobad.h @ 42

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

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