source: trunk/src/sinobad.h @ 6

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

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