source: soft/giet_vm/applications/ocean/jacobcalc2.c @ 807

Last change on this file since 807 was 799, checked in by alain, 9 years ago

Introducing a port of ocean application without macros.
File size: 11.1 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
20#include "decs.h"
21
22////////////////////////////
23void jacobcalc2( double ****x, 
24                 double ****y, 
25                 double ****z, 
26                 long psiindex, 
27                 long procid, 
28                 long firstrow, 
29                 long lastrow, 
30                 long firstcol, 
31                 long lastcol )
32{
33    double f1;
34    double f2;
35    double f3;
36    double f4;
37    double f5;
38    double f6;
39    double f7;
40    double f8;
41    long iindex;
42    long indexp1;
43    long indexm1;
44    long im1;
45    long ip1;
46    long i;
47    long j;
48    long jj;
49    double **t2a;
50    double **t2b;
51    double **t2c;
52    double *t1a;
53    double *t1b;
54    double *t1c;
55    double *t1d;
56    double *t1e;
57    double *t1f;
58    double *t1g;
59
60    t2a = z[procid][psiindex];
61    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) 
62    {
63        t2a[0][0] = 0.0;
64    }
65    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) 
66    {
67        t2a[im - 1][0] = 0.0;
68    }
69    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) 
70    {
71        t2a[0][jm - 1] = 0.0;
72    }
73    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) 
74    {
75        t2a[im - 1][jm - 1] = 0.0;
76    }
77
78    t2a = x[procid][psiindex];
79    jj = gps[procid]->neighbors[UPLEFT];
80    if (jj != -1) {
81        t2a[0][0] = x[jj][psiindex][im - 2][jm - 2];
82    }
83    jj = gps[procid]->neighbors[UPRIGHT];
84    if (jj != -1) {
85        t2a[0][jm - 1] = x[jj][psiindex][im - 2][1];
86    }
87    jj = gps[procid]->neighbors[DOWNLEFT];
88    if (jj != -1) {
89        t2a[im - 1][0] = x[jj][psiindex][1][jm - 2];
90    }
91    jj = gps[procid]->neighbors[DOWNRIGHT];
92    if (jj != -1) {
93        t2a[im - 1][jm - 1] = x[jj][psiindex][1][1];
94    }
95
96    t2a = y[procid][psiindex];
97    jj = gps[procid]->neighbors[UPLEFT];
98    if (jj != -1) {
99        t2a[0][0] = y[jj][psiindex][im - 2][jm - 2];
100    }
101    jj = gps[procid]->neighbors[UPRIGHT];
102    if (jj != -1) {
103        t2a[0][jm - 1] = y[jj][psiindex][im - 2][1];
104    }
105    jj = gps[procid]->neighbors[DOWNLEFT];
106    if (jj != -1) {
107        t2a[im - 1][0] = y[jj][psiindex][1][jm - 2];
108    }
109    jj = gps[procid]->neighbors[DOWNRIGHT];
110    if (jj != -1) {
111        t2a[im - 1][jm - 1] = y[jj][psiindex][1][1];
112    }
113
114    t2a = x[procid][psiindex];
115    if (gps[procid]->neighbors[UP] == -1) {
116        jj = gps[procid]->neighbors[LEFT];
117        if (jj != -1) {
118            t2a[0][0] = x[jj][psiindex][0][jm - 2];
119        } else {
120            jj = gps[procid]->neighbors[DOWN];
121            if (jj != -1) {
122                t2a[im - 1][0] = x[jj][psiindex][1][0];
123            }
124        }
125        jj = gps[procid]->neighbors[RIGHT];
126        if (jj != -1) {
127            t2a[0][jm - 1] = x[jj][psiindex][0][1];
128        } else {
129            jj = gps[procid]->neighbors[DOWN];
130            if (jj != -1) {
131                t2a[im - 1][jm - 1] = x[jj][psiindex][1][jm - 1];
132            }
133        }
134    } else if (gps[procid]->neighbors[DOWN] == -1) {
135        jj = gps[procid]->neighbors[LEFT];
136        if (jj != -1) {
137            t2a[im - 1][0] = x[jj][psiindex][im - 1][jm - 2];
138        } else {
139            jj = gps[procid]->neighbors[UP];
140            if (jj != -1) {
141                t2a[0][0] = x[jj][psiindex][im - 2][0];
142            }
143        }
144        jj = gps[procid]->neighbors[RIGHT];
145        if (jj != -1) {
146            t2a[im - 1][jm - 1] = x[jj][psiindex][im - 1][1];
147        } else {
148            jj = gps[procid]->neighbors[UP];
149            if (jj != -1) {
150                t2a[0][jm - 1] = x[jj][psiindex][im - 2][jm - 1];
151            }
152        }
153    } else if (gps[procid]->neighbors[LEFT] == -1) {
154        jj = gps[procid]->neighbors[UP];
155        if (jj != -1) {
156            t2a[0][0] = x[jj][psiindex][im - 2][0];
157        }
158        jj = gps[procid]->neighbors[DOWN];
159        if (jj != -1) {
160            t2a[im - 1][0] = x[jj][psiindex][1][0];
161        }
162    } else if (gps[procid]->neighbors[RIGHT] == -1) {
163        jj = gps[procid]->neighbors[UP];
164        if (jj != -1) {
165            t2a[0][jm - 1] = x[jj][psiindex][im - 2][jm - 1];
166        }
167        jj = gps[procid]->neighbors[DOWN];
168        if (jj != -1) {
169            t2a[im - 1][jm - 1] = x[jj][psiindex][1][jm - 1];
170        }
171    }
172
173    t2a = y[procid][psiindex];
174    if (gps[procid]->neighbors[UP] == -1) {
175        jj = gps[procid]->neighbors[LEFT];
176        if (jj != -1) {
177            t2a[0][0] = y[jj][psiindex][0][jm - 2];
178        } else {
179            jj = gps[procid]->neighbors[DOWN];
180            if (jj != -1) {
181                t2a[im - 1][0] = y[jj][psiindex][1][0];
182            }
183        }
184        jj = gps[procid]->neighbors[RIGHT];
185        if (jj != -1) {
186            t2a[0][jm - 1] = y[jj][psiindex][0][1];
187        } else {
188            jj = gps[procid]->neighbors[DOWN];
189            if (jj != -1) {
190                t2a[im - 1][jm - 1] = y[jj][psiindex][1][jm - 1];
191            }
192        }
193    } else if (gps[procid]->neighbors[DOWN] == -1) {
194        jj = gps[procid]->neighbors[LEFT];
195        if (jj != -1) {
196            t2a[im - 1][0] = y[jj][psiindex][im - 1][jm - 2];
197        } else {
198            jj = gps[procid]->neighbors[UP];
199            if (jj != -1) {
200                t2a[0][0] = y[jj][psiindex][im - 2][0];
201            }
202        }
203        jj = gps[procid]->neighbors[RIGHT];
204        if (jj != -1) {
205            t2a[im - 1][jm - 1] = y[jj][psiindex][im - 1][1];
206        } else {
207            jj = gps[procid]->neighbors[UP];
208            if (jj != -1) {
209                t2a[0][jm - 1] = y[jj][psiindex][im - 2][jm - 1];
210            }
211        }
212    } else if (gps[procid]->neighbors[LEFT] == -1) {
213        jj = gps[procid]->neighbors[UP];
214        if (jj != -1) {
215            t2a[0][0] = y[jj][psiindex][im - 2][0];
216        }
217        jj = gps[procid]->neighbors[DOWN];
218        if (jj != -1) {
219            t2a[im - 1][0] = y[jj][psiindex][1][0];
220        }
221    } else if (gps[procid]->neighbors[RIGHT] == -1) {
222        jj = gps[procid]->neighbors[UP];
223        if (jj != -1) {
224            t2a[0][jm - 1] = y[jj][psiindex][im - 2][jm - 1];
225        }
226        jj = gps[procid]->neighbors[DOWN];
227        if (jj != -1) {
228            t2a[im - 1][jm - 1] = y[jj][psiindex][1][jm - 1];
229        }
230    }
231
232    t2a = y[procid][psiindex];
233    j = gps[procid]->neighbors[UP];
234    if (j != -1) {
235        t1a = (double *) t2a[0];
236        t1b = (double *) y[j][psiindex][im - 2];
237        for (i = 1; i <= lastcol; i++) {
238            t1a[i] = t1b[i];
239        }
240    }
241    j = gps[procid]->neighbors[DOWN];
242    if (j != -1) {
243        t1a = (double *) t2a[im - 1];
244        t1b = (double *) y[j][psiindex][1];
245        for (i = 1; i <= lastcol; i++) {
246            t1a[i] = t1b[i];
247        }
248    }
249    j = gps[procid]->neighbors[LEFT];
250    if (j != -1) {
251        t2b = y[j][psiindex];
252        for (i = 1; i <= lastrow; i++) {
253            t2a[i][0] = t2b[i][jm - 2];
254        }
255    }
256    j = gps[procid]->neighbors[RIGHT];
257    if (j != -1) {
258        t2b = y[j][psiindex];
259        for (i = 1; i <= lastrow; i++) {
260            t2a[i][jm - 1] = t2b[i][1];
261        }
262    }
263
264    t2a = x[procid][psiindex];
265    j = gps[procid]->neighbors[UP];
266    if (j != -1) {
267        t1a = (double *) t2a[0];
268        t1b = (double *) x[j][psiindex][im - 2];
269        for (i = 1; i <= lastcol; i++) {
270            t1a[i] = t1b[i];
271        }
272    }
273    j = gps[procid]->neighbors[DOWN];
274    if (j != -1) {
275        t1a = (double *) t2a[im - 1];
276        t1b = (double *) x[j][psiindex][1];
277        for (i = 1; i <= lastcol; i++) {
278            t1a[i] = t1b[i];
279        }
280    }
281    j = gps[procid]->neighbors[LEFT];
282    if (j != -1) {
283        t2b = x[j][psiindex];
284        for (i = 1; i <= lastrow; i++) {
285            t2a[i][0] = t2b[i][jm - 2];
286        }
287    }
288    j = gps[procid]->neighbors[RIGHT];
289    if (j != -1) {
290        t2b = x[j][psiindex];
291        for (i = 1; i <= lastrow; i++) {
292            t2a[i][jm - 1] = t2b[i][1];
293        }
294    }
295
296    t2a = x[procid][psiindex];
297    t2b = y[procid][psiindex];
298    t2c = z[procid][psiindex];
299    for (i = firstrow; i <= lastrow; i++) {
300        ip1 = i + 1;
301        im1 = i - 1;
302        t1a = (double *) t2a[i];
303        t1b = (double *) t2b[i];
304        t1c = (double *) t2c[i];
305        t1d = (double *) t2b[ip1];
306        t1e = (double *) t2b[im1];
307        t1f = (double *) t2a[ip1];
308        t1g = (double *) t2a[im1];
309        for (iindex = firstcol; iindex <= lastcol; iindex++) {
310            indexp1 = iindex + 1;
311            indexm1 = iindex - 1;
312            f1 = (t1b[indexm1] + t1d[indexm1] - t1b[indexp1] - t1d[indexp1]) * (t1f[iindex] - t1a[iindex]);
313            f2 = (t1e[indexm1] + t1b[indexm1] - t1e[indexp1] - t1b[indexp1]) * (t1a[iindex] - t1g[iindex]);
314            f3 = (t1d[iindex] + t1d[indexp1] - t1e[iindex] - t1e[indexp1]) * (t1a[indexp1] - t1a[iindex]);
315            f4 = (t1d[indexm1] + t1d[iindex] - t1e[indexm1] - t1e[iindex]) * (t1a[iindex] - t1a[indexm1]);
316            f5 = (t1d[iindex] - t1b[indexp1]) * (t1f[indexp1] - t1a[iindex]);
317            f6 = (t1b[indexm1] - t1e[iindex]) * (t1a[iindex] - t1g[indexm1]);
318            f7 = (t1b[indexp1] - t1e[iindex]) * (t1g[indexp1] - t1a[iindex]);
319            f8 = (t1d[iindex] - t1b[indexm1]) * (t1a[iindex] - t1f[indexm1]);
320
321            t1c[iindex] = factjacob * (f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8);
322
323        }
324    }
325
326    if (gps[procid]->neighbors[UP] == -1) 
327    {
328        t1c = (double *) t2c[0];
329        for (j = firstcol; j <= lastcol; j++) 
330        {
331            t1c[j] = 0.0;
332        }
333    }
334    if (gps[procid]->neighbors[DOWN] == -1) 
335    {
336        t1c = (double *) t2c[im - 1];
337        for (j = firstcol; j <= lastcol; j++) 
338        {
339            t1c[j] = 0.0;
340        }
341    }
342    if (gps[procid]->neighbors[LEFT] == -1) 
343    {
344        for (j = firstrow; j <= lastrow; j++) 
345        {
346            t2c[j][0] = 0.0;
347        }
348    }
349    if (gps[procid]->neighbors[RIGHT] == -1) 
350    {
351        for (j = firstrow; j <= lastrow; j++) 
352        {
353            t2c[j][jm - 1] = 0.0;
354        }
355    }
356} // end jacobclc2()
Note: See TracBrowser for help on using the repository browser.