source: soft/giet_vm/applications/ocean/jacobcalc.C @ 667

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

ocean: fix app broken by r589
File size: 10.2 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 jacobcalc(double ***x, double ***y, double ***z, 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 = (double **) z[pid];
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 = (double **) x[pid];
71    jj = gp[pid].neighbors[UPLEFT];
72    if (jj != -1) {
73        t2a[0][0] = x[jj][im - 2][jm - 2];
74    }
75    jj = gp[pid].neighbors[UPRIGHT];
76    if (jj != -1) {
77        t2a[0][jm - 1] = x[jj][im - 2][1];
78    }
79    jj = gp[pid].neighbors[DOWNLEFT];
80    if (jj != -1) {
81        t2a[im - 1][0] = x[jj][1][jm - 2];
82    }
83    jj = gp[pid].neighbors[DOWNRIGHT];
84    if (jj != -1) {
85        t2a[im - 1][jm - 1] = x[jj][1][1];
86    }
87
88    t2a = (double **) y[pid];
89    jj = gp[pid].neighbors[UPLEFT];
90    if (jj != -1) {
91        t2a[0][0] = y[jj][im - 2][jm - 2];
92    }
93    jj = gp[pid].neighbors[UPRIGHT];
94    if (jj != -1) {
95        t2a[0][jm - 1] = y[jj][im - 2][1];
96    }
97    jj = gp[pid].neighbors[DOWNLEFT];
98    if (jj != -1) {
99        t2a[im - 1][0] = y[jj][1][jm - 2];
100    }
101    jj = gp[pid].neighbors[DOWNRIGHT];
102    if (jj != -1) {
103        t2a[im - 1][jm - 1] = y[jj][1][1];
104    }
105
106    t2a = (double **) x[pid];
107    if (gp[pid].neighbors[UP] == -1) {
108        jj = gp[pid].neighbors[LEFT];
109        if (jj != -1) {
110            t2a[0][0] = x[jj][0][jm - 2];
111        } else {
112            jj = gp[pid].neighbors[DOWN];
113            if (jj != -1) {
114                t2a[im - 1][0] = x[jj][1][0];
115            }
116        }
117        jj = gp[pid].neighbors[RIGHT];
118        if (jj != -1) {
119            t2a[0][jm - 1] = x[jj][0][1];
120        } else {
121            jj = gp[pid].neighbors[DOWN];
122            if (jj != -1) {
123                t2a[im - 1][jm - 1] = x[jj][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][im - 1][jm - 2];
130        } else {
131            jj = gp[pid].neighbors[UP];
132            if (jj != -1) {
133                t2a[0][0] = x[jj][im - 2][0];
134            }
135        }
136        jj = gp[pid].neighbors[RIGHT];
137        if (jj != -1) {
138            t2a[im - 1][jm - 1] = x[jj][im - 1][1];
139        } else {
140            jj = gp[pid].neighbors[UP];
141            if (jj != -1) {
142                t2a[0][jm - 1] = x[jj][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][im - 2][0];
149        }
150        jj = gp[pid].neighbors[DOWN];
151        if (jj != -1) {
152            t2a[im - 1][0] = x[jj][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][im - 2][jm - 1];
158        }
159        jj = gp[pid].neighbors[DOWN];
160        if (jj != -1) {
161            t2a[im - 1][jm - 1] = x[jj][1][jm - 1];
162        }
163    }
164
165    t2a = (double **) y[pid];
166    if (gp[pid].neighbors[UP] == -1) {
167        jj = gp[pid].neighbors[LEFT];
168        if (jj != -1) {
169            t2a[0][0] = y[jj][0][jm - 2];
170        } else {
171            jj = gp[pid].neighbors[DOWN];
172            if (jj != -1) {
173                t2a[im - 1][0] = y[jj][1][0];
174            }
175        }
176        jj = gp[pid].neighbors[RIGHT];
177        if (jj != -1) {
178            t2a[0][jm - 1] = y[jj][0][1];
179        } else {
180            jj = gp[pid].neighbors[DOWN];
181            if (jj != -1) {
182                t2a[im - 1][jm - 1] = y[jj][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][im - 1][jm - 2];
189        } else {
190            jj = gp[pid].neighbors[UP];
191            if (jj != -1) {
192                t2a[0][0] = y[jj][im - 2][0];
193            }
194        }
195        jj = gp[pid].neighbors[RIGHT];
196        if (jj != -1) {
197            t2a[im - 1][jm - 1] = y[jj][im - 1][1];
198        } else {
199            jj = gp[pid].neighbors[UP];
200            if (jj != -1) {
201                t2a[0][jm - 1] = y[jj][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][im - 2][0];
208        }
209        jj = gp[pid].neighbors[DOWN];
210        if (jj != -1) {
211            t2a[im - 1][0] = y[jj][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][im - 2][jm - 1];
217        }
218        jj = gp[pid].neighbors[DOWN];
219        if (jj != -1) {
220            t2a[im - 1][jm - 1] = y[jj][1][jm - 1];
221        }
222    }
223
224    j = gp[pid].neighbors[UP];
225    if (j != -1) {
226        t1a = (double *) t2a[0];
227        t1b = (double *) y[j][im - 2];
228        for (i = 1; i <= lastcol; i++) {
229            t1a[i] = t1b[i];
230        }
231    }
232    j = gp[pid].neighbors[DOWN];
233    if (j != -1) {
234        t1a = (double *) t2a[im - 1];
235        t1b = (double *) y[j][1];
236        for (i = 1; i <= lastcol; i++) {
237            t1a[i] = t1b[i];
238        }
239    }
240    j = gp[pid].neighbors[LEFT];
241    if (j != -1) {
242        t2b = (double **) y[j];
243        for (i = 1; i <= lastrow; i++) {
244            t2a[i][0] = t2b[i][jm - 2];
245        }
246    }
247    j = gp[pid].neighbors[RIGHT];
248    if (j != -1) {
249        t2b = (double **) y[j];
250        for (i = 1; i <= lastrow; i++) {
251            t2a[i][jm - 1] = t2b[i][1];
252        }
253    }
254
255    t2a = (double **) x[pid];
256    j = gp[pid].neighbors[UP];
257    if (j != -1) {
258        t1a = (double *) t2a[0];
259        t1b = (double *) x[j][im - 2];
260        for (i = 1; i <= lastcol; i++) {
261            t1a[i] = t1b[i];
262        }
263    }
264    j = gp[pid].neighbors[DOWN];
265    if (j != -1) {
266        t1a = (double *) t2a[im - 1];
267        t1b = (double *) x[j][1];
268        for (i = 1; i <= lastcol; i++) {
269            t1a[i] = t1b[i];
270        }
271    }
272    j = gp[pid].neighbors[LEFT];
273    if (j != -1) {
274        t2b = (double **) x[j];
275        for (i = 1; i <= lastrow; i++) {
276            t2a[i][0] = t2b[i][jm - 2];
277        }
278    }
279    j = gp[pid].neighbors[RIGHT];
280    if (j != -1) {
281        t2b = (double **) x[j];
282        for (i = 1; i <= lastrow; i++) {
283            t2a[i][jm - 1] = t2b[i][1];
284        }
285    }
286
287    t2a = (double **) x[pid];
288    t2b = (double **) y[pid];
289    t2c = (double **) z[pid];
290    for (i = firstrow; i <= lastrow; i++) {
291        ip1 = i + 1;
292        im1 = i - 1;
293        t1a = (double *) t2a[i];
294        t1b = (double *) t2b[i];
295        t1c = (double *) t2c[i];
296        t1d = (double *) t2b[ip1];
297        t1e = (double *) t2b[im1];
298        t1f = (double *) t2a[ip1];
299        t1g = (double *) t2a[im1];
300        for (iindex = firstcol; iindex <= lastcol; iindex++) {
301            indexp1 = iindex + 1;
302            indexm1 = iindex - 1;
303            f1 = (t1b[indexm1] + t1d[indexm1] - t1b[indexp1] - t1d[indexp1]) * (t1f[iindex] - t1a[iindex]);
304            f2 = (t1e[indexm1] + t1b[indexm1] - t1e[indexp1] - t1b[indexp1]) * (t1a[iindex] - t1g[iindex]);
305            f3 = (t1d[iindex] + t1d[indexp1] - t1e[iindex] - t1e[indexp1]) * (t1a[indexp1] - t1a[iindex]);
306            f4 = (t1d[indexm1] + t1d[iindex] - t1e[indexm1] - t1e[iindex]) * (t1a[iindex] - t1a[indexm1]);
307            f5 = (t1d[iindex] - t1b[indexp1]) * (t1f[indexp1] - t1a[iindex]);
308            f6 = (t1b[indexm1] - t1e[iindex]) * (t1a[iindex] - t1g[indexm1]);
309            f7 = (t1b[indexp1] - t1e[iindex]) * (t1g[indexp1] - t1a[iindex]);
310            f8 = (t1d[iindex] - t1b[indexm1]) * (t1a[iindex] - t1f[indexm1]);
311
312            t1c[iindex] = factjacob * (f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8);
313        }
314    }
315
316    if (gp[pid].neighbors[UP] == -1) {
317        t1c = (double *) t2c[0];
318        for (j = firstcol; j <= lastcol; j++) {
319            t1c[j] = 0.0;
320        }
321    }
322    if (gp[pid].neighbors[DOWN] == -1) {
323        t1c = (double *) t2c[im - 1];
324        for (j = firstcol; j <= lastcol; j++) {
325            t1c[j] = 0.0;
326        }
327    }
328    if (gp[pid].neighbors[LEFT] == -1) {
329        for (j = firstrow; j <= lastrow; j++) {
330            t2c[j][0] = 0.0;
331        }
332    }
333    if (gp[pid].neighbors[RIGHT] == -1) {
334        for (j = firstrow; j <= lastrow; j++) {
335            t2c[j][jm - 1] = 0.0;
336        }
337    }
338
339}
Note: See TracBrowser for help on using the repository browser.