/////////////////////////////////////////////////////////////////////////////// // define, globales, fonctions perso (obligatoire et autres) ... // mohamed.berrada@upmc.fr /////////////////////////////////////////////////////////////////////////////// #include"../include/ncutil.h"//netcdf #define key_zco #include"../include/phy_cst_file.h" #include"../include/namelist.h" //#define PATH_NCFILES "/usr/home/mblod/neuro/nemo/nemo_opa/TRY/modipsl/config/GYRE/EXP00" #define PATH_NCFILES "../data_in" #include"../include/meshmask.h" #ifdef K_FILTER #ifdef K_FZIMP void luzimp_init(); // initialise tab_luzimp #define tab_luzimp( i, j, k , t ) (tab_luzimp[(i)*(NY*NZ)+(j)*(NZ)+(k)][t]) // factorisation LU double tab_luzimp[NZ*NY*NX][2]; // tableau factorisation LU #endif #define ta_eb( i, j, k ) (ta_eb[(i)*(NY*NZ)+(j)*(NZ)+(k)]) // ebauche double ta_eb[NZ*NY*NX]; // ebauche void load_eb (int argc, char *argv[]); // load "ebauche" from current variable ta_c #endif // For PPCA /*double shfs_ta[NZ*NY*NX][NPCA]; double moy_ta[NZ*NY*NX]; double stdev_ta[NPCA]; #define shfs_ta( i, j, k ,n ) (shfs_ta[(i)*(NY*NZ)+(j)*(NZ)+(k)][(n)]) // shape functions #define moy_ta( i, j, k ) (moy_ta[(i)*(NY*NZ)+(j)*(NZ)+(k)]) // mean void load_shape_func (int argc, char *argv[]); // load PCA axes from file with mean void load_stdev_pac (int argc, char *argv[]); // load stdev from file void load_mean (int argc, char *argv[]); // load mean */ ////#define bmask( i, j ) (tmask( i, j ,0)) // ! elliptic equation is written at t-point int tniter[TA+TU]; // nbre d'iterations pour backward de solsor_dynspg_flt.h int ndiag_rst=0; // diagnostic restart if =1 //Arrays for neuler=1 (neuler=0::Euler 1st time step) #define tb_neuler1( i, j, k) (tb_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)]) #define sb_neuler1( i, j, k) (sb_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)]) #define ub_neuler1( i, j, k) (ub_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)]) #define vb_neuler1( i, j, k) (vb_neuler1[(k)*(NY*NX)+(j)*(NX)+(i)]) #define sshb_neuler1( i, j) (sshb_neuler1[(j)*(NX)+(i)]) #define rotn_at_TU( i, j, k) (rotn_at_TU[(k)*(NY*NX)+(j)*(NX)+(i)]) #define hdivn_at_TU( i, j, k) (hdivn_at_TU[(k)*(NY*NX)+(j)*(NX)+(i)]) #define gcx_at_TU( i, j) (gcx_at_TU[(j)*(NX)+(i)]) #define masque_obs( i, j) (masque_obs[(j)*(NX)+(i)]) double tb_neuler1[NZ*NY*NX]; // double sb_neuler1[NZ*NY*NX]; // double ub_neuler1[NZ*NY*NX]; // double vb_neuler1[NZ*NY*NX]; // double sshb_neuler1[NY*NX]; // double rotn_at_TU[NZ*NY*NX]; // double hdivn_at_TU[NZ*NY*NX]; // double gcx_at_TU[NY*NX]; // short masque_obs[NY*NX]; // void true_target_in_tab(int argc, char *argv[]); // charge YS_*a_c(0,i,j,k) dans true_*a_c(i,j,k) void normsup_alla_c(); // diagnostic of minimization double gdsr[NZ]; int nksr; void xtraqsr_init();// init of solar radiation penetration # define rdttra( k ) rdt double r2dt; double atfp1 = 1. - 2. * atfp;// !: asselin time filter coeff. (atfp1= 1-2*atfp) //------->zdfmxl double avt_c = 5.e-4; // ! Kz criterion for the turbocline depth double rho_c = 0.01;// ! density criterion for mixed layer depth double ppgphi0 = 29.; // latitude of first raw column T-point //latitude for the Coriolis or Beta parameter //mode de sauvegarde (==0 3D, ==1 2D surface ==2 un seul point) short savemode=0; short isave=23; short jsave=16; short ksave=0; int n_zdfexp=3; int jphgr_msh=5;//pour ff :beta-plane with regular grid-spacing and rotated domain (GYRE config) // a faire /*nyear = ndastp / 10000 : current year nmonth : current month of the year nyear nday_year : current day of the year nyea ndastp : = nyear*10000 + nmonth*100 + nday*/ #ifdef K_FILTER double pdth_fil=100.;//7200.; double pdtz_fil=100.;//100.;//7200.; double fsaht_fil= 5.e+02;//1000.; double avt_fil=1.2e-05;//1.2e-5; //double lscale= 212000.; // length-scale double norm_fil; //coef de normalisation void normfil();// calcul du coef de normalisation norm_fil void xset_fdt_fco();// set fdt fcoef #endif void cal_ff(); void xdisplay() ; void xistate_gyre(); void xistate_yao_file(char * pth_file_yao); // a faire void xdom_init(); void xflt_rst(); void xsolmat_init(); void xistate_init(int argc, char *argv[]); void xrst_save(int argc, char *argv[]); void xinit_masque_obs(); ////////////////////////////////////////////////////////////// // Les fonctions OBLIGATOIRES //_____________________________________________________________________________ void appli_start (int argc, char *argv[]) { // permet si besoin de prendre la main des le debut de l'appl; printf("////////////////////////////////////////////////////////////////////////\n"); printf("// NEMO/YAO PROJECT //\n"); printf("// M. Berrada 02-2009 J. Brajard 01-2012 //\n"); printf("// LOCEAN-IPSL.UPMC (Paris 6) //\n"); printf("//====================================================================//\n"); printf("\n"); #ifdef K_FILTER printf("option FILTER enable\n"); #else printf("option FILTER disable\n"); #endif #ifdef K_SOLSORYAO printf("option SOLSORYAO enable\n"); #else printf("option SOLSORYAO disable\n"); #endif #ifdef K_OPTIMORDER printf("option OPTIMORDER enable\n"); #else printf("option OPTIMORDER disable\n"); #endif //Vérification du type de réel utilisé //A modifier - pas génial #ifdef YFLOAT #ifndef YYFLOAT printf("Incoherent real type between O_REAL option and ncutil.h\n"); exit(1); #endif #endif #ifdef YDOUBLE #ifndef YYDOUBLE printf("Incoherent real type between O_REAL option and ncutil.h\n"); exit(1); #endif #endif for (int i=0;i>> current forward time : kt=%lf \n",kt); if(Yt==TU+1 && neuler==0) r2dt=rdt; else r2dt=2.*rdt; } //_____________________________________________________________________________ void forward_after (int ctrp) { // permet d'intervenir si besoin après le forward //printf("neuler=%d\n",neuler); } //_____________________________________________________________________________ void backward_before (int ctrp) { // permet d'intervenir si besoin avant le backward if( ndiag_rst==1){ printf("backward_before:: ndiag_rst will be =0\n"); exit(99); } if(Yt==TU+1 && neuler==0) r2dt=rdt; else r2dt=2.*rdt; } //_____________________________________________________________________________ void backward_after (int ctrp) { // permet d'intervenir si besoin apres le backward } //_____________________________________________________________________________ short select_io(int indic, char *nmmod, int sortie, int iaxe, int jaxe, int kaxe, int pdt, YREAL *val) { // Pour faire des selections sur les fonctions d'entrees sorties de Yao ou autre en generale. // Doit retourner 1 si l'element dont les caracteristiques sont passes en parametre doit // etre retenu, pris, selectionne ; et 0 sinon // indic permet de savoir de quelle instance provient l'appel : // =YIO_LOADSTATE => appel via loadstate // =YIO_SAVESTATE => appel via savestate // =YIO_LOADOBS => appel via loadobs // =YIO_OUTOOBS => appel via outoobs if(indic==YIO_SAVESTATE && !strcmp(nmmod,"tb") && savemode==1){ if(kaxe==0 ) return(1); else return(0); } if(indic==YIO_SAVESTATE && !strcmp(nmmod,"tb") && savemode==2){ if(iaxe==isave && jaxe==jsave && kaxe==ksave) return(1); else return(0); } if(indic==YIO_SAVESTATE && (!strcmp(nmmod,"tb") || !strcmp(nmmod,"sb")) && savemode==3) return(masque_obs(iaxe,jaxe)); /* if(indic==YIO_OUTOOBS){ if(iaxe==9 && jaxe==9 ) return(1); else return(0); }*/ return(1); } //======================== Fin des foncts obligatoires ========================= //============================================================================== //Les autres fonctions //_____________________________________________________________________________ void xdisplay() { normsup_alla_c(); } //_____________________________________________________________________________ #ifdef K_FILTER //_____________________________________________________________________________ void normfil(){ double norm_zfil=1.; double norm_lfil=1.; norm_lfil= 2.*sqrt(PI*fsaht_fil*pdth_fil*2.*(OFTL-1.));// schema explicite norm_zfil= 2.*sqrt(PI*avt_fil*pdtz_fil*2.*(OFTZ-1.));// schema explicite #ifdef K_FZIMP // schema implicite double fact=1.; int M=2*(OFTZ-1); for(int i=1;igam tab_luzimp( i, j, k , 1 ) = 0.;// 1-->bet } double a0,a1; double a2,b2,c2; for(int j=1;jbet #endif tab_luzimp( i, j, 0 , 0 ) = 0.;// 0-->gam tab_luzimp( i, j, 0 , 1 ) = b2;// 1-->bet for(int k=1;kNZ-1) nksr=NZ-1; } } //_____________________________________________________________________________ void cal_ff() //domhgr.F90 Coriolis parameter { if(jphgr_msh==5){// beta-plane and rotated domain (gyre configuration) double zbeta = 2.*omega*cos(rad*ppgphi0)/ra ;// beta at latitude ppgphi0 double zphi0 = 15.;//latitude of the first row F-points double zf0 = 2.*omega*sin(rad*zphi0);// compute f0 1st point south for(int i=0;i