source: soft/giet_vm/applications/ocean/jacobcalc2.C @ 769

Last change on this file since 769 was 598, checked in by guerin, 9 years ago

ocean: fix app broken by r589
File size: 10.5 KB
Line 
1/*************************************************************************/
2/*                                                                       */
3/*  Copyright (c) 1994 Stanford University                               */
4/*                                                                       */
5/*  All rights reserved.                                                 */
6/*                                                                       */
7/*  Permission is given to use, copy, and modify this software for any   */
8/*  non-commercial purpose as long as this copyright notice is not       */
9/*  removed.  All other uses, including redistribution in whole or in    */
10/*  part, are forbidden without prior written permission.                */
11/*                                                                       */
12/*  This software is provided with absolutely no warranty and no         */
13/*  support.                                                             */
14/*                                                                       */
15/*************************************************************************/
16
17/* Does the arakawa jacobian calculation (of the x and y matrices,
18   putting the results in the z matrix) for a subblock. */
19
20EXTERN_ENV
21
22#include <stdio.h>
23#include <math.h>
24
25#include "decs.h"
26
27void jacobcalc2(double ****x, double ****y, double ****z, long psiindex, long pid, long firstrow, long lastrow, long firstcol, long lastcol)
28{
29    double f1;
30    double f2;
31    double f3;
32    double f4;
33    double f5;
34    double f6;
35    double f7;
36    double f8;
37    long iindex;
38    long indexp1;
39    long indexm1;
40    long im1;
41    long ip1;
42    long i;
43    long j;
44    long jj;
45    double **t2a;
46    double **t2b;
47    double **t2c;
48    double *t1a;
49    double *t1b;
50    double *t1c;
51    double *t1d;
52    double *t1e;
53    double *t1f;
54    double *t1g;
55
56    t2a = z[pid][psiindex];
57    if ((gp[pid].neighbors[UP] == -1) && (gp[pid].neighbors[LEFT] == -1)) {
58        t2a[0][0] = 0.0;
59    }
60    if ((gp[pid].neighbors[DOWN] == -1) && (gp[pid].neighbors[LEFT] == -1)) {
61        t2a[im - 1][0] = 0.0;
62    }
63    if ((gp[pid].neighbors[UP] == -1) && (gp[pid].neighbors[RIGHT] == -1)) {
64        t2a[0][jm - 1] = 0.0;
65    }
66    if ((gp[pid].neighbors[DOWN] == -1) && (gp[pid].neighbors[RIGHT] == -1)) {
67        t2a[im - 1][jm - 1] = 0.0;
68    }
69
70    t2a = x[pid][psiindex];
71    jj = gp[pid].neighbors[UPLEFT];
72    if (jj != -1) {
73        t2a[0][0] = x[jj][psiindex][im - 2][jm - 2];
74    }
75    jj = gp[pid].neighbors[UPRIGHT];
76    if (jj != -1) {
77        t2a[0][jm - 1] = x[jj][psiindex][im - 2][1];
78    }
79    jj = gp[pid].neighbors[DOWNLEFT];
80    if (jj != -1) {
81        t2a[im - 1][0] = x[jj][psiindex][1][jm - 2];
82    }
83    jj = gp[pid].neighbors[DOWNRIGHT];
84    if (jj != -1) {
85        t2a[im - 1][jm - 1] = x[jj][psiindex][1][1];
86    }
87
88    t2a = y[pid][psiindex];
89    jj = gp[pid].neighbors[UPLEFT];
90    if (jj != -1) {
91        t2a[0][0] = y[jj][psiindex][im - 2][jm - 2];
92    }
93    jj = gp[pid].neighbors[UPRIGHT];
94    if (jj != -1) {
95        t2a[0][jm - 1] = y[jj][psiindex][im - 2][1];
96    }
97    jj = gp[pid].neighbors[DOWNLEFT];
98    if (jj != -1) {
99        t2a[im - 1][0] = y[jj][psiindex][1][jm - 2];
100    }
101    jj = gp[pid].neighbors[DOWNRIGHT];
102    if (jj != -1) {
103        t2a[im - 1][jm - 1] = y[jj][psiindex][1][1];
104    }
105
106    t2a = x[pid][psiindex];
107    if (gp[pid].neighbors[UP] == -1) {
108        jj = gp[pid].neighbors[LEFT];
109        if (jj != -1) {
110            t2a[0][0] = x[jj][psiindex][0][jm - 2];
111        } else {
112            jj = gp[pid].neighbors[DOWN];
113            if (jj != -1) {
114                t2a[im - 1][0] = x[jj][psiindex][1][0];
115            }
116        }
117        jj = gp[pid].neighbors[RIGHT];
118        if (jj != -1) {
119            t2a[0][jm - 1] = x[jj][psiindex][0][1];
120        } else {
121            jj = gp[pid].neighbors[DOWN];
122            if (jj != -1) {
123                t2a[im - 1][jm - 1] = x[jj][psiindex][1][jm - 1];
124            }
125        }
126    } else if (gp[pid].neighbors[DOWN] == -1) {
127        jj = gp[pid].neighbors[LEFT];
128        if (jj != -1) {
129            t2a[im - 1][0] = x[jj][psiindex][im - 1][jm - 2];
130        } else {
131            jj = gp[pid].neighbors[UP];
132            if (jj != -1) {
133                t2a[0][0] = x[jj][psiindex][im - 2][0];
134            }
135        }
136        jj = gp[pid].neighbors[RIGHT];
137        if (jj != -1) {
138            t2a[im - 1][jm - 1] = x[jj][psiindex][im - 1][1];
139        } else {
140            jj = gp[pid].neighbors[UP];
141            if (jj != -1) {
142                t2a[0][jm - 1] = x[jj][psiindex][im - 2][jm - 1];
143            }
144        }
145    } else if (gp[pid].neighbors[LEFT] == -1) {
146        jj = gp[pid].neighbors[UP];
147        if (jj != -1) {
148            t2a[0][0] = x[jj][psiindex][im - 2][0];
149        }
150        jj = gp[pid].neighbors[DOWN];
151        if (jj != -1) {
152            t2a[im - 1][0] = x[jj][psiindex][1][0];
153        }
154    } else if (gp[pid].neighbors[RIGHT] == -1) {
155        jj = gp[pid].neighbors[UP];
156        if (jj != -1) {
157            t2a[0][jm - 1] = x[jj][psiindex][im - 2][jm - 1];
158        }
159        jj = gp[pid].neighbors[DOWN];
160        if (jj != -1) {
161            t2a[im - 1][jm - 1] = x[jj][psiindex][1][jm - 1];
162        }
163    }
164
165    t2a = y[pid][psiindex];
166    if (gp[pid].neighbors[UP] == -1) {
167        jj = gp[pid].neighbors[LEFT];
168        if (jj != -1) {
169            t2a[0][0] = y[jj][psiindex][0][jm - 2];
170        } else {
171            jj = gp[pid].neighbors[DOWN];
172            if (jj != -1) {
173                t2a[im - 1][0] = y[jj][psiindex][1][0];
174            }
175        }
176        jj = gp[pid].neighbors[RIGHT];
177        if (jj != -1) {
178            t2a[0][jm - 1] = y[jj][psiindex][0][1];
179        } else {
180            jj = gp[pid].neighbors[DOWN];
181            if (jj != -1) {
182                t2a[im - 1][jm - 1] = y[jj][psiindex][1][jm - 1];
183            }
184        }
185    } else if (gp[pid].neighbors[DOWN] == -1) {
186        jj = gp[pid].neighbors[LEFT];
187        if (jj != -1) {
188            t2a[im - 1][0] = y[jj][psiindex][im - 1][jm - 2];
189        } else {
190            jj = gp[pid].neighbors[UP];
191            if (jj != -1) {
192                t2a[0][0] = y[jj][psiindex][im - 2][0];
193            }
194        }
195        jj = gp[pid].neighbors[RIGHT];
196        if (jj != -1) {
197            t2a[im - 1][jm - 1] = y[jj][psiindex][im - 1][1];
198        } else {
199            jj = gp[pid].neighbors[UP];
200            if (jj != -1) {
201                t2a[0][jm - 1] = y[jj][psiindex][im - 2][jm - 1];
202            }
203        }
204    } else if (gp[pid].neighbors[LEFT] == -1) {
205        jj = gp[pid].neighbors[UP];
206        if (jj != -1) {
207            t2a[0][0] = y[jj][psiindex][im - 2][0];
208        }
209        jj = gp[pid].neighbors[DOWN];
210        if (jj != -1) {
211            t2a[im - 1][0] = y[jj][psiindex][1][0];
212        }
213    } else if (gp[pid].neighbors[RIGHT] == -1) {
214        jj = gp[pid].neighbors[UP];
215        if (jj != -1) {
216            t2a[0][jm - 1] = y[jj][psiindex][im - 2][jm - 1];
217        }
218        jj = gp[pid].neighbors[DOWN];
219        if (jj != -1) {
220            t2a[im - 1][jm - 1] = y[jj][psiindex][1][jm - 1];
221        }
222    }
223
224    t2a = y[pid][psiindex];
225    j = gp[pid].neighbors[UP];
226    if (j != -1) {
227        t1a = (double *) t2a[0];
228        t1b = (double *) y[j][psiindex][im - 2];
229        for (i = 1; i <= lastcol; i++) {
230            t1a[i] = t1b[i];
231        }
232    }
233    j = gp[pid].neighbors[DOWN];
234    if (j != -1) {
235        t1a = (double *) t2a[im - 1];
236        t1b = (double *) y[j][psiindex][1];
237        for (i = 1; i <= lastcol; i++) {
238            t1a[i] = t1b[i];
239        }
240    }
241    j = gp[pid].neighbors[LEFT];
242    if (j != -1) {
243        t2b = y[j][psiindex];
244        for (i = 1; i <= lastrow; i++) {
245            t2a[i][0] = t2b[i][jm - 2];
246        }
247    }
248    j = gp[pid].neighbors[RIGHT];
249    if (j != -1) {
250        t2b = y[j][psiindex];
251        for (i = 1; i <= lastrow; i++) {
252            t2a[i][jm - 1] = t2b[i][1];
253        }
254    }
255
256    t2a = x[pid][psiindex];
257    j = gp[pid].neighbors[UP];
258    if (j != -1) {
259        t1a = (double *) t2a[0];
260        t1b = (double *) x[j][psiindex][im - 2];
261        for (i = 1; i <= lastcol; i++) {
262            t1a[i] = t1b[i];
263        }
264    }
265    j = gp[pid].neighbors[DOWN];
266    if (j != -1) {
267        t1a = (double *) t2a[im - 1];
268        t1b = (double *) x[j][psiindex][1];
269        for (i = 1; i <= lastcol; i++) {
270            t1a[i] = t1b[i];
271        }
272    }
273    j = gp[pid].neighbors[LEFT];
274    if (j != -1) {
275        t2b = x[j][psiindex];
276        for (i = 1; i <= lastrow; i++) {
277            t2a[i][0] = t2b[i][jm - 2];
278        }
279    }
280    j = gp[pid].neighbors[RIGHT];
281    if (j != -1) {
282        t2b = x[j][psiindex];
283        for (i = 1; i <= lastrow; i++) {
284            t2a[i][jm - 1] = t2b[i][1];
285        }
286    }
287
288    t2a = x[pid][psiindex];
289    t2b = y[pid][psiindex];
290    t2c = z[pid][psiindex];
291    for (i = firstrow; i <= lastrow; i++) {
292        ip1 = i + 1;
293        im1 = i - 1;
294        t1a = (double *) t2a[i];
295        t1b = (double *) t2b[i];
296        t1c = (double *) t2c[i];
297        t1d = (double *) t2b[ip1];
298        t1e = (double *) t2b[im1];
299        t1f = (double *) t2a[ip1];
300        t1g = (double *) t2a[im1];
301        for (iindex = firstcol; iindex <= lastcol; iindex++) {
302            indexp1 = iindex + 1;
303            indexm1 = iindex - 1;
304            f1 = (t1b[indexm1] + t1d[indexm1] - t1b[indexp1] - t1d[indexp1]) * (t1f[iindex] - t1a[iindex]);
305            f2 = (t1e[indexm1] + t1b[indexm1] - t1e[indexp1] - t1b[indexp1]) * (t1a[iindex] - t1g[iindex]);
306            f3 = (t1d[iindex] + t1d[indexp1] - t1e[iindex] - t1e[indexp1]) * (t1a[indexp1] - t1a[iindex]);
307            f4 = (t1d[indexm1] + t1d[iindex] - t1e[indexm1] - t1e[iindex]) * (t1a[iindex] - t1a[indexm1]);
308            f5 = (t1d[iindex] - t1b[indexp1]) * (t1f[indexp1] - t1a[iindex]);
309            f6 = (t1b[indexm1] - t1e[iindex]) * (t1a[iindex] - t1g[indexm1]);
310            f7 = (t1b[indexp1] - t1e[iindex]) * (t1g[indexp1] - t1a[iindex]);
311            f8 = (t1d[iindex] - t1b[indexm1]) * (t1a[iindex] - t1f[indexm1]);
312
313            t1c[iindex] = factjacob * (f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8);
314
315        }
316    }
317
318    if (gp[pid].neighbors[UP] == -1) {
319        t1c = (double *) t2c[0];
320        for (j = firstcol; j <= lastcol; j++) {
321            t1c[j] = 0.0;
322        }
323    }
324    if (gp[pid].neighbors[DOWN] == -1) {
325        t1c = (double *) t2c[im - 1];
326        for (j = firstcol; j <= lastcol; j++) {
327            t1c[j] = 0.0;
328        }
329    }
330    if (gp[pid].neighbors[LEFT] == -1) {
331        for (j = firstrow; j <= lastrow; j++) {
332            t2c[j][0] = 0.0;
333        }
334    }
335    if (gp[pid].neighbors[RIGHT] == -1) {
336        for (j = firstrow; j <= lastrow; j++) {
337            t2c[j][jm - 1] = 0.0;
338        }
339    }
340
341}
Note: See TracBrowser for help on using the repository browser.