source: CONFIG_DEVT/IPSLCM6.2.2_ENSEMBLES/oasis3-mct/pyoasis/examples/10-grid/C/writer.c @ 6204

Last change on this file since 6204 was 6204, checked in by aclsce, 2 years ago

First import of IPSLCM6.2.2_work_ENSEMBLES working configuration

File size: 13.7 KB
Line 
1#include <stdio.h>
2#include <math.h>
3#include "mpi.h"
4#include "netcdf.h"
5#include "oasis_c.h"
6
7double maxval(int m, int n, double a[m][n]) {
8  int c, d;
9  double maximum = a[0][0];
10
11  for (c = 0; c < m; c++)
12    for (d = 0; d < n; d++)
13      if (a[c][d] > maximum)
14        maximum = a[c][d];
15  return maximum;
16}
17
18double minval(int m, int n, double a[m][n]) {
19  int c, d;
20  double minimum = a[0][0];
21
22  for (c = 0; c < m; c++)
23    for (d = 0; d < n; d++)
24      if (a[c][d] < minimum)
25        minimum = a[c][d];
26  return minimum;
27}
28
29int main(int argc, char *argv[])
30{
31  char *comp_name = "writer";
32  fprintf(stdout,"Component name: %s\n", comp_name);
33  fflush(stdout);
34  int comp_id;
35  OASIS_CHECK_ERR(oasis_c_init_comp(&comp_id, comp_name, OASIS_COUPLED));
36
37  MPI_Comm localcomm;
38  OASIS_CHECK_ERR(oasis_c_get_localcomm(&localcomm));
39  int comm_size;
40  OASIS_CHECK_MPI_ERR(MPI_Comm_size(localcomm, &comm_size));
41  int comm_rank;
42  OASIS_CHECK_MPI_ERR(MPI_Comm_rank(localcomm, &comm_rank));
43  fprintf(stdout,"Writer: rank = %d of %d\n",comm_rank,comm_size);
44  fflush(stdout);
45
46  const int nx_loc = 18;
47  const int ny_loc = 18;
48  const int n_crn  =  4;
49  int nx_global, ny_global;
50  double dx,dy;
51  int i, j;
52  const double dp_conv = atan(1.)/45.0;
53
54  int part_params[OASIS_Box_Params];
55  part_params[OASIS_Strategy] = OASIS_Box;
56  part_params[OASIS_Offset]   = comm_rank * nx_loc;
57  part_params[OASIS_SizeX]    = nx_loc;
58  part_params[OASIS_SizeY]    = ny_loc;
59  part_params[OASIS_LdX]      = comm_size * nx_loc;
60  int part_id;
61  OASIS_CHECK_ERR(oasis_c_def_partition(&part_id, OASIS_Box_Params,
62                                        part_params, OASIS_No_Gsize,
63                                        OASIS_No_Name));
64
65  OASIS_CHECK_ERR(oasis_c_start_grids_writing());
66
67  {
68    /* pyoa */
69    nx_global = comm_size * nx_loc;
70    ny_global = ny_loc;
71    dx = 360.0/(double) nx_global;
72    dy = 180.0/(double) ny_global;
73
74    double lon[ny_loc][nx_loc];
75    double lat[ny_loc][nx_loc];
76    for ( i = 0; i<nx_loc; i++) {
77      for ( j = 0; j<ny_loc; j++) {
78        lon[j][i] = -180.0 + comm_rank*nx_loc*dx + (double) i * dx + dx/2. ;
79        lat[j][i] = -90.0 + (double) j * dy + dy/2.;
80      }
81    }
82    OASIS_CHECK_ERR(oasis_c_write_grid("pyoa", nx_global, ny_global,
83                                       nx_loc, ny_loc, &lon[0][0], &lat[0][0], part_id));
84
85    double clo[n_crn][ny_loc][nx_loc];
86    double cla[n_crn][ny_loc][nx_loc];
87    for ( i = 0; i<nx_loc; i++) {
88      for ( j = 0; j<ny_loc; j++) {
89        clo[3][j][i] = clo[0][j][i] = -180.0 + comm_rank*nx_loc*dx + (double) i * dx;
90        clo[2][j][i] = clo[1][j][i] = -180.0 + comm_rank*nx_loc*dx + (double) (i+1) * dx;
91        cla[1][j][i] = cla[0][j][i] = -90.0 + (double) j * dy;
92        cla[3][j][i] = cla[2][j][i] = -90.0 + (double) (j+1) * dy;
93      }
94    }
95    OASIS_CHECK_ERR(oasis_c_write_corner("pyoa", nx_global, ny_global, n_crn,
96                                         nx_loc, ny_loc, &clo[0][0][0], &cla[0][0][0], part_id));
97
98    int imsk[ny_loc][nx_loc];
99    for ( i = 0; i<nx_loc; i++) {
100      for ( j = 0; j<ny_loc; j++) {
101        imsk[j][i] = 0;
102      }
103    }
104    switch (comm_rank) {
105    case 0:
106      for (i = 4; i<6; i++) for (j = 2; j<16; j++) imsk[j][i] = 1;
107      for (i = 6; i<11; i++) imsk[15][i] = imsk[14][i] = imsk[9][i] = imsk[8][i] = 1;
108      imsk[15][11] = imsk[14][11] = imsk[13][11] = imsk[10][11] = imsk[9][11] = imsk[8][11] = 1;
109      for (j = 9; j<15; j++) imsk[j][12] = 1;
110      for (j = 10; j<14; j++) imsk[j][13] = 1;
111      break;
112    case 1:
113      for (j = 14; j<16; j++) imsk[j][3] = imsk[j][14] = 1;
114      for (j = 12; j<16; j++) imsk[j][4] = imsk[j][13] = 1;
115      for (j = 11; j<15; j++) imsk[j][5] = imsk[j][12] = 1;
116      for (j = 10; j<13; j++) imsk[j][6] = imsk[j][11] = 1;
117      for (j =  9; j<12; j++) imsk[j][7] = imsk[j][10] = 1;
118      for (i =  8; i<10; i++) for (j = 2; j<11; j++) imsk[j][i] = 1;
119      break;
120    case 2:
121      for (j =  4; j<14; j++) imsk[j][4] = imsk[j][13] = 1;
122      for (j =  3; j<15; j++) imsk[j][5] = imsk[j][12] = 1;
123      for (j =  2; j< 5; j++) imsk[j][6] = imsk[j][11] = 1;
124      for (j = 13; j<16; j++) imsk[j][6] = imsk[j][11] = 1;
125      for (i =  7; i<11; i++) for (j = 2; j<4; j++) imsk[j][i] = 1;
126      for (i =  7; i<11; i++) for (j = 14; j<16; j++) imsk[j][i] = 1;
127      break;
128    }
129    OASIS_CHECK_ERR(oasis_c_write_mask("pyoa", nx_global, ny_global,
130                                       nx_loc, ny_loc, &imsk[0][0], part_id, "STFC"));
131
132    double frac[ny_loc][nx_loc];
133    for ( i = 0; i<nx_loc; i++) {
134      for ( j = 0; j<ny_loc; j++) {
135        frac[j][i] = 1.0;
136        if (imsk[j][i] == 1) frac[j][i] = 0.;
137      }
138    }
139    OASIS_CHECK_ERR(oasis_c_write_frac("pyoa", nx_global, ny_global,
140                                       nx_loc, ny_loc, &frac[0][0], part_id, "STFC"));
141
142    double area[ny_loc][nx_loc];
143    for ( i = 0; i<nx_loc; i++) {
144      for ( j = 0; j<ny_loc; j++) {
145        area[j][i] = dp_conv * fabs(sin(cla[2][j][i]*dp_conv)-sin(cla[0][j][i]*dp_conv)) *
146          fabs(clo[1][j][i]-clo[0][j][i]);
147      }
148    }
149    OASIS_CHECK_ERR(oasis_c_write_area("pyoa", nx_global, ny_global,
150                                       nx_loc, ny_loc, &area[0][0], part_id));
151
152    double angle[ny_loc][nx_loc];
153    for ( i = 0; i<nx_loc; i++) for ( j = 0; j<ny_loc; j++) angle[j][i] = 0.;
154    OASIS_CHECK_ERR(oasis_c_write_angle("pyoa", nx_global, ny_global,
155                                        nx_loc, ny_loc, &angle[0][0], part_id));
156
157  }
158
159  {
160    /* mono */
161    nx_global = 180;
162    ny_global = 90;
163    dx = 360.0/(double) nx_global;
164    dy = 180.0/(double) ny_global;
165
166    double lon[ny_global][nx_global];
167    double lat[ny_global][nx_global];
168    for ( i = 0; i<nx_global; i++) {
169      for ( j = 0; j<ny_global; j++) {
170        lon[j][i] = (double) i * dx + dx/2. ;
171        lat[j][i] = -90.0 + (double) j * dy + dy/2.;
172      }
173    }
174    OASIS_CHECK_ERR(oasis_c_write_grid("mono", nx_global, ny_global,
175                                       nx_global, ny_global, &lon[0][0], &lat[0][0],
176                                       OASIS_No_Part));
177
178    double clo[n_crn][ny_global][nx_global];
179    double cla[n_crn][ny_global][nx_global];
180    for ( i = 0; i<nx_global; i++) {
181      for ( j = 0; j<ny_global; j++) {
182        clo[3][j][i] = clo[0][j][i] = (double) i * dx;
183        clo[2][j][i] = clo[1][j][i] = (double) (i+1) * dx;
184        cla[1][j][i] = cla[0][j][i] = -90.0 + (double) j * dy;
185        cla[3][j][i] = cla[2][j][i] = -90.0 + (double) (j+1) * dy;
186      }
187    }
188    OASIS_CHECK_ERR(oasis_c_write_corner("mono", nx_global, ny_global, n_crn,
189                                         nx_global, ny_global, &clo[0][0][0], &cla[0][0][0],
190                                         OASIS_No_Part));
191
192    int imsk[ny_global][nx_global];
193    for ( i = 0; i<nx_global; i++) {
194      for ( j = 0; j<ny_global; j++) {
195        imsk[j][i] = 0;
196        if ((lon[j][i]-180.)*(lon[j][i]-180.)+(lat[j][i]*lat[j][i]) < 900.) imsk[j][i] = 1;
197      }
198    }
199    OASIS_CHECK_ERR(oasis_c_write_mask("mono", nx_global, ny_global,
200                                       nx_global, ny_global, &imsk[0][0], OASIS_No_Part,
201                                       OASIS_No_Companion));
202    double frac[ny_global][nx_global];
203    for ( i = 0; i<nx_global; i++) {
204      for ( j = 0; j<ny_global; j++) {
205        frac[j][i] = 1.0;
206        if (imsk[j][i] == 1) frac[j][i] = 0.;
207      }
208    }
209    OASIS_CHECK_ERR(oasis_c_write_frac("mono", nx_global, ny_global,
210                                       nx_global, ny_global, &frac[0][0], OASIS_No_Part,
211                                       OASIS_No_Companion));
212
213    double area[ny_global][nx_global];
214    for ( i = 0; i<nx_global; i++) {
215      for ( j = 0; j<ny_global; j++) {
216        area[j][i] = dp_conv * fabs(sin(cla[2][j][i]*dp_conv)-sin(cla[0][j][i]*dp_conv)) *
217          fabs(clo[1][j][i]-clo[0][j][i]);
218      }
219    }
220    OASIS_CHECK_ERR(oasis_c_write_area("mono", nx_global, ny_global,
221                                       nx_global, ny_global, &area[0][0], OASIS_No_Part));
222
223    double angle[ny_global][nx_global];
224    for ( i = 0; i<nx_global; i++) for ( j = 0; j<ny_global; j++) angle[j][i] = 0.;
225    OASIS_CHECK_ERR(oasis_c_write_angle("mono", nx_global, ny_global,
226                                        nx_global, ny_global, &angle[0][0], OASIS_No_Part));
227  }
228
229  OASIS_CHECK_ERR(oasis_c_terminate_grids_writing());
230
231  OASIS_CHECK_ERR(oasis_c_enddef());
232
233  if ( comm_rank == 0 ) {
234    int status;
235    int ncid, varid;
236    nx_global = 180;
237    ny_global = 90;
238
239
240    /* Open grids.nc check for min and max of lat and lon */
241    if (nc_open("grids.nc", 0, &ncid) != NC_NOERR)
242      OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in opening grids.nc",
243                                    __FILE__, __LINE__));
244
245    {
246      double lon[ny_loc][comm_size*nx_loc];
247      double lat[ny_loc][comm_size*nx_loc];
248      if (nc_inq_varid(ncid, "pyoa.lat", &varid) != NC_NOERR)
249        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in getting pyoa.lat id",__FILE__, __LINE__));
250      if (nc_get_var(ncid, varid, &lat[0][0]) != NC_NOERR)
251        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in getting pyoa.lat values",__FILE__, __LINE__));
252      if (maxval(ny_loc,comm_size*nx_loc,lat) > 90.0)
253        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in pyoa.lat max value",__FILE__, __LINE__));
254      if (minval(ny_loc,comm_size*nx_loc,lat) < -90.0)
255        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in pyoa.lat min value",__FILE__, __LINE__));
256      if (nc_inq_varid(ncid, "pyoa.lon", &varid) != NC_NOERR)
257        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in getting pyoa.lat id",__FILE__, __LINE__));
258      if (nc_get_var(ncid, varid, &lon[0][0]) != NC_NOERR)
259        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in getting pyoa.lat values",__FILE__, __LINE__));
260      if (maxval(ny_loc,comm_size*nx_loc,lon) > 180.0)
261        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in pyoa.lon max value",__FILE__, __LINE__));
262      if (minval(ny_loc,comm_size*nx_loc,lon) < -180.0)
263        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in pyoa.lon min value",__FILE__, __LINE__));
264
265    }
266    {
267      double lon[ny_global][nx_global];
268      double lat[ny_global][nx_global];
269      if (nc_inq_varid(ncid, "mono.lat", &varid) != NC_NOERR)
270        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in getting mono.lat id",__FILE__, __LINE__));
271      if (nc_get_var(ncid, varid, &lat[0][0]) != NC_NOERR)
272        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in getting mono.lat values",__FILE__, __LINE__));
273      if (maxval(ny_global,nx_global,lat) > 90.0)
274        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in mono.lat max value",__FILE__, __LINE__));
275      if (minval(ny_global,nx_global,lat) < -90.0)
276        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in mono.lat min value",__FILE__, __LINE__));
277      if (nc_inq_varid(ncid, "mono.lon", &varid) != NC_NOERR)
278        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in getting mono.lat id",__FILE__, __LINE__));
279      if (nc_get_var(ncid, varid, &lon[0][0]) != NC_NOERR)
280        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in getting mono.lat values",__FILE__, __LINE__));
281      if (maxval(ny_global,nx_global,lon) > 360.0)
282        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in mono.lon max value",__FILE__, __LINE__));
283      if (minval(ny_global,nx_global,lon) < 0.0)
284        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in mono.lon min value",__FILE__, __LINE__));
285
286    }
287    if (nc_close(ncid))
288      OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in closing grids.nc",__FILE__, __LINE__));
289
290    /* Open masks.nc and check for the hole in mono.msk */
291    if (nc_open("masks.nc", 0, &ncid) != NC_NOERR)
292      OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in opening masks.nc",__FILE__, __LINE__));
293
294    {
295      int imsk[ny_global][nx_global];
296      for ( i = 0; i<nx_global; i++) for ( j = 0; j<ny_global; j++) imsk[j][i] = 0;
297      if (nc_inq_varid(ncid, "mono.msk", &varid) != NC_NOERR)
298        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in getting mono.msk id",__FILE__, __LINE__));
299      if (nc_get_var(ncid, varid, &imsk[0][0]) != NC_NOERR)
300        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in getting mono.msk values",__FILE__, __LINE__));
301      double dist;
302      for ( i = 0; i<nx_global; i++)
303        for ( j = 0; j<ny_global; j++) {
304          dist = sqrt(pow((double)(i - nx_global/2),2.0) + pow((double) (j - ny_global/2),2.0));
305          if ( dist < 14.0 ) {
306            if ( imsk[j][i] != 1 )
307              OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in mono.msk value",__FILE__, __LINE__));
308          } else if ( dist > 16.0 ) {
309            if ( imsk[j][i] != 0 )
310              OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in mono.msk value",__FILE__, __LINE__));
311          }
312        }
313    }
314    if (nc_close(ncid))
315      OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in closing masks.nc",__FILE__, __LINE__));
316
317    /* Open areas.nc and check for the area on the sphere */
318    if (nc_open("areas.nc", 0, &ncid) != NC_NOERR)
319      OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in opening areas.nc",__FILE__, __LINE__));
320
321    {
322      double area[ny_loc][comm_size*nx_loc];
323      for ( i = 0; i<comm_size*nx_loc; i++) for ( j = 0; j<ny_loc; j++) area[j][i] = 0.0;
324      if (nc_inq_varid(ncid, "pyoa.srf", &varid) != NC_NOERR)
325        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in getting pyoa.srf id",__FILE__, __LINE__));
326      if (nc_get_var(ncid, varid, &area[0][0]) != NC_NOERR)
327        OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in getting pyoa.srf values",__FILE__, __LINE__));
328
329      const double epsilon = 1.e-6;
330      for ( i = 0; i<comm_size*nx_loc; i++)
331        for ( j = 0; j<ny_loc; j++) {
332          if (fabs(area[j][i] - dp_conv *
333                   fabs(sin((-90.+(double)(j+1) * 180./ny_loc)*dp_conv)-
334                        sin((-90.+(double)(j)   * 180./ny_loc)*dp_conv)) *
335                   fabs(360./(double)(comm_size*nx_loc))) > epsilon)
336            OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in pyoa.srf value",__FILE__, __LINE__));
337        }
338    }
339    if (nc_close(ncid))
340      OASIS_CHECK_ERR(oasis_c_abort(comp_id, comp_name, "Writer: Error in closing areas.nc",__FILE__, __LINE__));
341
342    /* Write out the successfully message if everything went well */
343    fprintf(stdout,"Writer: successfully ended writing grids\n");
344    fflush(stdout);
345  }
346
347  OASIS_CHECK_ERR(oasis_c_terminate());
348
349}
Note: See TracBrowser for help on using the repository browser.