1 /*
2  *  Copyright (C) 2003, Northwestern University and Argonne National Laboratory
3  *  See COPYRIGHT notice in top-level directory.
4  */
5 /* $Id: pres_temp_4D_wr.c 2744 2016-12-28 16:25:22Z wkliao $ */
6 
7 
8 /*
9    This is an example program which writes some 4D pressure and
10    temperatures. It is intended to illustrate the use of the netCDF
11    C API. The companion program pres_temp_4D_rd.c shows how
12    to read the netCDF data file created by this program.
13 
14    This program is part of the netCDF tutorial:
15    http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-tutorial
16 
17    Full documentation of the netCDF C API can be found at:
18    http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-c
19 
20    $Id: pres_temp_4D_wr.c 2744 2016-12-28 16:25:22Z wkliao $
21 */
22 
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <string.h>
26 #include <libgen.h> /* basename() */
27 #include <pnetcdf.h>
28 #include <mpi.h>
29 #include <testutils.h>
30 
31 /* This is the name of the data file we will create. */
32 #define FILE_NAME "pres_temp_4D.nc"
33 
34 /* We are writing 4D data, a 2 x 6 x 12 lvl-lat-lon grid, with 2
35    timesteps of data. */
36 #define NDIMS 4
37 #define NLAT 6
38 #define NLON 12
39 #define LAT_NAME "latitude"
40 #define LON_NAME "longitude"
41 #define NREC 2
42 #define REC_NAME "time"
43 #define LVL_NAME "level"
44 #define NLVL 4
45 
46 /* Names of things. */
47 #define PRES_NAME "pressure"
48 #define TEMP_NAME "temperature"
49 #define UNITS "units"
50 #define DEGREES_EAST "degrees_east"
51 #define DEGREES_NORTH "degrees_north"
52 
53 /* These are used to construct some example data. */
54 #define SAMPLE_PRESSURE 900
55 #define SAMPLE_TEMP 9.0
56 #define START_LAT 25.0
57 #define START_LON -125.0
58 
59 /* For the units attributes. */
60 #define UNITS "units"
61 #define PRES_UNITS "hPa"
62 #define TEMP_UNITS "celsius"
63 #define LAT_UNITS "degrees_north"
64 #define LON_UNITS "degrees_east"
65 #define MAX_ATT_LEN 80
66 
67 /* Handle errors by printing an error message and exiting with a
68  * non-zero status. */
69 #define ERR(e) {printf("Error: %s\n", nc_strerror(e)); return 2;}
70 
71 #define CHECK_ERR { \
72     if (err != NC_NOERR) { \
73         nerrs++; \
74         printf("Error: %s at line %d: %s\n", __FILE__,__LINE__,ncmpi_strerror(err)); \
75     } \
76 }
77 
78 
79 int
main(int argc,char ** argv)80 main(int argc, char ** argv)
81 {
82    /* IDs for the netCDF file, dimensions, and variables. */
83    int nprocs, rank, nerrs=0;
84    int ncid;
85    int lon_dimid, lat_dimid, lvl_dimid, rec_dimid;
86    int lat_varid, lon_varid, pres_varid, temp_varid;
87    int dimids[NDIMS];
88 
89    /* The start and count arrays will tell the netCDF library where to
90       write our data. */
91    MPI_Offset start[NDIMS], count[NDIMS];
92 
93    /* Program variables to hold the data we will write out. We will only
94       need enough space to hold one timestep of data; one record. */
95    float **pres_out; /* [NLVL/nprocs][NLAT][NLON] */
96    float **temp_out; /* [NLVL/nprocs][NLAT][NLON] */
97 
98    /* These program variables hold the latitudes and longitudes. */
99    float lats[NLAT], lons[NLON];
100 
101    /* Loop indexes. */
102    int lvl, lat, lon, rec, i = 0;
103 
104    /* Error handling. */
105    int err;
106 
107    char *filename=FILE_NAME;
108 
109    MPI_Init(&argc, &argv);
110    MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
111    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
112 
113    if (argc > 2) {
114        if (!rank) printf("Usage: %s [filename]\n",argv[0]);
115        MPI_Finalize();
116        return 0;
117    }
118    if (argc == 2) filename = argv[1];
119 
120    if (rank == 0) {
121        char *cmd_str = (char*)malloc(strlen(argv[0]) + 256);
122        sprintf(cmd_str, "*** TESTING C   %s for writing file", basename(argv[0]));
123        printf("%-66s ------ ", cmd_str);
124        free(cmd_str);
125    }
126 
127    /* Create some pretend data. If this wasn't an example program, we
128     * would have some real data to write, for example, model
129     * output. */
130    for (lat = 0; lat < NLAT; lat++)
131       lats[lat] = START_LAT + 5.*lat;
132    for (lon = 0; lon < NLON; lon++)
133       lons[lon] = START_LON + 5.*lon;
134 
135    /* Create the file. */
136    err = ncmpi_create(MPI_COMM_WORLD, filename, NC_CLOBBER, MPI_INFO_NULL, &ncid);
137 
138    CHECK_ERR
139 
140    /* Define the dimensions. The record dimension is defined to have
141     * unlimited length - it can grow as needed. In this example it is
142     * the time dimension.*/
143    err = ncmpi_def_dim(ncid, LVL_NAME, NLVL, &lvl_dimid);
144    CHECK_ERR
145    err = ncmpi_def_dim(ncid, LAT_NAME, NLAT, &lat_dimid);
146    CHECK_ERR
147    err = ncmpi_def_dim(ncid, LON_NAME, NLON, &lon_dimid);
148    CHECK_ERR
149    err = ncmpi_def_dim(ncid, REC_NAME, NC_UNLIMITED, &rec_dimid);
150    CHECK_ERR
151 
152    /* Define the coordinate variables. We will only define coordinate
153       variables for lat and lon.  Ordinarily we would need to provide
154       an array of dimension IDs for each variable's dimensions, but
155       since coordinate variables only have one dimension, we can
156       simply provide the address of that dimension ID (&lat_dimid) and
157       similarly for (&lon_dimid). */
158    err = ncmpi_def_var(ncid, LAT_NAME, NC_FLOAT, 1, &lat_dimid, &lat_varid);
159    CHECK_ERR
160    err = ncmpi_def_var(ncid, LON_NAME, NC_FLOAT, 1, &lon_dimid, &lon_varid);
161    CHECK_ERR
162 
163    /* Assign units attributes to coordinate variables. */
164    err = ncmpi_put_att_text(ncid, lat_varid, UNITS,
165 				 strlen(DEGREES_NORTH), DEGREES_NORTH);
166    CHECK_ERR
167    err = ncmpi_put_att_text(ncid, lon_varid, UNITS,
168 				 strlen(DEGREES_EAST), DEGREES_EAST);
169    CHECK_ERR
170 
171    /* The dimids array is used to pass the dimids of the dimensions of
172       the netCDF variables. Both of the netCDF variables we are
173       creating share the same four dimensions. In C, the
174       unlimited dimension must come first on the list of dimids. */
175    dimids[0] = rec_dimid;
176    dimids[1] = lvl_dimid;
177    dimids[2] = lat_dimid;
178    dimids[3] = lon_dimid;
179 
180    /* Define the netCDF variables for the pressure and temperature
181     * data. */
182    err = ncmpi_def_var(ncid, PRES_NAME, NC_FLOAT, NDIMS, dimids, &pres_varid);
183    CHECK_ERR
184    err = ncmpi_def_var(ncid, TEMP_NAME, NC_FLOAT, NDIMS, dimids, &temp_varid);
185    CHECK_ERR
186 
187    /* Assign units attributes to the netCDF variables. */
188    err = ncmpi_put_att_text(ncid, pres_varid, UNITS,
189 				 strlen(PRES_UNITS), PRES_UNITS);
190    CHECK_ERR
191    err = ncmpi_put_att_text(ncid, temp_varid, UNITS,
192 				 strlen(TEMP_UNITS), TEMP_UNITS);
193    CHECK_ERR
194 
195    /* End define mode. */
196    err = ncmpi_enddef(ncid);
197    CHECK_ERR
198 
199    err = ncmpi_begin_indep_data(ncid);
200    /* Write the coordinate variable data. This will put the latitudes
201       and longitudes of our data grid into the netCDF file. */
202    if (rank == 0) {
203        err = ncmpi_put_var_float(ncid, lat_varid, &lats[0]);
204        CHECK_ERR
205        err = ncmpi_put_var_float(ncid, lon_varid, &lons[0]);
206        CHECK_ERR
207    }
208    err = ncmpi_end_indep_data(ncid);
209    CHECK_ERR
210 
211    /* These settings tell netcdf to write one timestep of data. (The
212      setting of start[0] inside the loop below tells netCDF which
213                     &data[0][0][0]);
214      timestep to write.) */
215    count[0] = 1;
216    count[2] = NLAT;
217    count[3] = NLON;
218    start[2] = 0;
219    start[3] = 0;
220 
221    /* divide NLVL dimension among processes */
222    count[1] = NLVL / nprocs;
223    start[1] = count[1] * rank;
224    if (rank < NLVL % nprocs) {
225        start[1] += rank;
226        count[1]++;
227    }
228    else {
229        start[1] += NLVL % nprocs;
230    }
231    if (count[1] == 0) start[1] = 0;
232 
233    /* allocate write buffers */
234    pres_out = (float**) malloc(count[1]*2 * sizeof(float*));
235    temp_out = pres_out + count[1];
236    if (count[1] > 0) {
237        pres_out[0] = (float*) malloc(count[1]*2 * NLAT*NLON * sizeof(float));
238        temp_out[0] = pres_out[0] + count[1] * NLAT*NLON;
239        for (i=1; i<count[1]; i++) {
240            pres_out[i] = pres_out[i-1] + NLAT*NLON;
241            temp_out[i] = temp_out[i-1] + NLAT*NLON;
242        }
243    }
244 
245    /* initialize write buffers */
246    i = (int)start[1] * NLAT * NLON;
247    for (lvl=0; lvl<count[1]; lvl++)
248       for (lat = 0; lat < NLAT; lat++)
249 	 for (lon = 0; lon < NLON; lon++) {
250 	    pres_out[lvl][lat*NLON + lon] = SAMPLE_PRESSURE + i;
251 	    temp_out[lvl][lat*NLON + lon] = SAMPLE_TEMP + i++;
252 	 }
253 
254    /* Write the pretend data. This will write our surface pressure and
255       surface temperature data. The arrays only hold one timestep worth
256       of data. We will just rewrite the same data for each timestep. In
257       a real application, the data would change between timesteps. */
258 
259    for (rec = 0; rec < NREC; rec++)
260    {
261       start[0] = rec;
262       err = ncmpi_put_vara_float_all(ncid, pres_varid, start, count, &pres_out[0][0]);
263       CHECK_ERR
264       err = ncmpi_put_vara_float_all(ncid, temp_varid, start, count, &temp_out[0][0]);
265       CHECK_ERR
266    }
267 
268    /* Close the file. */
269    err = ncmpi_close(ncid);
270    CHECK_ERR
271 
272    if (count[1] > 0) free(pres_out[0]);
273    free(pres_out);
274 
275     /* check if there is any malloc residue */
276     MPI_Offset malloc_size, sum_size;
277     err = ncmpi_inq_malloc_size(&malloc_size);
278     if (err == NC_NOERR) {
279         MPI_Reduce(&malloc_size, &sum_size, 1, MPI_OFFSET, MPI_SUM, 0, MPI_COMM_WORLD);
280         if (rank == 0 && sum_size > 0)
281             printf("heap memory allocated by PnetCDF internally has %lld bytes yet to be freed\n",
282                    sum_size);
283     }
284 
285     MPI_Allreduce(MPI_IN_PLACE, &nerrs, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
286     if (rank == 0) {
287         if (nerrs) printf(FAIL_STR,nerrs);
288         else       printf(PASS_STR);
289     }
290 
291     MPI_Finalize();
292 
293     return nerrs;
294 }
295