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 | #ifdef K_FILTER |
---|
16 | #ifdef K_FZIMP |
---|
17 | void luzimp_init(); // initialise tab_luzimp |
---|
18 | #define tab_luzimp( i, j, k , t ) (tab_luzimp[(i)*(NY*NZ)+(j)*(NZ)+(k)][t]) // factorisation LU |
---|
19 | double tab_luzimp[NZ*NY*NX][2]; // tableau factorisation LU |
---|
20 | #endif |
---|
21 | #define ta_eb( i, j, k ) (ta_eb[(i)*(NY*NZ)+(j)*(NZ)+(k)]) // ebauche |
---|
22 | double ta_eb[NZ*NY*NX]; // ebauche |
---|
23 | void load_eb (int argc, char *argv[]); // load "ebauche" from current variable ta_c |
---|
24 | #endif |
---|
25 | // For PPCA |
---|
26 | /*double shfs_ta[NZ*NY*NX][NPCA]; |
---|
27 | double moy_ta[NZ*NY*NX]; |
---|
28 | double stdev_ta[NPCA]; |
---|
29 | #define shfs_ta( i, j, k ,n ) (shfs_ta[(i)*(NY*NZ)+(j)*(NZ)+(k)][(n)]) // shape functions |
---|
30 | #define moy_ta( i, j, k ) (moy_ta[(i)*(NY*NZ)+(j)*(NZ)+(k)]) // mean |
---|
31 | void load_shape_func (int argc, char *argv[]); // load PCA axes from file with mean |
---|
32 | void load_stdev_pac (int argc, char *argv[]); // load stdev from file |
---|
33 | void load_mean (int argc, char *argv[]); // load mean |
---|
34 | */ |
---|
35 | ////#define bmask( i, j ) (tmask( i, j ,0)) // ! elliptic equation is written at t-point |
---|
36 | int tniter[TA+TU]; // nbre d'iterations pour backward de solsor_dynspg_flt.h |
---|
37 | int ndiag_rst=0; // diagnostic restart if =1 |
---|
38 | |
---|
39 | //Arrays for neuler=1 (neuler=0::Euler 1st time step) |
---|
40 | #define tb_neuler1( i, j, k) (tb_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)]) |
---|
41 | #define sb_neuler1( i, j, k) (sb_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)]) |
---|
42 | #define ub_neuler1( i, j, k) (ub_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)]) |
---|
43 | #define vb_neuler1( i, j, k) (vb_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)]) |
---|
44 | #define sshb_neuler1( i, j) (sshb_neuler1[(j)*(NX)+(i)]) |
---|
45 | #define rotn_at_TU( i, j, k) (rotn_at_TU[(k)*(NY*NX)+(j)*(NX)+(i)]) |
---|
46 | #define hdivn_at_TU( i, j, k) (hdivn_at_TU[(k)*(NY*NX)+(j)*(NX)+(i)]) |
---|
47 | #define gcx_at_TU( i, j) (gcx_at_TU[(j)*(NX)+(i)]) |
---|
48 | double tb_neuler1[NZ*NY*NX]; // |
---|
49 | double sb_neuler1[NZ*NY*NX]; // |
---|
50 | double ub_neuler1[NZ*NY*NX]; // |
---|
51 | double vb_neuler1[NZ*NY*NX]; // |
---|
52 | double sshb_neuler1[NY*NX]; // |
---|
53 | double rotn_at_TU[NZ*NY*NX]; // |
---|
54 | double hdivn_at_TU[NZ*NY*NX]; // |
---|
55 | double gcx_at_TU[NY*NX]; // |
---|
56 | |
---|
57 | void true_target_in_tab(int argc, char *argv[]); // charge YS_*a_c(0,i,j,k) dans true_*a_c(i,j,k) |
---|
58 | void normsup_alla_c(); // diagnostic of minimization |
---|
59 | double gdsr[NZ]; |
---|
60 | int nksr; |
---|
61 | void xtraqsr_init();// init of solar radiation penetration |
---|
62 | |
---|
63 | |
---|
64 | # define rdttra( k ) rdt |
---|
65 | double r2dt; |
---|
66 | |
---|
67 | double atfp1 = 1. - 2. * atfp;// !: asselin time filter coeff. (atfp1= 1-2*atfp) |
---|
68 | //------->zdfmxl |
---|
69 | double avt_c = 5.e-4; // ! Kz criterion for the turbocline depth |
---|
70 | double rho_c = 0.01;// ! density criterion for mixed layer depth |
---|
71 | |
---|
72 | double ppgphi0 = 29.; // latitude of first raw column T-point |
---|
73 | //latitude for the Coriolis or Beta parameter |
---|
74 | |
---|
75 | //mode de sauvegarde (==0 3D, ==1 2D surface ==2 un seul point) |
---|
76 | short savemode=0; |
---|
77 | short isave=23; |
---|
78 | short jsave=16; |
---|
79 | short ksave=0; |
---|
80 | |
---|
81 | int n_zdfexp=3; |
---|
82 | int jphgr_msh=5;//pour ff :beta-plane with regular grid-spacing and rotated domain (GYRE config) |
---|
83 | // a faire |
---|
84 | /*nyear = ndastp / 10000 : current year |
---|
85 | nmonth : current month of the year nyear |
---|
86 | nday_year : current day of the year nyea |
---|
87 | ndastp : = nyear*10000 + nmonth*100 + nday*/ |
---|
88 | |
---|
89 | #ifdef K_FILTER |
---|
90 | double pdth_fil=100.;//7200.; |
---|
91 | double pdtz_fil=100.;//100.;//7200.; |
---|
92 | double fsaht_fil= 5.e+02;//1000.; |
---|
93 | double avt_fil=1.2e-05;//1.2e-5; |
---|
94 | //double lscale= 212000.; // length-scale |
---|
95 | double norm_fil; //coef de normalisation |
---|
96 | void normfil();// calcul du coef de normalisation norm_fil |
---|
97 | void xset_fdt_fco();// set fdt fcoef |
---|
98 | #endif |
---|
99 | |
---|
100 | void cal_ff(); |
---|
101 | void xdisplay() ; |
---|
102 | void xistate_gyre(); |
---|
103 | void xistate_yao_file(char * pth_file_yao); // a faire |
---|
104 | void xdom_init(); |
---|
105 | void xflt_rst(); |
---|
106 | void xsolmat_init(); |
---|
107 | void xistate_init(int argc, char *argv[]); |
---|
108 | void xrst_save(int argc, char *argv[]); |
---|
109 | |
---|
110 | ////////////////////////////////////////////////////////////// |
---|
111 | // Les fonctions OBLIGATOIRES |
---|
112 | //_____________________________________________________________________________ |
---|
113 | void appli_start (int argc, char *argv[]) |
---|
114 | { // permet si besoin de prendre la main des le debut de l'appl; |
---|
115 | |
---|
116 | printf("////////////////////////////////////////////////////////////////////////\n"); |
---|
117 | printf("// NEMO/YAO PROJECT //\n"); |
---|
118 | printf("// M. Berrada 02-2009 //\n"); |
---|
119 | printf("// LOCEAN-IPSL.UPMC (Paris 6) //\n"); |
---|
120 | printf("//====================================================================//\n"); |
---|
121 | printf("\n"); |
---|
122 | #ifdef K_FILTER |
---|
123 | printf("option FILTER enable\n"); |
---|
124 | #else |
---|
125 | printf("option FILTER disable\n"); |
---|
126 | #endif |
---|
127 | #ifdef K_SOLSORYAO |
---|
128 | printf("option SOLSORYAO enable\n"); |
---|
129 | #else |
---|
130 | printf("option SOLSORYAO disable\n"); |
---|
131 | #endif |
---|
132 | #ifdef K_OPTIMORDER |
---|
133 | printf("option OPTIMORDER enable\n"); |
---|
134 | #else |
---|
135 | printf("option OPTIMORDER disable\n"); |
---|
136 | #endif |
---|
137 | |
---|
138 | //Vérification du type de réel utilisé |
---|
139 | //A modifier - pas génial |
---|
140 | #ifdef YFLOAT |
---|
141 | #ifndef YYFLOAT |
---|
142 | printf("Incoherent real type between O_REAL option and ncutil.h\n"); |
---|
143 | exit(1); |
---|
144 | #endif |
---|
145 | #endif |
---|
146 | #ifdef YDOUBLE |
---|
147 | #ifndef YYDOUBLE |
---|
148 | printf("Incoherent real type between O_REAL option and ncutil.h\n"); |
---|
149 | exit(1); |
---|
150 | #endif |
---|
151 | #endif |
---|
152 | |
---|
153 | |
---|
154 | for (int i=0;i<TA+TU;i++) tniter[i]=0;//nmax;// pour backward de solsor_dynspg_flt.h |
---|
155 | phy_cst(); |
---|
156 | xinit_mesh_mask_nc();// lire le fichier mesh_mask.nc |
---|
157 | cal_ff(); |
---|
158 | xdom_init(); |
---|
159 | #ifdef K_FILTER |
---|
160 | #ifdef K_FZIMP |
---|
161 | luzimp_init(); |
---|
162 | #endif |
---|
163 | normfil(); |
---|
164 | #endif |
---|
165 | xsolmat_init(); |
---|
166 | xflt_rst(); |
---|
167 | // xistate_init(); basculer dans .i |
---|
168 | xdisplay(); |
---|
169 | xtraqsr_init(); |
---|
170 | YDispTime=1; |
---|
171 | } |
---|
172 | //____________________________________________________________________________ |
---|
173 | int test=0; |
---|
174 | void before_it (int nit) // permet d'intervenir si besoin avant une iteration |
---|
175 | { |
---|
176 | // static int test=0; |
---|
177 | FILE *p; |
---|
178 | char* c="a"; |
---|
179 | if(!test) c="w"; |
---|
180 | test=1; |
---|
181 | p=fopen("cost.dat",c); |
---|
182 | fprintf(p,"%23.16e\n",YTotalCost); |
---|
183 | fclose(p); |
---|
184 | //kt=nit000+(Yt-TU); |
---|
185 | } |
---|
186 | //_____________________________________________________________________________ |
---|
187 | void cost_function (int pdt) |
---|
188 | { // fonction de cout (si les standards ne conviennent pas) |
---|
189 | } |
---|
190 | //_____________________________________________________________________________ |
---|
191 | void adjust_target () |
---|
192 | { // fonction d'ajustement (si la standard ne convient pas) |
---|
193 | } |
---|
194 | //_____________________________________________________________________________ |
---|
195 | void after_it (int nit) |
---|
196 | { |
---|
197 | xdisplay(); |
---|
198 | } |
---|
199 | //_____________________________________________________________________________ |
---|
200 | void forward_before (int ctrp) |
---|
201 | { // permet d'intervenir si besoin avant le forward |
---|
202 | kt=nit000 + (Yt-TU); |
---|
203 | if(YDispTime) printf(" >>> current forward time : kt=%lf \n",kt); |
---|
204 | if(Yt==TU+1 && neuler==0) |
---|
205 | r2dt=rdt; |
---|
206 | else |
---|
207 | r2dt=2.*rdt; |
---|
208 | } |
---|
209 | //_____________________________________________________________________________ |
---|
210 | void forward_after (int ctrp) |
---|
211 | { |
---|
212 | // permet d'intervenir si besoin aprÚs le forward |
---|
213 | //printf("neuler=%d\n",neuler); |
---|
214 | } |
---|
215 | //_____________________________________________________________________________ |
---|
216 | void backward_before (int ctrp) |
---|
217 | { // permet d'intervenir si besoin avant le backward |
---|
218 | |
---|
219 | if( ndiag_rst==1){ |
---|
220 | printf("backward_before:: ndiag_rst will be =0\n"); |
---|
221 | exit(99); |
---|
222 | } |
---|
223 | if(Yt==TU+1 && neuler==0) |
---|
224 | r2dt=rdt; |
---|
225 | else |
---|
226 | r2dt=2.*rdt; |
---|
227 | } |
---|
228 | //_____________________________________________________________________________ |
---|
229 | void backward_after (int ctrp) |
---|
230 | { // permet d'intervenir si besoin apres le backward |
---|
231 | } |
---|
232 | //_____________________________________________________________________________ |
---|
233 | short select_io(int indic, char *nmmod, int sortie, int iaxe, int jaxe, int kaxe, int pdt, YREAL *val) |
---|
234 | { // Pour faire des selections sur les fonctions d'entrees sorties de Yao ou autre en generale. |
---|
235 | // Doit retourner 1 si l'element dont les caracteristiques sont passes en parametre doit |
---|
236 | // etre retenu, pris, selectionne ; et 0 sinon |
---|
237 | // indic permet de savoir de quelle instance provient l'appel : |
---|
238 | // =YIO_LOADSTATE => appel via loadstate |
---|
239 | // =YIO_SAVESTATE => appel via savestate |
---|
240 | // =YIO_LOADOBS => appel via loadobs |
---|
241 | // =YIO_OUTOOBS => appel via outoobs |
---|
242 | |
---|
243 | if(indic==YIO_SAVESTATE && !strcmp(nmmod,"tb") && savemode==1){ |
---|
244 | if(kaxe==0 ) |
---|
245 | return(1); |
---|
246 | else |
---|
247 | return(0); |
---|
248 | } |
---|
249 | if(indic==YIO_SAVESTATE && !strcmp(nmmod,"tb") && savemode==2){ |
---|
250 | if(iaxe==isave && jaxe==jsave && kaxe==ksave) |
---|
251 | return(1); |
---|
252 | else |
---|
253 | return(0); |
---|
254 | } |
---|
255 | |
---|
256 | |
---|
257 | /* if(indic==YIO_OUTOOBS){ |
---|
258 | if(iaxe==9 && jaxe==9 ) |
---|
259 | return(1); |
---|
260 | else |
---|
261 | return(0); |
---|
262 | }*/ |
---|
263 | return(1); |
---|
264 | } |
---|
265 | //======================== Fin des foncts obligatoires ========================= |
---|
266 | //============================================================================== |
---|
267 | |
---|
268 | //Les autres fonctions |
---|
269 | |
---|
270 | //_____________________________________________________________________________ |
---|
271 | void xdisplay() |
---|
272 | { |
---|
273 | normsup_alla_c(); |
---|
274 | } |
---|
275 | |
---|
276 | //_____________________________________________________________________________ |
---|
277 | #ifdef K_FILTER |
---|
278 | //_____________________________________________________________________________ |
---|
279 | void normfil(){ |
---|
280 | double norm_zfil=1.; |
---|
281 | double norm_lfil=1.; |
---|
282 | norm_lfil= 2.*sqrt(PI*fsaht_fil*pdth_fil*2.*(OFTL-1.));// schema explicite |
---|
283 | norm_zfil= 2.*sqrt(PI*avt_fil*pdtz_fil*2.*(OFTZ-1.));// schema explicite |
---|
284 | #ifdef K_FZIMP // schema implicite |
---|
285 | double fact=1.; |
---|
286 | int M=2*(OFTZ-1); |
---|
287 | for(int i=1;i<M;i++) |
---|
288 | fact = 4.*fact * (double) i / (double) (M+i-1); |
---|
289 | norm_zfil=2.*fact*sqrt(avt_fil*pdtz_fil); |
---|
290 | #endif |
---|
291 | |
---|
292 | #if defined (K_FILZNOL) |
---|
293 | norm_lfil=1.; |
---|
294 | #endif |
---|
295 | #if defined (K_FILLNOZ) |
---|
296 | norm_zfil=1.; |
---|
297 | #endif |
---|
298 | |
---|
299 | norm_fil= norm_zfil*norm_lfil; |
---|
300 | printf("norm_fil = %23.16e\n",norm_fil); |
---|
301 | } |
---|
302 | //_____________________________________________________________________________ |
---|
303 | void xset_fdt_fco(int argc, char *argv[]){// set fdt fcoef |
---|
304 | if (argc!=3){ |
---|
305 | printf("\n inappropriate arguments choice in command xset_fdt_fco\n"); |
---|
306 | exit(99); |
---|
307 | } |
---|
308 | if(!strcmp(argv[1],"th")){ |
---|
309 | pdth_fil=atof(argv[2]); |
---|
310 | } |
---|
311 | else if(!strcmp(argv[1],"tz")){ |
---|
312 | pdtz_fil=atof(argv[2]); |
---|
313 | } |
---|
314 | else if(!strcmp(argv[1],"kh")){ |
---|
315 | fsaht_fil=atof(argv[2]); |
---|
316 | } |
---|
317 | else if(!strcmp(argv[1],"kz")){ |
---|
318 | avt_fil=atof(argv[2]); |
---|
319 | } |
---|
320 | else{ |
---|
321 | printf("\n inappropriate first argument choice in command xset_fdt_fco\n"); |
---|
322 | exit(99); |
---|
323 | } |
---|
324 | } |
---|
325 | //_____________________________________________________________________________ |
---|
326 | void init_nu(int argc, char *argv[]){ |
---|
327 | if (argc==1){ |
---|
328 | for(int tt=0;tt<24;tt++) |
---|
329 | for(int k=0;k<NZ;k++) |
---|
330 | for(int j=0;j<NY;j++) |
---|
331 | for(int i=0;i<NX;i++){ |
---|
332 | YS_nu(0,i,j,k)=YS_nu(0,i,j,k)+0.1*exp(-((i+1.-(3.*NX)/4.)*(i+1.-(3.*NX)/4.)+ |
---|
333 | (j+1.-(3.*NY)/4.)*(j+1.-(3.*NY)/4.))/4.) |
---|
334 | *(1.-exp(-double(k+1.-NZ)/NZ))*(1.-exp(-double(tt)/24.)); |
---|
335 | } |
---|
336 | } |
---|
337 | else if (argc==5){ |
---|
338 | for(int k=0;k<NZ;k++) |
---|
339 | for(int j=0;j<NY;j++) |
---|
340 | for(int i=0;i<NX;i++){ |
---|
341 | YS_nu(0,i,j,k)=0.; |
---|
342 | } |
---|
343 | |
---|
344 | float v=atof(argv[1]); |
---|
345 | int ji=atoi(argv[2])-1; |
---|
346 | int jj=atoi(argv[3])-1; |
---|
347 | int jk=atoi(argv[4])-1; |
---|
348 | YS_nu(0,ji,jj,jk)=v; |
---|
349 | } |
---|
350 | else{ |
---|
351 | printf("\n inappropriate arguments choice in command init_nu\n"); |
---|
352 | exit(99); |
---|
353 | } |
---|
354 | } |
---|
355 | //_____________________________________________________________________________ |
---|
356 | #ifdef K_FZIMP |
---|
357 | //_____________________________________________________________________________ |
---|
358 | void luzimp_init(){ |
---|
359 | for(int j=0;j<NY;j++) |
---|
360 | for(int i=0;i<NX;i++) |
---|
361 | for(int k=0;k<NZ;k++){ // init a 0 |
---|
362 | tab_luzimp( i, j, k , 0 ) = 0.;// 0-->gam |
---|
363 | tab_luzimp( i, j, k , 1 ) = 0.;// 1-->bet |
---|
364 | } |
---|
365 | double a0,a1; |
---|
366 | double a2,b2,c2; |
---|
367 | for(int j=1;j<NY-1;j++) |
---|
368 | for(int i=1;i<NX-1;i++){ |
---|
369 | a1=avt_fil*pdtz_fil/fse3w(i,j,1); |
---|
370 | b2= 1.-a1/fse3t(i,j,0); |
---|
371 | #ifndef K_NEUM |
---|
372 | a0=avt_fil*pdtz_fil/fse3w(i,j,0); |
---|
373 | b2 -= (a0/fse3t(i,j,0));// 1-->bet |
---|
374 | #endif |
---|
375 | tab_luzimp( i, j, 0 , 0 ) = 0.;// 0-->gam |
---|
376 | tab_luzimp( i, j, 0 , 1 ) = b2;// 1-->bet |
---|
377 | for(int k=1;k<NZ-2;k++){ |
---|
378 | a0=avt_fil*pdtz_fil/fse3w(i,j,k); |
---|
379 | a1=avt_fil*pdtz_fil/fse3w(i,j,k+1); |
---|
380 | a2 = a0/fse3t(i,j,k); |
---|
381 | b2 = 1.-a1/fse3t(i,j,k)-a0/fse3t(i,j,k); |
---|
382 | c2 = a0/fse3t(i,j,k-1); |
---|
383 | tab_luzimp( i, j, k , 0 ) = c2/tab_luzimp( i, j, k-1 , 1); |
---|
384 | tab_luzimp( i, j, k , 1 ) = b2 - a2 * tab_luzimp( i, j, k , 0 ); |
---|
385 | } |
---|
386 | a0=avt_fil*pdtz_fil/fse3w(i,j,NZ-2); |
---|
387 | a2 = a0/fse3t(i,j,NZ-2); |
---|
388 | b2 = 1.-a0/fse3t(i,j,NZ-2); |
---|
389 | #ifndef K_NEUM |
---|
390 | a1=avt_fil*pdtz_fil/fse3w(i,j,NZ-1); |
---|
391 | b2 -= a1/fse3t(i,j,NZ-2); |
---|
392 | #endif |
---|
393 | c2 = a0/fse3t(i,j,NZ-3); |
---|
394 | |
---|
395 | tab_luzimp( i, j, NZ-2 , 0 ) = c2/tab_luzimp( i, j, NZ-3 , 1); |
---|
396 | tab_luzimp( i, j, NZ-2 , 1 ) = b2 - a2 * tab_luzimp( i, j, NZ-2 , 0 ); |
---|
397 | } |
---|
398 | for(int k=0;k<NZ;k++){ // init a 0 |
---|
399 | printf("gam(%d)=%f \t",k+1,tab_luzimp( 23, 15, k , 0 )); |
---|
400 | printf("bet(%d)=%f \n",k+1,tab_luzimp( 23, 15, k , 1 )); |
---|
401 | } |
---|
402 | } |
---|
403 | #endif |
---|
404 | #endif |
---|
405 | //_____________________________________________________________________________ |
---|
406 | #define hu( i, j ) (hu[(j)*(NX)+(i)]) |
---|
407 | #define hv( i, j ) (hv[(j)*(NX)+(i)]) |
---|
408 | #define hur( i, j ) (hur[(j)*(NX)+(i)]) |
---|
409 | #define hvr( i, j ) (hvr[(j)*(NX)+(i)]) |
---|
410 | |
---|
411 | double hu[NY*NX]; // xdom_init |
---|
412 | double hv[NY*NX]; // xdom_init |
---|
413 | double hur[NY*NX]; // xdom_init |
---|
414 | double hvr[NY*NX]; // xdom_init |
---|
415 | |
---|
416 | void xdom_init(){ |
---|
417 | for(int j=0;j<NY;j++) |
---|
418 | for(int i=0;i<NX;i++){ |
---|
419 | hu(i,j) =0.; |
---|
420 | hv(i,j) =0.; |
---|
421 | for(int k=0;k<NZ;k++){ |
---|
422 | hu(i,j) = hu(i,j) + fse3u(i,j,k) * umask(i,j,k); |
---|
423 | hv(i,j) = hv(i,j) + fse3v(i,j,k) * vmask(i,j,k); |
---|
424 | } |
---|
425 | |
---|
426 | // ! Inverse of the local depth |
---|
427 | hur(i,j) = fse3u(i,j,0); //! Lower bound : thickness of the first model level |
---|
428 | hvr(i,j) = fse3v(i,j,0); |
---|
429 | for(int k=1;k<NZ;k++){ // ! Sum of the vertical scale factors |
---|
430 | hur(i,j) = hur(i,j) + fse3u(i,j,k) * umask(i,j,k); |
---|
431 | hvr(i,j) = hvr(i,j) + fse3v(i,j,k) * vmask(i,j,k); |
---|
432 | } |
---|
433 | // ! Cte and mask the inverse of the local depth |
---|
434 | hur(i,j) = 1. / hur(i,j) * umask(i,j,0); |
---|
435 | hvr(i,j) = 1. / hvr(i,j) * vmask(i,j,0); |
---|
436 | |
---|
437 | } |
---|
438 | } |
---|
439 | |
---|
440 | //_____________________________________________________________________________ |
---|
441 | #define gcp( i, j, k ,t ) (gcp[(j)*(NX)+(i)][(k)][(t)]) // a t0 |
---|
442 | #define gcdmat( i, j ,t ) (gcdmat[(j)*(NX)+(i)][(t)]) |
---|
443 | #define gcdprc( i, j ,t ) (gcdprc[(j)*(NX)+(i)][(t)]) |
---|
444 | #define gccd( i, j ,t ) (gccd[(j)*(NX)+(i)][(t)]) |
---|
445 | |
---|
446 | #define gcx( i, j ) (gcx[(i)][(j)]) |
---|
447 | double gcx[NX][NY]; // xsolsor |
---|
448 | #define G_gcx( i , j ) (G_gcx[(j)*(NX)+(i)]) |
---|
449 | double G_gcx[NY*NX]; // xsolsor |
---|
450 | #define G_gcx0( i , j ) (G_gcx0[(j)*(NX)+(i)]) |
---|
451 | double G_gcx0[NY*NX]; // xsolsor |
---|
452 | #define Ytb_gcx( i , j ) (Ytb_gcx[(j)*(NX)+(i)]) |
---|
453 | double Ytb_gcx[NY*NX]; // xsolsor |
---|
454 | #define gcr( i, j ) (gcr[(j)*(NX)+(i)]) |
---|
455 | double gcr[NY*NX]; // xsolsor |
---|
456 | |
---|
457 | double gcp[NY*NX][4][2]; // xsolmat_init |
---|
458 | double gcdmat[NY*NX][2]; // xsolmat_init |
---|
459 | double gcdprc[NY*NX][2]; // xsolmat_init |
---|
460 | double gccd[NY*NX][2]; // xsolmat_init |
---|
461 | |
---|
462 | int nn_rstssh=0; |
---|
463 | void xflt_rst(){ |
---|
464 | for(int j=0;j<NY;j++) |
---|
465 | for(int i=0;i<NX;i++){ |
---|
466 | YS_gcx_dynspg_flt(0,i,j,TU-1)=0.; YS_gcx_dynspg_flt(0,i,j,TU)=0.; |
---|
467 | YS_gcx2(0,i,j,TU-1)=0.; YS_gcx2(0,i,j,TU)=0.; |
---|
468 | gcx(i,j)=0.; |
---|
469 | } |
---|
470 | if( nn_rstssh == 1 ) { |
---|
471 | for(int j=0;j<NY;j++) |
---|
472 | for(int i=0;i<NX;i++){ |
---|
473 | YS_sshb(0,i,j,TU)=0.; YS_sshn(0,i,j,TU)=0.; |
---|
474 | } |
---|
475 | } |
---|
476 | } |
---|
477 | |
---|
478 | //_____________________________________________________________________________ |
---|
479 | void xsolmat_init(){ |
---|
480 | // ! initialize to zero |
---|
481 | double z2dt_t0,z2dt_t1;//zcoef_t0 = 0.e0,zcoef_t1 = 0.e0; |
---|
482 | for(int j=0;j<NY;j++) |
---|
483 | for(int i=0;i<NX;i++){// pour neuler=0 |
---|
484 | gcp(i,j,0,0) = 0.e0; gcp(i,j,0,1) = 0.e0; |
---|
485 | gcp(i,j,1,0) = 0.e0; gcp(i,j,1,1) = 0.e0; |
---|
486 | gcp(i,j,2,0) = 0.e0; gcp(i,j,2,1) = 0.e0; |
---|
487 | gcp(i,j,3,0) = 0.e0; gcp(i,j,3,1) = 0.e0; |
---|
488 | |
---|
489 | gcdprc(i,j,0) = 0.e0; gcdprc(i,j,1) = 0.e0; |
---|
490 | gcdmat(i,j,0) = 0.e0; gcdmat(i,j,1) = 0.e0; |
---|
491 | } |
---|
492 | z2dt_t0 = rdt; |
---|
493 | z2dt_t1 = 2.*rdt; |
---|
494 | double zcoef_t0,zcoefs_t0,zcoefw_t0,zcoefe_t0,zcoefn_t0; |
---|
495 | double zcoef_t1,zcoefs_t1,zcoefw_t1,zcoefe_t1,zcoefn_t1; |
---|
496 | // ! defined the coefficients for free surface elliptic system |
---|
497 | for(int j=1;j<NY-1;j++) |
---|
498 | for(int i=1;i<NX-1;i++){ |
---|
499 | zcoef_t0 = z2dt_t0 * z2dt_t0 * grav * rnu * bmask(i,j); |
---|
500 | zcoefs_t0 = -zcoef_t0 * hv(i ,j-1) * e1v(i ,j-1) / e2v(i ,j-1);// ! south coefficient |
---|
501 | zcoefw_t0 = -zcoef_t0 * hu(i-1,j ) * e2u(i-1,j ) / e1u(i-1,j );// ! west coefficient |
---|
502 | zcoefe_t0 = -zcoef_t0 * hu(i ,j ) * e2u(i ,j ) / e1u(i ,j );// ! east coefficient |
---|
503 | zcoefn_t0 = -zcoef_t0 * hv(i ,j ) * e1v(i ,j ) / e2v(i ,j );// ! north coefficient |
---|
504 | zcoef_t1 = z2dt_t1 * z2dt_t1 * grav * rnu * bmask(i,j); |
---|
505 | zcoefs_t1 = -zcoef_t1 * hv(i ,j-1) * e1v(i ,j-1) / e2v(i ,j-1);// ! south coefficient |
---|
506 | zcoefw_t1 = -zcoef_t1 * hu(i-1,j ) * e2u(i-1,j ) / e1u(i-1,j );// ! west coefficient |
---|
507 | zcoefe_t1 = -zcoef_t1 * hu(i ,j ) * e2u(i ,j ) / e1u(i ,j );// ! east coefficient |
---|
508 | zcoefn_t1 = -zcoef_t1 * hv(i ,j ) * e1v(i ,j ) / e2v(i ,j );// ! north coefficient |
---|
509 | |
---|
510 | gcp(i,j,0,0) = zcoefs_t0; gcp(i,j,0,1) = zcoefs_t1; |
---|
511 | gcp(i,j,1,0) = zcoefw_t0; gcp(i,j,1,1) = zcoefw_t1; |
---|
512 | gcp(i,j,2,0) = zcoefe_t0; gcp(i,j,2,1) = zcoefe_t1; |
---|
513 | gcp(i,j,3,0) = zcoefn_t0; gcp(i,j,3,1) = zcoefn_t1; |
---|
514 | gcdmat(i,j,0) = e1t(i,j) * e2t(i,j) * bmask(i,j) // ! diagonal coefficient |
---|
515 | - zcoefs_t0 -zcoefw_t0 -zcoefe_t0 -zcoefn_t0; |
---|
516 | gcdmat(i,j,1) = e1t(i,j) * e2t(i,j) * bmask(i,j) // ! diagonal coefficient |
---|
517 | - zcoefs_t1 -zcoefw_t1 -zcoefe_t1 -zcoefn_t1; |
---|
518 | |
---|
519 | } |
---|
520 | // ! SOR and PCG solvers |
---|
521 | for(int j=0;j<NY;j++) |
---|
522 | for(int i=0;i<NX;i++){ |
---|
523 | if( bmask(i,j) != 0 ){ |
---|
524 | gcdprc(i,j,0) = 1.e0 / gcdmat(i,j,0); |
---|
525 | gcdprc(i,j,1) = 1.e0 / gcdmat(i,j,1); |
---|
526 | } |
---|
527 | } |
---|
528 | |
---|
529 | for(int j=0;j<NY;j++) |
---|
530 | for(int i=0;i<NX;i++){ |
---|
531 | 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); |
---|
532 | 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); |
---|
533 | 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); |
---|
534 | 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); |
---|
535 | gccd(i,j,0) = sor * gcp(i,j,1,0); gccd(i,j,1) = sor * gcp(i,j,1,1); |
---|
536 | } |
---|
537 | } |
---|
538 | |
---|
539 | //_____________________________________________________________________________ |
---|
540 | void xistate_init(int argc, char *argv[]){ |
---|
541 | int itest=atoi(argv[1]); |
---|
542 | if(itest==1){// from restart file |
---|
543 | printf("start from rest_file\n "); |
---|
544 | neuler=0;//=1;// Restart from a file // Set time-step indicator at nit000 (leap-frog) |
---|
545 | // ifdef YE_*a_c alors *a=*b=*a_c |
---|
546 | char rest_file[200]; |
---|
547 | if (argc!=3) |
---|
548 | { |
---|
549 | printf("\nncf file for restart not specified in command xistate_init\n"); |
---|
550 | exit(1); |
---|
551 | } |
---|
552 | else |
---|
553 | { |
---|
554 | strcpy(rest_file,argv[2]); |
---|
555 | } |
---|
556 | //char const *rest_file = PATH_NCFILES"/rst3Yao.nc";//GYRE_00000003_restart.nc"; |
---|
557 | int rest_file_id; |
---|
558 | rest_file_id=Ouvre_nc(rest_file); |
---|
559 | xistate_rest_file(rest_file_id); |
---|
560 | nc_close( rest_file_id); |
---|
561 | //printf("kt=%f, ndastp=%f, adatrj=%f",kt,ndastp,adatrj); |
---|
562 | |
---|
563 | for(int i=0;i<NX;i++) |
---|
564 | for(int j=0;j<NY;j++){ |
---|
565 | for(int k=0;k<NZ;k++){ |
---|
566 | YS_ua_c(0,i,j,k)=YS_ua(0,i,j,k,TU); |
---|
567 | YS_va_c(0,i,j,k)=YS_va(0,i,j,k,TU); |
---|
568 | YS_ta_c(0,i,j,k)=YS_ta(0,i,j,k,TU); |
---|
569 | YS_sa_c(0,i,j,k)=YS_sa(0,i,j,k,TU); |
---|
570 | |
---|
571 | tb_neuler1(i,j,k)=YS_tb(0,i,j,k,TU); |
---|
572 | sb_neuler1(i,j,k)=YS_sb(0,i,j,k,TU); |
---|
573 | ub_neuler1(i,j,k)=YS_ub(0,i,j,k,TU); |
---|
574 | vb_neuler1(i,j,k)=YS_vb(0,i,j,k,TU); |
---|
575 | rotn_at_TU( i, j, k)=YS_rotn(0,i,j,k,TU); |
---|
576 | hdivn_at_TU( i, j, k)=YS_hdivn(0,i,j,k,TU); |
---|
577 | } |
---|
578 | YS_sshn_c(0,i,j)=YS_sshn(0,i,j,TU); |
---|
579 | sshb_neuler1(i,j)=YS_sshb(0,i,j,TU); |
---|
580 | gcx_at_TU( i, j)=YS_gcx2(0,i,j,TU); |
---|
581 | } |
---|
582 | } |
---|
583 | else if (itest==2){ |
---|
584 | neuler=0;// |
---|
585 | xistate_gyre(); // GYRE configuration : start from pre-defined temperature |
---|
586 | // and salinity fields |
---|
587 | } |
---|
588 | else if (itest==3){ |
---|
589 | neuler=0;// |
---|
590 | xistate_yao_file("data_in/file_yao/"); // GYRE configuration : start from pre-defined temperature |
---|
591 | // and salinity fields |
---|
592 | } |
---|
593 | } |
---|
594 | |
---|
595 | //_____________________________________________________________________________ |
---|
596 | void xeos_init(){ |
---|
597 | } |
---|
598 | |
---|
599 | //_____________________________________________________________________________ |
---|
600 | void xistate_gyre(){ |
---|
601 | printf("start from gyre analytic initialization\n"); |
---|
602 | for(int i=0;i<NX;i++){ |
---|
603 | for(int j=0;j<NY;j++){ |
---|
604 | YS_sshn(0,i,j,TU)=0.; |
---|
605 | YS_sshb(0,i,j,TU)=0.; |
---|
606 | YS_gcx2(0,i,j,TU)=0.; |
---|
607 | YS_gcx2(0,i,j,TU-1)=0.; |
---|
608 | for(int k=0;k<NZ;k++){ |
---|
609 | YS_ua_c(0,i,j,k)=0.; |
---|
610 | YS_ua(0,i,j,k,TU)=0.; |
---|
611 | YS_ub(0,i,j,k,TU)=0.; |
---|
612 | YS_va_c(0,i,j,k)=0.; |
---|
613 | YS_va(0,i,j,k,TU)=0.; |
---|
614 | YS_vb(0,i,j,k,TU)=0.; |
---|
615 | YS_rotn(0,i,j,k,TU)=0.; |
---|
616 | YS_rotn(0,i,j,k,TU-1)=0.; |
---|
617 | YS_hdivn(0,i,j,k,TU)=0.; |
---|
618 | YS_hdivn(0,i,j,k,TU-1)=0.; |
---|
619 | YS_ta_c(0,i,j,k) = (16.-12.*tanh((fsdept(i,j,k)-400.)/700.)) |
---|
620 | * (-tanh((500.-fsdept(i,j,k))/150.)+1.)/2. |
---|
621 | + (15.*(1.-tanh((fsdept(i,j,k)-50.)/1500.))-1.4*tanh((fsdept(i,j,k)-100.)/100.) |
---|
622 | + 7.*(1500.-fsdept(i,j,k))/1500.) |
---|
623 | *(-tanh((fsdept(i,j,k)-500.)/150.)+1.)/2.; |
---|
624 | YS_ta_c(0,i,j,k) = YS_ta_c(0,i,j,k)* tmask(i,j,k); |
---|
625 | YS_ta(0,i,j,k,TU) = YS_ta_c(0,i,j,k); |
---|
626 | YS_tb(0,i,j,k,TU) = YS_ta_c(0,i,j,k); |
---|
627 | |
---|
628 | YS_sa_c(0,i,j,k) = (36.25-1.13*tanh((fsdept(i,j,k)-305.)/460.)) |
---|
629 | *(-tanh((500.-fsdept(i,j,k))/150.)+1.)/2. |
---|
630 | +(35.55+1.25*(5000.-fsdept(i,j,k))/5000. |
---|
631 | -1.62*tanh((fsdept(i,j,k)-60.)/650.) |
---|
632 | +0.2*tanh((fsdept(i,j,k)-35.)/100.) |
---|
633 | +0.2*tanh((fsdept(i,j,k)-1000.)/5000.)) |
---|
634 | *(-tanh((fsdept(i,j,k)-500.)/150.)+1.)/2.; |
---|
635 | YS_sa_c(0,i,j,k) = YS_sa_c(0,i,j,k)*tmask(i,j,k); |
---|
636 | YS_sa(0,i,j,k,TU) = YS_sa_c(0,i,j,k); |
---|
637 | YS_sb(0,i,j,k,TU) = YS_sa_c(0,i,j,k); |
---|
638 | |
---|
639 | rotn_at_TU( i, j, k)=YS_rotn(0,i,j,k,TU); |
---|
640 | hdivn_at_TU( i, j, k)=YS_hdivn(0,i,j,k,TU); |
---|
641 | |
---|
642 | } |
---|
643 | YS_sshn_c(0,i,j)=YS_sshn(0,i,j,TU); |
---|
644 | gcx_at_TU( i, j)=YS_gcx2(0,i,j,TU); |
---|
645 | } |
---|
646 | } |
---|
647 | } |
---|
648 | |
---|
649 | void xchangesavemode(int argc, char *argv[]){ |
---|
650 | short newmode=atoi(argv[1]); |
---|
651 | savemode=newmode; |
---|
652 | if (newmode==2) { |
---|
653 | isave=atoi(argv[2]); |
---|
654 | jsave=atoi(argv[3]); |
---|
655 | ksave=atoi(argv[4]); |
---|
656 | printf("\n savemode : 1 point (Yi=%d,Yj=%d,Yk=%d)\n",isave,jsave,ksave); |
---|
657 | } |
---|
658 | |
---|
659 | } |
---|
660 | |
---|
661 | //_____________________________________________________________________________ |
---|
662 | void xistate_yao_file(char * pth_file_yao){ // a faire |
---|
663 | /* printf("start from yao file\n"); |
---|
664 | FILE *pf; |
---|
665 | char *file_yao; |
---|
666 | sprintf(file_yao,"%stb.dat",pth_file_yao); |
---|
667 | pf=fopen(file_yao,'r'); |
---|
668 | fscanf(p,"%lf",); |
---|
669 | fclose(p); |
---|
670 | for(int j=0;j<NY;j++) |
---|
671 | for(int i=0;i<NX;i++){ |
---|
672 | for(int k=0;k<NZ;k++){ |
---|
673 | YS_ua_c(0,i,j,k)=0.; |
---|
674 | YS_va_c(0,i,j,k)=0.; |
---|
675 | YS_ta_c(0,i,j,k) =; |
---|
676 | YS_sa_c(0,i,j,k) =; |
---|
677 | } |
---|
678 | YS_sshn_c(0,i,j)=; |
---|
679 | }*/ |
---|
680 | } |
---|
681 | //_____________________________________________________________________________ |
---|
682 | void xtraqsr_init(){// init of solar radiation penetration |
---|
683 | //! z-coordinate with or without partial step : same w-level everywhere inside the ocean |
---|
684 | for(int k=0;k<NZ;k++) gdsr[k] = 0.e0; |
---|
685 | for(int k=0;k<NZ;k++) { |
---|
686 | double zdp1 = -gdepw_0[k]; //#define fsdepw( i , j , k ) gdepw_0[k] |
---|
687 | gdsr[k] = ro0cpr*(rabs*exp(zdp1/xsi1)+(1.-rabs)*exp(zdp1/xsi2)); |
---|
688 | if(gdsr[k]<= 1.e-10) break; |
---|
689 | } |
---|
690 | int indic = 0; |
---|
691 | for(int k=0;k<NZ;k++) { |
---|
692 | if( gdsr[k]<=1.e-15 && indic == 0){ |
---|
693 | gdsr[k] = 0.e0; |
---|
694 | nksr = k+1; |
---|
695 | indic = 1; |
---|
696 | } |
---|
697 | if(nksr>NZ-1) nksr=NZ-1; |
---|
698 | } |
---|
699 | } |
---|
700 | |
---|
701 | //_____________________________________________________________________________ |
---|
702 | void cal_ff() //domhgr.F90 Coriolis parameter |
---|
703 | { |
---|
704 | if(jphgr_msh==5){// beta-plane and rotated domain (gyre configuration) |
---|
705 | double zbeta = 2.*omega*cos(rad*ppgphi0)/ra ;// beta at latitude ppgphi0 |
---|
706 | double zphi0 = 15.;//latitude of the first row F-points |
---|
707 | double zf0 = 2.*omega*sin(rad*zphi0);// compute f0 1st point south |
---|
708 | for(int i=0;i<NX;i++) |
---|
709 | for(int j=0;j<NY;j++) |
---|
710 | ff(i,j) = (zf0+zbeta*fabs(gphif(i,j)-zphi0)*rad*ra);// f = f0 +beta* y ( y=0 at south) |
---|
711 | } |
---|
712 | } |
---|
713 | //______________________________________________________________________________ |
---|
714 | #ifdef K_FILTER |
---|
715 | void load_eb (int argc, char *argv[]){ // load "ebauche" |
---|
716 | // on doit presiser la variable exple : load_eb ta_c |
---|
717 | if (!strcmp(argv[1],"ta_c")){ |
---|
718 | for(int i=0;i<NX;i++) |
---|
719 | for(int j=0;j<NY;j++) |
---|
720 | for(int k=0;k<NZ;k++){ |
---|
721 | ta_eb( i, j, k )=0.0; |
---|
722 | } |
---|
723 | for(int i=0;i<NX;i++) |
---|
724 | for(int j=0;j<NY;j++) |
---|
725 | for(int k=0;k<NZ;k++){ |
---|
726 | ta_eb( i, j, k )=YS_ta_c(0,i,j,k); |
---|
727 | } |
---|
728 | |
---|
729 | } |
---|
730 | else{ |
---|
731 | printf( "function load_eb:: argument incorrect\n"); |
---|
732 | exit(99);// abort program |
---|
733 | } |
---|
734 | } |
---|
735 | #endif |
---|
736 | |
---|
737 | |
---|
738 | //_____________________________________________________________________________ |
---|
739 | /*void load_shape_func (int argc, char *argv[]) // load PCA axes from file |
---|
740 | { |
---|
741 | FILE *ffp; |
---|
742 | if ((ffp=fopen( argv[2], "r")) == NULL){ |
---|
743 | printf( "error: could not find file of shape functions '%s'\n", argv[2]); |
---|
744 | exit(99);// abort program |
---|
745 | } |
---|
746 | if (!strcmp(argv[1],"ta_c")){ |
---|
747 | for(int i=0;i<NX;i++) |
---|
748 | for(int j=0;j<NY;j++) |
---|
749 | for(int k=0;k<NZ;k++){ |
---|
750 | for(int n=0;n<NPCA;n++) shfs_ta( i, j, k ,n )=0.0; |
---|
751 | } |
---|
752 | |
---|
753 | int ctrl=-1; |
---|
754 | while (ctrl<NZ*NY*NX-1){//!feof(ffp)) { |
---|
755 | ctrl=ctrl+1; |
---|
756 | for(int n=0;n<NPCA;n++) fscanf(ffp,"%lf",&shfs_ta[ctrl][n]); |
---|
757 | } |
---|
758 | } |
---|
759 | else{ |
---|
760 | printf( "load_shape_func:: pas encore fait\n"); |
---|
761 | exit(99);// abort program |
---|
762 | } |
---|
763 | |
---|
764 | fclose(ffp); |
---|
765 | |
---|
766 | |
---|
767 | //////////////// |
---|
768 | }*/ |
---|
769 | |
---|
770 | //______________________________________________________________________________ |
---|
771 | /*void load_stdev_pca (int argc, char *argv[]) // load stdev from file |
---|
772 | { |
---|
773 | FILE *ffp; |
---|
774 | if ((ffp=fopen( argv[2], "r")) == NULL){ |
---|
775 | printf( "error: could not find file of stdev '%s'\n", argv[1]); |
---|
776 | exit(99);// abort program |
---|
777 | } |
---|
778 | if (!strcmp(argv[1],"ta_c")){ |
---|
779 | for(int n=0;n<NPCA;n++) fscanf(ffp,"%lf",&stdev_ta[n]); |
---|
780 | } |
---|
781 | else{ |
---|
782 | printf( "load_shape_func:: pas encore fait\n"); |
---|
783 | exit(99);// abort program |
---|
784 | } |
---|
785 | fclose(ffp); |
---|
786 | |
---|
787 | |
---|
788 | //////////////// |
---|
789 | }*/ |
---|
790 | //______________________________________________________________________________ |
---|
791 | /*void load_mean (int argc, char *argv[]){ // load mean to generate obs |
---|
792 | if (!strcmp(argv[2],"ta_c")){ |
---|
793 | if(atoi(argv[1])==0){ |
---|
794 | for(int i=0;i<NX;i++) |
---|
795 | for(int j=0;j<NY;j++) |
---|
796 | for(int k=0;k<NZ;k++){ |
---|
797 | moy_ta( i, j, k )=0.0; |
---|
798 | } |
---|
799 | for(int i=0;i<NX;i++) |
---|
800 | for(int j=0;j<NY;j++) |
---|
801 | for(int k=0;k<NZ;k++){ |
---|
802 | moy_ta( i, j, k )=YS_ta_c(0,i,j,k); |
---|
803 | } |
---|
804 | } |
---|
805 | else if(atoi(argv[1])==1){ |
---|
806 | FILE *ffp; |
---|
807 | if ((ffp=fopen( argv[3], "r")) == NULL){ |
---|
808 | printf( "error: could not find file of shape functions '%s'\n", argv[3]); |
---|
809 | exit(99);// abort program |
---|
810 | } |
---|
811 | for(int i=0;i<NX;i++) |
---|
812 | for(int j=0;j<NY;j++) |
---|
813 | for(int k=0;k<NZ;k++){ |
---|
814 | moy_ta( i, j, k )=0.0; |
---|
815 | } |
---|
816 | |
---|
817 | int ctrl=-1; |
---|
818 | while (ctrl<NZ*NY*NX-1){//!feof(ffp)) { |
---|
819 | ctrl=ctrl+1; |
---|
820 | fscanf(ffp,"%lf",&moy_ta[ctrl]); |
---|
821 | } |
---|
822 | } |
---|
823 | } |
---|
824 | else{ |
---|
825 | printf( "load_shape_func:: pas encore fait\n"); |
---|
826 | exit(99);// abort program |
---|
827 | } |
---|
828 | }*/ |
---|
829 | //______________________________________________________________________________ |
---|
830 | |
---|
831 | //---------------------------------------------- |
---|
832 | #if defined (YE_ta_c) || defined (YE_pca_ta) || defined (YE_nu) |
---|
833 | #define true_tac( i, j, k ) (true_tac[(k)*(NY*NX)+(j)*(NX)+(i)]) |
---|
834 | double true_tac[NZ*NY*NX]; |
---|
835 | #endif |
---|
836 | |
---|
837 | #ifdef YE_sa_c |
---|
838 | #define true_sac( i, j, k ) (true_sac[(k)*(NY*NX)+(j)*(NX)+(i)]) |
---|
839 | double true_sac[NZ*NY*NX]; |
---|
840 | #endif |
---|
841 | |
---|
842 | #ifdef YE_ua_c |
---|
843 | #define true_uac( i, j, k ) (true_uac[(k)*(NY*NX)+(j)*(NX)+(i)]) |
---|
844 | double true_uac[NZ*NY*NX]; |
---|
845 | #endif |
---|
846 | |
---|
847 | #ifdef YE_va_c |
---|
848 | #define true_vac( i, j, k ) (true_vac[(k)*(NY*NX)+(j)*(NX)+(i)]) |
---|
849 | double true_vac[NZ*NY*NX]; |
---|
850 | #endif |
---|
851 | //---------------------------------------------- |
---|
852 | void true_target_in_tab(int argc, char *argv[]){ |
---|
853 | int i,j,k; |
---|
854 | if(!strcmp(argv[1],"ta_c")){ |
---|
855 | #if defined (YE_ta_c) || defined (YE_pca_ta) || defined (YE_nu) |
---|
856 | for(i=0;i<NX;i++) |
---|
857 | for(j=0;j<NY;j++) |
---|
858 | for(k=0;k<NZ;k++) |
---|
859 | true_tac( i, j , k ) =YS_ta_c(0,i,j,k); |
---|
860 | #endif |
---|
861 | } |
---|
862 | |
---|
863 | if(!strcmp(argv[1],"sa_c")){ |
---|
864 | #ifdef YE_sa_c |
---|
865 | for(i=0;i<NX;i++) |
---|
866 | for(j=0;j<NY;j++) |
---|
867 | for(k=0;k<NZ;k++) |
---|
868 | true_sac( i, j , k ) =YS_sa_c(0,i,j,k); |
---|
869 | #endif |
---|
870 | } |
---|
871 | |
---|
872 | if(!strcmp(argv[1],"ua_c")){ |
---|
873 | #ifdef YE_ua_c |
---|
874 | for(i=0;i<NX;i++) |
---|
875 | for(j=0;j<NY;j++) |
---|
876 | for(k=0;k<NZ;k++) |
---|
877 | true_uac( i, j , k ) =YS_ua_c(0,i,j,k); |
---|
878 | #endif |
---|
879 | } |
---|
880 | |
---|
881 | if(!strcmp(argv[1],"va_c")){ |
---|
882 | #ifdef YE_va_c |
---|
883 | for(i=0;i<NX;i++) |
---|
884 | for(j=0;j<NY;j++) |
---|
885 | for(k=0;k<NZ;k++) |
---|
886 | true_vac( i, j , k ) =YS_va_c(0,i,j,k); |
---|
887 | #endif |
---|
888 | } |
---|
889 | |
---|
890 | } |
---|
891 | int iii=28,jjj=12; |
---|
892 | //--------------------------- |
---|
893 | void normsup_alla_c(){ |
---|
894 | double sum=0.,sumn=0.; |
---|
895 | int i,j,k; |
---|
896 | |
---|
897 | #if defined (YE_ta_c) || defined (YE_pca_ta) || defined (YE_nu) |
---|
898 | sum=0.;sumn=0.; |
---|
899 | for(i=0;i<NX;i++) |
---|
900 | for(j=0;j<NY;j++) |
---|
901 | // if(i==iii && j==jjj){ |
---|
902 | for(k=0;k<NZ;k++){ |
---|
903 | 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)); |
---|
904 | 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)); |
---|
905 | }//} |
---|
906 | sumn=sqrt(sumn); |
---|
907 | printf("norm_sup_diff_ta_c : %23.16e\n",sum); |
---|
908 | printf("norm_2_diff_ta_c : %23.16e\n",sumn); |
---|
909 | #endif |
---|
910 | |
---|
911 | #ifdef YE_sa_c |
---|
912 | sum=0.;sumn=0.; |
---|
913 | for(i=0;i<NX;i++) |
---|
914 | for(j=0;j<NY;j++) |
---|
915 | //if(i==iii && j==jjj){ |
---|
916 | for(k=0;k<NZ;k++){ |
---|
917 | 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)); |
---|
918 | }//} |
---|
919 | printf("norm_sup_diff_sa_c : %23.16e\n",sum); |
---|
920 | #endif |
---|
921 | |
---|
922 | #ifdef YE_ua_c |
---|
923 | sum=0.;sumn=0.; |
---|
924 | for(i=0;i<NX;i++) |
---|
925 | for(j=0;j<NY;j++) |
---|
926 | for(k=0;k<NZ;k++){ |
---|
927 | 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)); |
---|
928 | } |
---|
929 | printf("norm_sup_diff_ua_c : %23.16e\n",sum); |
---|
930 | #endif |
---|
931 | |
---|
932 | #ifdef YE_va_c |
---|
933 | sum=0.;sumn=0.; |
---|
934 | for(i=0;i<NX;i++) |
---|
935 | for(j=0;j<NY;j++) |
---|
936 | for(k=0;k<NZ;k++){ |
---|
937 | 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)); |
---|
938 | } |
---|
939 | printf("norm_sup_diff_va_c : %23.16e\n",sum); |
---|
940 | #endif |
---|
941 | |
---|
942 | } |
---|
943 | |
---|
944 | //__________________________________________________________________ |
---|
945 | void xcomparYF(); |
---|
946 | void xcomparYF(){ |
---|
947 | //comparer Fortran vs Yao |
---|
948 | int i,j,k,t,ks=1,ite=0; |
---|
949 | double v,sum; |
---|
950 | FILE *pf; |
---|
951 | char* fua="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_pua.dat"; |
---|
952 | if((pf=fopen(fua,"r"))==NULL){ |
---|
953 | printf("file in function xcomparYF\n"); |
---|
954 | exit(99); |
---|
955 | } |
---|
956 | sum=0.; |
---|
957 | while(!feof(pf)){ |
---|
958 | fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v); |
---|
959 | if(t==TA-3 && k==ks){ |
---|
960 | 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)); |
---|
961 | ite=1; |
---|
962 | } |
---|
963 | } |
---|
964 | if(ite==0) |
---|
965 | printf("xcomparYF: problem dans ua\n"); |
---|
966 | else |
---|
967 | printf("norm_sup(uaF-uaY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum); |
---|
968 | fclose(pf); |
---|
969 | |
---|
970 | char* fva="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_pva.dat"; |
---|
971 | if((pf=fopen(fva,"r"))==NULL){ |
---|
972 | printf("file in function xcomparYF\n"); |
---|
973 | exit(99); |
---|
974 | } |
---|
975 | |
---|
976 | ite=0; |
---|
977 | sum=0.; |
---|
978 | while(feof(pf)==0){ |
---|
979 | fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v); |
---|
980 | if(t==TA-3 && k==ks){ |
---|
981 | 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)); |
---|
982 | ite=1; |
---|
983 | } |
---|
984 | } |
---|
985 | if(ite==0) |
---|
986 | printf("xcomparYF: problem dans va\n"); |
---|
987 | else |
---|
988 | printf("norm_sup(vaF-vaY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum); |
---|
989 | fclose(pf); |
---|
990 | |
---|
991 | char* fwa="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_pwn.dat"; |
---|
992 | if((pf=fopen(fwa,"r"))==NULL){ |
---|
993 | printf("file in function xcomparYF\n"); |
---|
994 | exit(99); |
---|
995 | } |
---|
996 | |
---|
997 | ite=0; |
---|
998 | sum=0.; |
---|
999 | while(!feof(pf)){ |
---|
1000 | fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v); |
---|
1001 | if(t==TA-3 && k==ks){ |
---|
1002 | 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)); |
---|
1003 | ite=1; |
---|
1004 | } |
---|
1005 | } |
---|
1006 | if(ite==0) |
---|
1007 | printf("xcomparYF: problem dans wa\n"); |
---|
1008 | else |
---|
1009 | printf("norm_sup(waF-waY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum); |
---|
1010 | fclose(pf); |
---|
1011 | |
---|
1012 | char* fta="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_pta.dat"; |
---|
1013 | if((pf=fopen(fta,"r"))==NULL){ |
---|
1014 | printf("file in function xcomparYF\n"); |
---|
1015 | exit(99); |
---|
1016 | } |
---|
1017 | |
---|
1018 | ite=0; |
---|
1019 | sum=0.; |
---|
1020 | while(!feof(pf)){ |
---|
1021 | fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v); |
---|
1022 | if(t==TA-3 && k==ks){ |
---|
1023 | 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)); |
---|
1024 | ite=1; |
---|
1025 | } |
---|
1026 | } |
---|
1027 | if(ite==0) |
---|
1028 | printf("xcomparYF: problem dans ta\n"); |
---|
1029 | else |
---|
1030 | printf("norm_sup(taF-taY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum); |
---|
1031 | fclose(pf); |
---|
1032 | |
---|
1033 | char* fsa="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_psa.dat"; |
---|
1034 | if((pf=fopen(fsa,"r"))==NULL){ |
---|
1035 | printf("file in function xcomparYF\n"); |
---|
1036 | exit(99); |
---|
1037 | } |
---|
1038 | |
---|
1039 | ite=0; |
---|
1040 | sum=0.; |
---|
1041 | while(!feof(pf)){ |
---|
1042 | fscanf(pf," %i %i %i %i %lf",&t,&i,&j,&k,&v); |
---|
1043 | if(t==TA-3 && k==ks){ |
---|
1044 | 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)); |
---|
1045 | ite=1; |
---|
1046 | } |
---|
1047 | } |
---|
1048 | if(ite==0) |
---|
1049 | printf("xcomparYF: problem dans sa\n"); |
---|
1050 | else |
---|
1051 | printf("norm_sup(saF-saY)(t=%i,k=%i) : %23.16e\n",TA,ks,sum); |
---|
1052 | fclose(pf); |
---|
1053 | |
---|
1054 | |
---|
1055 | char* fsshn="data_in/comparF_in/traadv_cen2/GYRE_sortie_mohamed_psshn.dat"; |
---|
1056 | if((pf=fopen(fsshn,"r"))==NULL){ |
---|
1057 | printf("file in function xcomparYF\n"); |
---|
1058 | exit(99); |
---|
1059 | } |
---|
1060 | |
---|
1061 | ite=0; |
---|
1062 | sum=0.; |
---|
1063 | while(!feof(pf)){ |
---|
1064 | fscanf(pf," %i %i %i %lf",&t,&i,&j,&v); |
---|
1065 | if(t==TA-3){ |
---|
1066 | 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)); |
---|
1067 | ite=1; |
---|
1068 | } |
---|
1069 | } |
---|
1070 | if(ite==0) |
---|
1071 | printf("xcomparYF: problem dans sshn\n"); |
---|
1072 | else |
---|
1073 | printf("norm_sup(sshnF-sshnY)(t=%i) : %23.16e\n",TA,sum); |
---|
1074 | fclose(pf); |
---|
1075 | |
---|
1076 | } |
---|
1077 | |
---|
1078 | //Gestion des fichiers netcdf |
---|
1079 | |
---|
1080 | void xrst_save(int argc, char *argv[]) { |
---|
1081 | //Sauvegarde un restart |
---|
1082 | if (argc!=2) |
---|
1083 | { |
---|
1084 | printf("\ncf file for saving restart not specified in command xrst_save\n"); |
---|
1085 | exit(1); |
---|
1086 | } |
---|
1087 | else |
---|
1088 | { |
---|
1089 | char rest_file[200]; |
---|
1090 | int rest_file_id; |
---|
1091 | int dimid[4]; |
---|
1092 | strcpy(rest_file,argv[1]); |
---|
1093 | rest_file_id=Ouvre_nc_write(rest_file); //Ouvre le fichier netcds pour l'écriture |
---|
1094 | write_rst_global_att(rest_file_id); //Ecrit les global attributes |
---|
1095 | define_dim(rest_file_id,dimid); //Définit les dimensions |
---|
1096 | printf("\n Save restart at kt=%f (Yt=%d)\n",kt,Yt); |
---|
1097 | write_rst_var(kt-nit000+TU,rest_file_id,dimid); //Ecrit les variables |
---|
1098 | nc_close( rest_file_id); //ferme le fichier |
---|
1099 | } |
---|
1100 | } |
---|
1101 | |
---|
1102 | void xwriteout (int argc, char *argv[]) { |
---|
1103 | //xwriteout t suff ncfile |
---|
1104 | //Ecrit dans le fichier ncfile les variables u,v,t,s,ssh avec le suffixe suff au pas de temps t |
---|
1105 | if (argc!=4) |
---|
1106 | { |
---|
1107 | printf("\n not enough arguments in writeout\n"); |
---|
1108 | exit(1); |
---|
1109 | } |
---|
1110 | |
---|
1111 | int t; |
---|
1112 | char suff[200]; |
---|
1113 | char ncfile[200]; |
---|
1114 | int ncfile_id; |
---|
1115 | t=atoi(argv[1]); |
---|
1116 | strcpy(suff,argv[2]); |
---|
1117 | strcpy(ncfile,argv[3]); |
---|
1118 | ncfile_id=Ouvre_nc_add(ncfile); |
---|
1119 | write_out_var(t,ncfile_id,suff); |
---|
1120 | nc_close(ncfile_id); |
---|
1121 | |
---|
1122 | } |
---|
1123 | void xwritegrad (int argc, char *argv[]) { |
---|
1124 | //xwritegrad t suff ncfile |
---|
1125 | //Ecrit dans le fichier ncfile les grad des variables u,v,t,s,ssh avec le suffixe suff au pas de temps t |
---|
1126 | if (argc!=4) |
---|
1127 | { |
---|
1128 | printf("\n not enough arguments in writeout\n"); |
---|
1129 | exit(1); |
---|
1130 | } |
---|
1131 | |
---|
1132 | int t; |
---|
1133 | char suff[200]; |
---|
1134 | char ncfile[200]; |
---|
1135 | int ncfile_id; |
---|
1136 | t=atoi(argv[1]); |
---|
1137 | strcpy(suff,argv[2]); |
---|
1138 | strcpy(ncfile,argv[3]); |
---|
1139 | ncfile_id=Ouvre_nc_add(ncfile); |
---|
1140 | write_grad_var(t,ncfile_id,suff); |
---|
1141 | nc_close(ncfile_id); |
---|
1142 | |
---|
1143 | } |
---|
1144 | |
---|
1145 | void xinitnc(int argc,char *argv[]) { |
---|
1146 | //xinit ncfile |
---|
1147 | //init le fichier netcdf pour sauvegarde |
---|
1148 | if (argc!=2) |
---|
1149 | { |
---|
1150 | printf("\ncf file for saving output not specified in command xinitnc\n"); |
---|
1151 | exit(1); |
---|
1152 | } |
---|
1153 | else |
---|
1154 | { |
---|
1155 | char ncfile[200]; |
---|
1156 | int nc_file_id; |
---|
1157 | int dimid[4]; |
---|
1158 | strcpy(ncfile,argv[1]); |
---|
1159 | nc_file_id=Ouvre_nc_write(ncfile); |
---|
1160 | define_dim(nc_file_id,dimid); |
---|
1161 | write_out_var_init(nc_file_id,dimid); |
---|
1162 | nc_close(nc_file_id); |
---|
1163 | } |
---|
1164 | |
---|
1165 | } |
---|
1166 | |
---|
1167 | void init_euler(int argc, char *argv[]){ |
---|
1168 | if (argc!=2) |
---|
1169 | { |
---|
1170 | printf("\nwrong number of argument in init_euler : specify 0 or 1\n"); |
---|
1171 | exit(1); |
---|
1172 | } |
---|
1173 | else |
---|
1174 | { |
---|
1175 | int n; |
---|
1176 | n=atoi(argv[1]); |
---|
1177 | if(n!=0 && n!=1) |
---|
1178 | { |
---|
1179 | printf("\nwrong argument in init_euler : specify 0 or 1 in command init_euler\n"); |
---|
1180 | exit(1); |
---|
1181 | } |
---|
1182 | neuler=n; |
---|
1183 | } |
---|
1184 | } |
---|
1185 | |
---|
1186 | void init_kt(int argc, char *argv[]){ |
---|
1187 | if (argc!=2) |
---|
1188 | { |
---|
1189 | printf("\nwrong number of argument in init_kt : specify a positive number\n"); |
---|
1190 | exit(1); |
---|
1191 | } |
---|
1192 | else |
---|
1193 | { |
---|
1194 | double t0; |
---|
1195 | t0=atof(argv[1]); |
---|
1196 | if(t0<0) |
---|
1197 | { |
---|
1198 | printf("\nwrong argument in init_kt : specify a postiive number\n"); |
---|
1199 | exit(1); |
---|
1200 | } |
---|
1201 | nit000=t0; |
---|
1202 | kt=nit000; |
---|
1203 | } |
---|
1204 | } |
---|
1205 | |
---|
1206 | /* Plus tard pour tester l'erreur modÚle. |
---|
1207 | void genere_obs(int argc, char *argv[]){ |
---|
1208 | //genere_obs Module increment t0 dt tend i j k |
---|
1209 | if (argc!=9) |
---|
1210 | { |
---|
1211 | printf("\n wrong number of argument in genere_obs \n"); |
---|
1212 | exit(1); |
---|
1213 | } |
---|
1214 | else |
---|
1215 | { |
---|
1216 | |
---|
1217 | } |
---|
1218 | |
---|
1219 | |
---|
1220 | |
---|
1221 | } |
---|
1222 | */ |
---|