[58] | 1 | /*************************************************************************** |
---|
| 2 | module classe filtre_lap.h - description |
---|
| 3 | ***************************************************************************/ |
---|
| 4 | // Mohamed Berrada |
---|
| 5 | // locean-ipsl.upmc, Paris, July 19, 2011 |
---|
| 6 | //=========================================================================== |
---|
| 7 | // methode forward |
---|
| 8 | #define dta_c( i, j, k) (dta_c[(k)*(NY*NX)+(j)*(NX)+(i)]) // l'increment de t calculer par le fitre |
---|
| 9 | #define ztu( i, j, k) (ztu[(k)*(NY*NX)+(j)*(NX)+(i)]) |
---|
| 10 | #define ztv( i, j, k) (ztv[(k)*(NY*NX)+(j)*(NX)+(i)]) |
---|
| 11 | #define zwy( i, j, k) (zwy[(k)*(NY*NX)+(j)*(NX)+(i)]) |
---|
| 12 | |
---|
| 13 | double dta_c[NZ*NY*NX]; // l'increment de t calculer par le fitre |
---|
| 14 | double ztu[NZ*NY*NX]; // |
---|
| 15 | double ztv[NZ*NY*NX]; // |
---|
| 16 | double zwy[NZ*NY*NX]; // |
---|
| 17 | |
---|
| 18 | forward () |
---|
| 19 | { |
---|
| 20 | int n,ji,jj,jk; |
---|
| 21 | double zabe1,zabe2,zbtr,zta; |
---|
| 22 | double zave3r; |
---|
| 23 | double fsahtu_fil; // coef diff hori |
---|
| 24 | double fsahtv_fil; // coef diff hori |
---|
| 25 | double avt_fil; // coef diffu verti |
---|
| 26 | double pdtv_fil, pdtz_fil; // pas du temps verti et horiz |
---|
| 27 | //boundarie conditions |
---|
| 28 | for(jk=0;jk<NZ;jk++) |
---|
| 29 | for(jj=0;jj<NY;jj++){ |
---|
| 30 | YS1_nu(0,jj,jk)=0.; YS1_nu(NX-1,jj,jk)=0.; |
---|
| 31 | } |
---|
| 32 | for(jk=0;jk<NZ;jk++) |
---|
| 33 | for(ji=0;ji<NX;ji++){ |
---|
| 34 | YS1_nu(ji,0,jk)=0.; YS1_nu(ji,NY-1,jk)=0.; |
---|
| 35 | } |
---|
| 36 | //end boundarie conditions |
---|
| 37 | |
---|
| 38 | // initialisation de l'increment et des tableaux temp |
---|
| 39 | for(jk=0;jk<NZ;jk++) |
---|
| 40 | for(jj=0;jj<NY;jj++) |
---|
| 41 | for(ji=0;ji<NX;ji++){ |
---|
| 42 | dta_c(ji,jj,jk)=YS1_nu(ji,jj,jk); |
---|
| 43 | ztu(ji,jj,jk)=0.; |
---|
| 44 | ztv(ji,jj,jk)=0.; |
---|
| 45 | zwy(ji,jj,jk)=0.; |
---|
| 46 | } |
---|
| 47 | |
---|
| 48 | |
---|
| 49 | // sqrt W^(-1/2) a faire ici |
---|
| 50 | |
---|
| 51 | // sqrt horiz diffu Lh^1/2 |
---|
| 52 | for(int n=0;n<MHs2;n++){ |
---|
| 53 | for(jk = 0; jk<NZ-1,jk++){ |
---|
| 54 | //! 1. First derivative (gradient) |
---|
| 55 | for( jj = 1; jj<NY-1;jj++){ |
---|
| 56 | for(ji = 1; ji<NX-1;ji++){ |
---|
| 57 | zabe1 = fsahtu_fil * umask(ji,jj,jk) * ze1ur(ji,jj); |
---|
| 58 | zabe2 = fsahtv_fil * vmask(ji,jj,jk) * ze2vr(ji,jj); |
---|
| 59 | ztu(ji,jj,jk) = zabe1 * ( dta_c(ji+1,jj ,jk) - dta_c(ji,jj,jk) ); |
---|
| 60 | ztv(ji,jj,jk) = zabe2 * ( dta_c(ji ,jj+1,jk) - dta_c(ji,jj,jk) ); |
---|
| 61 | } |
---|
| 62 | } |
---|
| 63 | //! 2. Second derivative (divergence) |
---|
| 64 | //! -------------------- |
---|
| 65 | for(jj = 1; jj<NY-1; jj++){ |
---|
| 66 | for(ji = 1; ji<NX-1;ji++){ |
---|
| 67 | zbtr = zbtr2(ji,jj); |
---|
| 68 | //! horizontal diffusive trends |
---|
| 69 | zta = pdth_fil*zbtr*(ztu(ji,jj,jk)- ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)); |
---|
| 70 | //! add it to the general tracer trends |
---|
| 71 | dta_c(ji,jj,jk) = dta_c(ji,jj,jk) + zta; |
---|
| 72 | } |
---|
| 73 | } |
---|
| 74 | } |
---|
| 75 | } |
---|
| 76 | |
---|
| 77 | // sqrt verti diffu Lv^1/2 |
---|
| 78 | for(n=0;n<MZs2;n++){ |
---|
| 79 | for(jk=1;jk<NZ;jk++){ //! first vertical derivative (gradient) |
---|
| 80 | for(jj = 1;jj<NY-1;jj++){ |
---|
| 81 | for(ji = 1;ji<NX-1;ji++){ |
---|
| 82 | zave3r = 1.e0 / fse3w(ji,jj,jk); //! 1 sur le pas veritcal au point w (Arakawa mesh) |
---|
| 83 | zwy(ji,jj,jk) = avt_fil * ( dta_c(ji,jj,jk-1) - dta_c(ji,jj,jk) ) * zave3r; |
---|
| 84 | } |
---|
| 85 | } |
---|
| 86 | } |
---|
| 87 | |
---|
| 88 | for(jk=0;jk<NZ-1;jk++){ //! second vertical derivative (divergence) |
---|
| 89 | for(jj=1;jj<NY-1;jj++){ |
---|
| 90 | for(ji=1;ji<NX-1;ji++){ |
---|
| 91 | ze3tr=1.e0 / fse3t(ji,jj,jk);// ! 1 sur le pas veritcal au point t |
---|
| 92 | zta = pdtz_fil* ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) * ze3tr; |
---|
| 93 | dta_c(ji,jj,jk) = dta_c(ji,jj,jk) + zta; |
---|
| 94 | } |
---|
| 95 | } |
---|
| 96 | } |
---|
| 97 | } |
---|
| 98 | // |
---|
| 99 | } |
---|
| 100 | |
---|
| 101 | //=========================================================================== |
---|
| 102 | // methode backward and linward |
---|
| 103 | backward () |
---|
| 104 | { |
---|
| 105 | |
---|
| 106 | if(Ycurward==BACKWARD){ |
---|
| 107 | } |
---|
| 108 | else if(Ycurward==LINWARD){ |
---|
| 109 | } |
---|
| 110 | } |
---|