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