source: branches/branche-mb/src/dta_lfexp.h @ 68

Last change on this file since 68 was 68, checked in by berrada, 13 years ago

commit final avec filtre
  • Property svn:eol-style set to native
File size: 3.6 KB
Line 
1
2/***************************************************************************
3                 module classe dta_lfexp.h  -  description
4***************************************************************************/
5// Mohamed Berrada
6// locean-ipsl.upmc, Paris, August 9, 2011
7//===========================================================================
8//                          methode forward
9forward (YREAL x1,YREAL x2,YREAL x3,YREAL x4,YREAL x5,YREAL x6)
10{ /*    1       from    nu                      1  i    j    k
11        2       from    dta_lfexp               1  i    j    k     t-1
12        3       from    dta_lfexp               1  i-1  j    k     t-1
13        4       from    dta_lfexp               1  i+1  j    k     t-1
14        5       from    dta_lfexp               1  i    j-1  k     t-1
15        6       from    dta_lfexp               1  i    j+1  k     t-1 */
16  //  if (Yi==5 && Yj==5)
17  //printf("dta_lfil : %d\n",Yt);
18
19  double zbtr = zbtr2(Yi,Yj);
20  if(Yt==1){
21    double wm12=sqrt(zbtr);//1./(106000.); //1/sqrt(ds1*ds2)
22    YS1=wm12*x1;
23  }
24  else{
25    if(Yi==0 || Yi==NX-1 || Yj==0 || Yj==NY-1 || Yk==NZ-1){
26      YS1=0.;
27    }
28    else{
29      double bi   = fsaht_fil * pdth_fil * umask(Yi,Yj,Yk) * ze1ur(Yi,Yj);
30      double bim1 = fsaht_fil * pdth_fil * umask(Yi-1,Yj,Yk) * ze1ur(Yi-1,Yj);
31      double bj   = fsaht_fil * pdth_fil * vmask(Yi,Yj,Yk) * ze2vr(Yi,Yj);
32      double bjm1 = fsaht_fil * pdth_fil * vmask(Yi,Yj-1,Yk) * ze2vr(Yi,Yj-1);
33
34      double tip1,tim1,tjp1,tjm1, tij;
35      //cdt au bord i
36      if (Yi==1){
37        tim1= 0.;
38#ifdef K_NEU
39        tim1= bim1 * x2;
40#endif
41      }
42      else{
43        tim1= bim1 * x3;
44      }
45      if (Yi==NX-2){
46        tip1= 0.;
47#ifdef K_NEU
48        tip1= bi * x2;
49#endif
50      }
51      else{
52        tip1= bi * x4;
53      }
54      //cdt au bord j
55      if (Yj==1){
56        tjm1= 0.;
57#ifdef K_NEU
58        tjm1= bjm1 * x2;
59#endif
60      }
61      else{
62        tjm1= bjm1 * x5;
63      }
64      if (Yj==NY-2){
65        tjp1= 0.;
66#ifdef K_NEU
67        tjp1= bj * x2;
68#endif
69      }
70      else{
71        tjp1= bj * x6;
72      }
73
74      tij= (1./zbtr-bi-bim1-bj-bjm1) * x2;
75      YS1 = zbtr*(tip1 + tim1 + tjp1 + tjm1 + tij);
76    }
77  }
78}
79
80//===========================================================================
81//                         methode  backward
82
83backward (YREAL x1,YREAL x2,YREAL x3,YREAL x4,YREAL x5,YREAL x6)
84{
85  YJ1I1=0.;     YJ1I2=0.;    YJ1I3=0.;    YJ1I4=0.;    YJ1I5=0.;    YJ1I6=0.;
86  double zbtr = zbtr2(Yi,Yj);
87  if(Yt==1){
88    double wm12=sqrt(zbtr);//1./(106000.); //1/sqrt(ds1*ds2)
89    YJ1I1=wm12;
90  }
91  else{
92    if(Yi==0 || Yi==NX-1 || Yj==0 || Yj==NY-1 || Yk==NZ-1){
93      YJ1I1=0.;     YJ1I2=0.;    YJ1I3=0.;    YJ1I4=0.;    YJ1I5=0.;    YJ1I6=0.;
94    }
95    else{
96      double bi   = fsaht_fil * pdth_fil * umask(Yi,Yj,Yk) * ze1ur(Yi,Yj);
97      double bim1 = fsaht_fil * pdth_fil * umask(Yi-1,Yj,Yk) * ze1ur(Yi-1,Yj);
98      double bj   = fsaht_fil * pdth_fil * vmask(Yi,Yj,Yk) * ze2vr(Yi,Yj);
99      double bjm1 = fsaht_fil * pdth_fil * vmask(Yi,Yj-1,Yk) * ze2vr(Yi,Yj-1);
100
101      double dtip1_x4=0.,dtim1_x3=0.,dtjp1_x6=0.,dtjm1_x5=0.;
102      double dtip1_x2=0.,dtim1_x2=0.,dtjp1_x2=0.,dtjm1_x2=0., dtij_x2=0.;
103      //cdt au bord i
104      if (Yi==1){
105        dtim1_x2= 0.;
106#ifdef K_NEU
107        dtim1_x2= bim1;
108#endif
109      }
110      else{
111        dtim1_x3= bim1;
112      }
113      if (Yi==NX-2){
114        dtip1_x2= 0.;
115#ifdef K_NEU
116        dtip1_x2= bi;
117#endif
118      }
119      else{
120        dtip1_x4= bi;
121      }
122      //cdt au bord j
123      if (Yj==1){
124        dtjm1_x2= 0.;
125#ifdef K_NEU
126        dtjm1_x2= bjm1;
127#endif
128      }
129      else{
130        dtjm1_x5= bjm1;
131      }
132      if (Yj==NY-2){
133        dtjp1_x2= 0.;
134#ifdef K_NEU
135        dtjp1_x2= bj;
136#endif
137      }
138      else{
139        dtjp1_x6= bj;
140      }
141      dtij_x2= (1./zbtr-bi-bim1-bj-bjm1);
142
143      YJ1I1 = 0.;
144      YJ1I2 = zbtr*(dtip1_x2 + dtim1_x2 + dtjp1_x2 + dtjm1_x2 + dtij_x2);
145      YJ1I3 = zbtr*(dtim1_x3);
146      YJ1I4 = zbtr*(dtip1_x4);
147      YJ1I5 = zbtr*(dtjm1_x5);
148      YJ1I6 = zbtr*(dtjp1_x6);
149    }
150  }
151}
Note: See TracBrowser for help on using the repository browser.