c-lamp_postproc: add_total_fields.c@d3122367777b

     1 #include <stdio.h>

     2 #include <stdlib.h>

     3 #include <string.h>

     4 #include <math.h>

     5 #include "netcdf.h"

     7 /*

     8  * Loop through all history tapes from the Community Land Model adding new

     9  * fields that are totals of the multi-level fields included in total_vars[]

    10  * below.

    11  */

    13 static char *total_vars[] = {

    14 	"SOILICE",

    15 	"SOILLIQ"

    16 };

    17 static int nmtotal_vars = sizeof(total_vars)/sizeof(*total_vars);

    18 static char *combo_vars[] = {

    19 	"LATENT",		/* FCTR + FCEV + FGEV */

    20 	"ET",			/* QVEGE + QVEGT + QSOIL */

    21 	"ALL_SKY_ALBEDO",	/* FSR / FSDS */

    22 	"BLACK_SKY_ALBEDO",	/* (FSRNDLN + FSRVDLN)/(FSDSNDLN + FSDSVDLN) */

    23 	"NETRAD",		/* FSA - FIRA */

    24 };

    25 static char *combo_long_names[] = {

    26 	"latent heat flux",

    27 	"evapotranspiration",

    28 	"surface all-sky albedo",

    29 	"surface black-sky albedo",

    30 	"net radiation"

    31 };

    32 static char *combo_units[] = {

    33 	"watt/m^2",

    34 	"m/s",

    35 	"unitless",

    36 	"unitless",

    37 	"watt/m^2"

    38 };

    39 static float combo_FillValues[] = {

    40 	1.e+36,

    41 	1.e+36,

    42 	1.e+36,

    43 	1.e+36,

    44 	1.e+36

    45 };

    46 static int nmcombo_vars = sizeof(combo_vars)/sizeof(*combo_vars);

    47 static char *varname_total_prefix = "TOTAL_", *longname_total_prefix = "total ";

    48 static char *time_name = "time", *lon_name = "lon", *lat_name = "lat";

    49 static char *long_name_name = "long_name", *units_name = "units",

    50 	*cell_method_name = "cell_method", *FillValue_name = "_FillValue",

    51 	*missing_value_name = "missing_value";

    52 static char *cell_method = "time: mean";

    54 void wrap_nc(int status)

    55 {

    56 	if (status != NC_NOERR) {

    57 		fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));

    58 		exit(1);

    59 	}

    61 	return;

    62 }

    64 void read_total_timeslice(int ncid, int varid, size_t nlevels, size_t d1,

    65 	size_t d2, char FillFlag, float FillValue, size_t tslice, float *var,

    66 	double *tvar)

    67 {

    68 	int i, j;

    69 	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];

    71 	/* Initialize tvar to FillValue */

    72 	for (i = 0; i < (d1 * d2); i++)

    73 		tvar[i] = (double)FillValue;

    76 	start[0] = tslice;

    77 	count[0] = 1;

    78 	start[2] = 0;

    79 	count[2] = d1;

    80 	start[3] = 0;

    81 	count[3] = d2;

    83 	/* March through the levels, totalling as we go */

    84 	for (j = 0; j < nlevels; j++) {

    85 		start[1] = j;

    86 		count[1] = 1;

    87 		wrap_nc(nc_get_vara_float(ncid, varid, start, count, var));

    88 		for (i = 0; i < (d1 * d2); i++) {

    89 			if (tvar[i] == (double)FillValue)

    90 				tvar[i] = (double)var[i];

    91 			else

    92 				tvar[i] += (double)var[i];

    93 		}

    94 	}

    96 	return;

    97 }

    99 void alloc_svars(size_t nsvars, size_t d1, size_t d2, float ***svarsp)

   100 {

   101 	int i;

   102 	float **svars;

   104 	if (!(*svarsp = (float **)malloc(nsvars * sizeof(float **)))) {

   105 		perror("alloc_svars:svarsp");

   106 		exit(2);

   107 	}

   108 	svars = *svarsp;

   109 	for (i = 0; i < nsvars; i++) {

   110 		if (!(svars[i] = (float *)malloc(d1 * d2 * sizeof(float)))) {

   111 			perror("alloc_svars:svars[i]");

   112 			exit(2);

   113 		}

   114 	}

   115 	return;

   116 }

   118 void free_svars(size_t nsvars, float **svars)

   119 {

   120 	int i;

   122 	for (i = 0; i < nsvars; i++)

   123 		free(svars[i]);

   124 	free(svars);

   125 	return;

   126 }

   128 void read_timeslice(int ncid, size_t d1, size_t d2, size_t tslice, char *name,

   129 	float *var)

   130 {

   131 	int varid;

   132 	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];

   134 	start[0] = tslice;

   135 	count[0] = 1;

   136 	start[1] = 0;

   137 	count[1] = d1;

   138 	start[2] = 0;

   139 	count[2] = d2;

   141 	/* Get variable id */

   142 	wrap_nc(nc_inq_varid(ncid, name, &varid));

   143 	/* Assume the correct size and read the variable */

   144 	wrap_nc(nc_get_vara_float(ncid, varid, start, count, var));

   146 	return;

   147 }

   149 void read_compute_timeslice(int ncid, size_t d1, size_t d2, size_t tslice,

   150 	int num, double *tvar)

   151 {

   152 	int i;

   153 	size_t nsvars;

   154 	float **svars; /* source variables */

   155 	double denom;

   157 	/* Initialize tvar to FillValue */

   158 	for (i = 0; i < (d1 * d2); i++)

   159 		tvar[i] = (double)combo_FillValues[num];

   161 	switch (num) {

   162 		case 0: /* LATENT */

   163 			nsvars = 3;

   164 			alloc_svars(nsvars, d1, d2, &svars);

   165 			/* read timeslices */

   166 			read_timeslice(ncid, d1, d2, tslice, "FCTR", svars[0]);

   167 			read_timeslice(ncid, d1, d2, tslice, "FCEV", svars[1]);

   168 			read_timeslice(ncid, d1, d2, tslice, "FGEV", svars[2]);

   169 			/* compute new variable */

   170 			for (i = 0; i < (d1 * d2); i++) {

   171 				if (svars[0][i] != combo_FillValues[num] &&

   172 					svars[1][i] != combo_FillValues[num] &&

   173 					svars[2][i] != combo_FillValues[num]) {

   174 					tvar[i] = (double)svars[0][i]

   175 						+ (double)svars[1][i]

   176 						+ (double)svars[2][i];

   177 				}

   178 			}

   179 			free_svars(nsvars, svars);

   180 			break;

   181 		case 1: /* ET */

   182 			nsvars = 3;

   183 			alloc_svars(nsvars, d1, d2, &svars);

   184 			/* read timeslices */

   185 			read_timeslice(ncid, d1, d2, tslice, "QVEGE", svars[0]);

   186 			read_timeslice(ncid, d1, d2, tslice, "QVEGT", svars[1]);

   187 			read_timeslice(ncid, d1, d2, tslice, "QSOIL", svars[2]);

   188 			/* compute new variable */

   189 			for (i = 0; i < (d1 * d2); i++) {

   190 				if (svars[0][i] != combo_FillValues[num] &&

   191 					svars[1][i] != combo_FillValues[num] &&

   192 					svars[2][i] != combo_FillValues[num]) {

   193 					tvar[i] = (double)svars[0][i]

   194 						+ (double)svars[1][i]

   195 						+ (double)svars[2][i];

   196 				}

   197 			}

   198 			free_svars(nsvars, svars);

   199 			break;

   200 		case 2: /* ALL_SKY_ALBEDO */

   201 			nsvars = 2;

   202 			alloc_svars(nsvars, d1, d2, &svars);

   203 			/* read timeslices */

   204 			read_timeslice(ncid, d1, d2, tslice, "FSR", svars[0]);

   205 			read_timeslice(ncid, d1, d2, tslice, "FSDS", svars[1]);

   206 			/* compute new variable */

   207 			for (i = 0; i < (d1 * d2); i++) {

   208 				if (svars[0][i] != combo_FillValues[num] &&

   209 					svars[1][i] != combo_FillValues[num] &&

   210 					(svars[1][i] > 1.e-35 || svars[1][i] < -1.e-35)) {

   211 					tvar[i] = (double)svars[0][i]

   212 						/ (double)svars[1][i];

   213 				}

   214 			}

   215 			free_svars(nsvars, svars);

   216 			break;

   217 		case 3: /* BLACK_SKY_ALBEDO */

   218 			nsvars = 4;

   219 			alloc_svars(nsvars, d1, d2, &svars);

   220 			/* read timeslices */

   221 			read_timeslice(ncid, d1, d2, tslice, "FSRNDLN", svars[0]);

   222 			read_timeslice(ncid, d1, d2, tslice, "FSRVDLN", svars[1]);

   223 			read_timeslice(ncid, d1, d2, tslice, "FSDSNDLN", svars[2]);

   224 			read_timeslice(ncid, d1, d2, tslice, "FSDSVDLN", svars[3]);

   225 			/* compute new variable */

   226 			for (i = 0; i < (d1 * d2); i++) {

   227 				if (svars[0][i] != combo_FillValues[num] &&

   228 					svars[1][i] != combo_FillValues[num] &&

   229 					svars[2][i] != combo_FillValues[num] &&

   230 					svars[3][i] != combo_FillValues[num]) {

   231 					denom = (double)svars[2][i] + (double)svars[3][i];

   232 					if (denom > 1.e-35 || denom < -1.e-35) {

   233 						tvar[i] = ((double)svars[0][i]

   234 							+ (double)svars[1][i]) /

   235 							((double)svars[2][i]

   236 							+ (double)svars[3][i]);

   237 					}

   238 				}

   239 			}

   240 			free_svars(nsvars, svars);

   241 			break;

   242 		case 4: /* NETRAD */

   243 			nsvars = 2;

   244 			alloc_svars(nsvars, d1, d2, &svars);

   245 			/* read timeslices */

   246 			read_timeslice(ncid, d1, d2, tslice, "FSA", svars[0]);

   247 			read_timeslice(ncid, d1, d2, tslice, "FIRA", svars[1]);

   248 			/* compute new variable */

   249 			for (i = 0; i < (d1 * d2); i++) {

   250 				if (svars[0][i] != combo_FillValues[num] &&

   251 					svars[1][i] != combo_FillValues[num]) {

   252 					tvar[i] = (double)svars[0][i]

   253 						- (double)svars[1][i];

   254 				}

   255 			}

   256 			free_svars(nsvars, svars);

   257 			break;

   258 		default:

   259 			fprintf(stderr, "Unknown combination variable %d!\n", num);

   260 			exit(3);

   261 	}

   263 	return;

   264 }

   266 void write_timeslice(int ncid, int varid, size_t d1, size_t d2,

   267 	size_t tslice, float *var)

   268 {

   269 	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];

   271 	start[0] = tslice;

   272 	count[0] = 1;

   273 	start[1] = 0;

   274 	count[1] = d1;

   275 	start[2] = 0;

   276 	count[2] = d2;

   278 	wrap_nc(nc_put_vara_float(ncid, varid, start, count, var));

   280 	return;

   281 }

   283 void usage(char *program)

   284 {

   285 	fprintf(stderr, "Usage: %s hist_file1.nc [ hist_file2.nc ... ]\n",

   286 		program);

   287 	fprintf(stderr, "Requires at least one history output file name\n");

   288 	exit(5);

   289 }

   291 int main(int argc, char **argv)

   292 {

   293 	char **ifnames;

   294 	int i, j, k, nifnames;

   295 	int ncid, ngdims, nvars, ngatts, unlimdimid, varid, ndims, natts,

   296 		dimids[NC_MAX_VAR_DIMS], new_ndims, new_dimids[NC_MAX_VAR_DIMS],

   297 		new_varid, lat_dimid, lon_dimid;

   298 	nc_type xtype, atype;

   299 	char name[NC_MAX_NAME+1], new_name[NC_MAX_NAME+1],

   300 		att_name[NC_MAX_NAME+1], *longname, *new_longname;

   301 	size_t tlen, ts, alen, nlevels, d1, d2;

   302 	float FillValue;

   303 	char FillFlag;

   304 	float *var;

   305 	double *tvar;

   307 #ifdef DEBUG

   308 	printf("Number of total variables (nmtotal_vars) = %d\n", nmtotal_vars);

   309 	printf("Number of combination variables (nmcombo_vars) = %d\n", nmcombo_vars);

   310 #endif /* DEBUG */

   312 	/* Require at least one parameter (the file on which to work) */

   313 	if (argc < 2)

   314 		usage(argv[0]);

   316 	/* Put the list of filenames in an array of strings */

   317 	if (!(ifnames = (char **)malloc((argc - 1) * sizeof(char *)))) {

   318 		perror("ifnames");

   319 		exit(2);

   320 	}

   321 	nifnames = 0;

   322 	for (i = 1; i < argc; i++) {

   323 		if (!(ifnames[nifnames++] = strdup(argv[i]))) {

   324 			perror("ifnames[nifnames++]");

   325 			exit(2);

   326 		}

   327 	}

   329 	/*

   330 	 * March through all input files adding new variables for totals

   331 	 * across levels.

   332 	 */

   333 	for (i = 0; i < nifnames; i++) {

   334 		printf("Processing %s\n", ifnames[i]);

   335 		/* Open file */

   336 		wrap_nc(nc_open(ifnames[i], NC_WRITE, &ncid));

   337 		/*

   338 		 * Inquire about number of dimensions, variables, global

   339 		 * attributes and the ID of the unlimited dimension

   340 		 */

   341 		wrap_nc(nc_inq(ncid, &ngdims, &nvars, &ngatts, &unlimdimid));

   342 		wrap_nc(nc_inq_dim(ncid, unlimdimid, name, &tlen));

   343 		if (strcmp(name, time_name)) {

   344 			fprintf(stderr, "%s is no the unlimited dimension for file %s!\n", time_name, ifnames[i]);

   345 			exit(4);

   346 		}

   347 		/* Retrieve lat and lon dimension IDs */

   348 		wrap_nc(nc_inq_dimid(ncid, lat_name, &lat_dimid));

   349 		wrap_nc(nc_inq_dimid(ncid, lon_name, &lon_dimid));

   350 		/* Put file into define mode */

   351 		wrap_nc(nc_redef(ncid));

   352 		/*

   353 		 * Define all of the new variables first to improve

   354 		 * performance.  This means some calls are made twice.

   355 		 */

   356 		/* First, add the new total variables */

   357 		for (j = 0; j < nmtotal_vars; j++) {

   358 			/* Figure out source variable information */

   359 			wrap_nc(nc_inq_varid(ncid, total_vars[j], &varid));

   360 			wrap_nc(nc_inq_var(ncid, varid, name, &xtype, &ndims, dimids, &natts));

   361 			/* Make sure it has the correct number of dimensions */

   362 			if (ndims != 4) {

   363 				fprintf(stderr, "Variable %s has %d dimensions!\n", name, ndims);

   364 				exit(4);

   365 			}

   366 			/* Make sure time is the first dimension */

   367 			if (dimids[0] != unlimdimid) {

   368 				fprintf(stderr, "First dimension of variable %s is not %s!\n", name, time_name);

   369 				exit(4);

   370 			}

   371 			/* Make sure it is a float type */

   372 			if (xtype != NC_FLOAT) {

   373 				fprintf(stderr, "Only NC_FLOAT type is supported presently!\n");

   374 				exit(4);

   375 			}

   376 			/* Set dimensions for new variable */

   377 			new_ndims = 3;

   378 			new_dimids[0] = unlimdimid;

   379 			new_dimids[1] = dimids[2];

   380 			new_dimids[2] = dimids[3];

   381 			/* Define new variable */

   382 			sprintf(new_name, "%s%s", varname_total_prefix, name);

   383 			printf("\tAdding variable %s\n", new_name);

   384 			wrap_nc(nc_def_var(ncid, new_name, xtype, new_ndims,

   385 				new_dimids, &new_varid));

   386 			/* Copy attributes from source variable to new one */

   387 			for (k = 0; k < natts; k++) {

   388 				wrap_nc(nc_inq_attname(ncid, varid, k, att_name));

   389 				if (!strcmp(att_name, "long_name")) {

   390 					wrap_nc(nc_inq_att(ncid, varid, att_name, &atype, &alen));

   391 					if (!(longname = (char *)malloc((alen+1)*sizeof(char)))) {

   392 						perror("longname");

   393 						exit(2);

   394 					}

   395 					wrap_nc(nc_get_att_text(ncid, varid, att_name, longname));

   396 					longname[alen] = '\0';

   397 					if (!(new_longname = (char *)malloc((strlen(longname_total_prefix)+alen+1)*sizeof(char)))) {

   398 						perror("new_longname");

   399 						exit(2);

   400 					}

   401 					sprintf(new_longname, "%s%s", longname_total_prefix, longname);

   402 					wrap_nc(nc_put_att_text(ncid, new_varid, att_name, strlen(new_longname), new_longname));

   403 					free(new_longname);

   404 					free(longname);

   405 				}

   406 				else

   407 					wrap_nc(nc_copy_att(ncid, varid, att_name, ncid, new_varid));

   408 			}

   409 		}

   410 		/* Second, add the new total variables */

   411 		for (j = 0; j < nmcombo_vars; j++) {

   412 			/* Set dimensions for new variable */

   413 			new_ndims = 3;

   414 			new_dimids[0] = unlimdimid;

   415 			new_dimids[1] = lat_dimid;

   416 			new_dimids[2] = lon_dimid;

   417 			/* Define new variable */

   418 			printf("\tAdding variable %s\n", combo_vars[j]);

   419 			wrap_nc(nc_def_var(ncid, combo_vars[j], NC_FLOAT,

   420 				new_ndims, new_dimids, &new_varid));

   421 			/* Add attributes */

   422 			wrap_nc(nc_put_att_text(ncid, new_varid, long_name_name,

   423 				strlen(combo_long_names[j]),

   424 				combo_long_names[j]));

   425 			wrap_nc(nc_put_att_text(ncid, new_varid, units_name,

   426 				strlen(combo_units[j]), combo_units[j]));

   427 			wrap_nc(nc_put_att_text(ncid, new_varid,

   428 				cell_method_name, strlen(cell_method),

   429 				cell_method));

   430 			wrap_nc(nc_put_att_float(ncid, new_varid,

   431 				FillValue_name, NC_FLOAT, 1,

   432 				&combo_FillValues[j]));

   433 			wrap_nc(nc_put_att_float(ncid, new_varid,

   434 				missing_value_name, NC_FLOAT, 1,

   435 				&combo_FillValues[j]));

   436 		}

   437 		/* Leave define mode */

   438 		wrap_nc(nc_enddef(ncid));

   439 		/* Now actually compute and write out the new total variables */

   440 		for (j = 0; j < nmtotal_vars; j++) {

   441 			/* Figure out source variable information */

   442 			wrap_nc(nc_inq_varid(ncid, total_vars[j], &varid));

   443 			wrap_nc(nc_inq_var(ncid, varid, name, &xtype, &ndims, dimids, &natts));

   444 			/* Check for _FillValue */

   445 			FillFlag = 0;

   446 			FillValue = 0.;

   447 			if (nc_inq_att(ncid, varid, FillValue_name, &atype, &alen) == NC_NOERR) {

   448 				if (atype == NC_FLOAT && alen) {

   449 					wrap_nc(nc_get_att_float(ncid, varid,

   450 						FillValue_name, &FillValue));

   451 					FillFlag = 1;

   452 				}

   453 			}

   454 			/* Set dimensions for new variable */

   455 			new_ndims = 3;

   456 			new_dimids[0] = unlimdimid;

   457 			new_dimids[1] = dimids[2];

   458 			new_dimids[2] = dimids[3];

   460 			sprintf(new_name, "%s%s", varname_total_prefix, name);

   461 			printf("\tComputing and writing %s\n", new_name);

   462 			/* Retrieve the new varid again */

   463 			wrap_nc(nc_inq_varid(ncid, new_name, &new_varid));

   464 			/*

   465 			 * Figure out dimensions and total size

   466 			 * for a single time slice

   467 			 */

   468 			wrap_nc(nc_inq_dimlen(ncid, dimids[1], &nlevels));

   469 			wrap_nc(nc_inq_dimlen(ncid, dimids[2], &d1));

   470 			wrap_nc(nc_inq_dimlen(ncid, dimids[3], &d2));

   472 			/*

   473 			 * Allocate space for reading in time slice and for

   474 			 * the sum

   475 			 */

   476 			if (!(var = (float *)malloc((d1*d2) * sizeof(float)))) {

   477 				perror("var");

   478 				exit(2);

   479 			}

   480 			if (!(tvar = (double *)malloc((d1*d2) * sizeof(double)))) {

   481 				perror("tvar");

   482 				exit(2);

   483 			}

   485 			/* Read timeslices, total, and write out new timeslices */

   486 			for (ts = 0; ts < tlen; ts++) {

   487 				read_total_timeslice(ncid, varid,

   488 					nlevels, d1, d2, FillFlag, FillValue,

   489 					ts, var, tvar);

   490 				for (k = 0; k < (d1*d2); k++)

   491 					var[k] = (float)tvar[k];

   492 				write_timeslice(ncid, new_varid, d1, d2,

   493 					ts, var);

   494 			}

   496 			free(var);

   497 			free(tvar);

   498 		}

   499 		/* Now actually compute and write out the new combo variables */

   500 		for (j = 0; j < nmcombo_vars; j++) {

   501 			/* Set dimensions for new variable */

   502 			new_ndims = 3;

   503 			new_dimids[0] = unlimdimid;

   504 			new_dimids[1] = lat_dimid;

   505 			new_dimids[2] = lon_dimid;

   507 			printf("\tComputing and writing %s\n", combo_vars[j]);

   508 			/* Retrieve the new varid again */

   509 			wrap_nc(nc_inq_varid(ncid, combo_vars[j], &new_varid));

   510 			/*

   511 			 * Figure out dimensions and total size

   512 			 * for a single time slice

   513 			 */

   514 			wrap_nc(nc_inq_dimlen(ncid, new_dimids[1], &d1));

   515 			wrap_nc(nc_inq_dimlen(ncid, new_dimids[2], &d2));

   517 			/*

   518 			 * Allocate space for reading in time slice and for

   519 			 * the sum

   520 			 */

   521 			if (!(var = (float *)malloc((d1*d2) * sizeof(float)))) {

   522 				perror("var");

   523 				exit(2);

   524 			}

   525 			if (!(tvar = (double *)malloc((d1*d2) * sizeof(double)))) {

   526 				perror("tvar");

   527 				exit(2);

   528 			}

   530 			/* Read timeslices, compute new variable, and write */

   531 			for (ts = 0; ts < tlen; ts++) {

   532 				read_compute_timeslice(ncid, d1, d2, ts,

   533 					j, tvar);

   534 				for (k = 0; k < (d1*d2); k++)

   535 					var[k] = (float)tvar[k];

   536 				write_timeslice(ncid, new_varid, d1, d2,

   537 					ts, var);

   538 			}

   540 			free(var);

   541 			free(tvar);

   542 		}

   543 		/* Close file */

   544 		wrap_nc(nc_close(ncid));

   545 	}

   547 	return 0;

   548 }

author	Forrest Hoffman <forrest@climatemodeling.org>
	Wed, 03 Oct 2007 11:23:02 -0400
changeset 3	d3122367777b
permissions	-rw-r--r--