var/convertor/update_msf/update_msf.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <netcdf.h>
   4
   5 /*
   6 To compile (after setting NETCDF environment variable):
   7 cc -o update_msf update_msf.c -I${NETCDF}/include -L${NETCDF}/lib -lnetcdf
   8
   9 To run:
  10 First, make a backup copy of the NetCDF file, since this program will modify the original.
  11 Then, simply run "update_msf <your netcdf file>".
  12 */
  13
  14 int main(int argc, char ** argv) {
  15
  16    int i, j;
  17    int ip;
  18    int istatus, ncid;
  19    int varid, msft_id, msfu_id, msfv_id, msftx_id, msfty_id, msfux_id, msfuy_id, msfvx_id, msfvy_id, msfvy_inv_id;
  20    int we_id, we_stag_id, sn_id, sn_stag_id;
  21    size_t we_len, we_stag_len, sn_len, sn_stag_len;
  22    int dimids[3];
  23    float * msft, * msfu, * msfv;
  24    char varname[1024];
  25
  26    if (argc != 2) {
  27       fprintf(stderr,"\nUsage: update_msf <file>\n\n");
  28       return 1;
  29    }
  30
  31    /* Try to open the file */
  32    istatus = nc_open(argv[1], NC_WRITE, &ncid);
  33    if (istatus == NC_NOERR) {
  34
  35       /* Get IDs of dimensions */
  36       istatus = nc_inq_dimid(ncid, "west_east", &we_id);
  37       istatus |= nc_inq_dimid(ncid, "south_north", &sn_id);
  38       istatus |= nc_inq_dimid(ncid, "west_east_stag", &we_stag_id);
  39       istatus |= nc_inq_dimid(ncid, "south_north_stag", &sn_stag_id);
  40       if (istatus != NC_NOERR) {
  41          fprintf(stderr,"Error: Could not get dimensions for map scale factor fields\n");
  42          return 1;
  43       }
  44
  45       /* Get lengths of the we and sn dimensions */
  46       istatus = nc_inq_dimlen(ncid, we_id, &we_len);
  47       istatus |= nc_inq_dimlen(ncid, sn_id, &sn_len);
  48       istatus |= nc_inq_dimlen(ncid, we_stag_id, &we_stag_len);
  49       istatus |= nc_inq_dimlen(ncid, sn_stag_id, &sn_stag_len);
  50       if (istatus != NC_NOERR) {
  51          fprintf(stderr,"Error: Could not get length of dimensions for map scale factor fields\n");
  52          return 1;
  53       }
  54
  55       msft = (float *)malloc(sizeof(float)*we_len*sn_len);
  56       msfu = (float *)malloc(sizeof(float)*we_stag_len*sn_len);
  57       msfv = (float *)malloc(sizeof(float)*we_len*sn_stag_len);
  58
  59       /* Read in MSFT */
  60       sprintf(varname,"MAPFAC_M");
  61       istatus = nc_inq_varid (ncid, varname, &varid);
  62       if (istatus != NC_NOERR) {
  63          fprintf(stderr,"Error: Could not find variable %s in %s!\n",varname,argv[1]);
  64          return 1;
  65       }
  66       istatus = nc_get_var_float(ncid, varid, msft);
  67       if (istatus != NC_NOERR) {
  68          fprintf(stderr,"Error: Could not read MSFT for %s in %s!\n",varname,argv[1]);
  69          free(msft);
  70          free(msfu);
  71          free(msfv);
  72          return 1;
  73       }
  74
  75       /* Read in MSFU */
  76       sprintf(varname,"MAPFAC_U");
  77       istatus = nc_inq_varid (ncid, varname, &varid);
  78       if (istatus != NC_NOERR) {
  79          fprintf(stderr,"Error: Could not find variable %s in %s!\n",varname,argv[1]);
  80          return 1;
  81       }
  82       istatus = nc_get_var_float(ncid, varid, msfu);
  83       if (istatus != NC_NOERR) {
  84          fprintf(stderr,"Error: Could not read MSFT for %s in %s!\n",varname,argv[1]);
  85          free(msft);
  86          free(msfu);
  87          free(msfv);
  88          return 1;
  89       }
  90
  91       /* Read in MSFV */
  92       sprintf(varname,"MAPFAC_V");
  93       istatus = nc_inq_varid (ncid, varname, &varid);
  94       if (istatus != NC_NOERR) {
  95          fprintf(stderr,"Error: Could not find variable %s in %s!\n",varname,argv[1]);
  96          return 1;
  97       }
  98       istatus = nc_get_var_float(ncid, varid, msfv);
  99       if (istatus != NC_NOERR) {
 100          fprintf(stderr,"Error: Could not read MSFT for %s in %s!\n",varname,argv[1]);
 101          free(msft);
 102          free(msfu);
 103          free(msfv);
 104          return 1;
 105       }
 106
 107       /* Write new fields back out to file */
 108       istatus = nc_redef(ncid);
 109
 110       dimids[0] = NC_UNLIMITED;
 111
 112       dimids[1] = sn_id;
 113       dimids[2] = we_id;
 114       istatus = nc_def_var(ncid, "MAPFAC_MX", NC_FLOAT, 3, dimids, &msftx_id);
 115       istatus = nc_def_var(ncid, "MAPFAC_MY", NC_FLOAT, 3, dimids, &msfty_id);
 116
 117       dimids[1] = sn_id;
 118       dimids[2] = we_stag_id;
 119       istatus = nc_def_var(ncid, "MAPFAC_UX", NC_FLOAT, 3, dimids, &msfux_id);
 120       istatus = nc_def_var(ncid, "MAPFAC_UY", NC_FLOAT, 3, dimids, &msfuy_id);
 121
 122       dimids[1] = sn_stag_id;
 123       dimids[2] = we_id;
 124       istatus = nc_def_var(ncid, "MAPFAC_VX", NC_FLOAT, 3, dimids, &msfvx_id);
 125       istatus = nc_def_var(ncid, "MAPFAC_VY", NC_FLOAT, 3, dimids, &msfvy_id);
 126       istatus = nc_def_var(ncid, "MF_VX_INV", NC_FLOAT, 3, dimids, &msfvy_inv_id);
 127
 128       ip = 104;
 129       istatus = nc_put_att_int(ncid, msftx_id, "FieldType", NC_INT, 1, &ip);
 130       istatus = nc_put_att_int(ncid, msfty_id, "FieldType", NC_INT, 1, &ip);
 131       istatus = nc_put_att_int(ncid, msfux_id, "FieldType", NC_INT, 1, &ip);
 132       istatus = nc_put_att_int(ncid, msfuy_id, "FieldType", NC_INT, 1, &ip);
 133       istatus = nc_put_att_int(ncid, msfvx_id, "FieldType", NC_INT, 1, &ip);
 134       istatus = nc_put_att_int(ncid, msfvy_id, "FieldType", NC_INT, 1, &ip);
 135       istatus = nc_put_att_int(ncid, msfvy_inv_id, "FieldType", NC_INT, 1, &ip);
 136
 137       istatus = nc_put_att_text(ncid, msftx_id, "MemoryOrder", 3, "XY ");
 138       istatus = nc_put_att_text(ncid, msfty_id, "MemoryOrder", 3, "XY ");
 139       istatus = nc_put_att_text(ncid, msfux_id, "MemoryOrder", 3, "XY ");
 140       istatus = nc_put_att_text(ncid, msfuy_id, "MemoryOrder", 3, "XY ");
 141       istatus = nc_put_att_text(ncid, msfvx_id, "MemoryOrder", 3, "XY ");
 142       istatus = nc_put_att_text(ncid, msfvy_id, "MemoryOrder", 3, "XY ");
 143       istatus = nc_put_att_text(ncid, msfvy_inv_id, "MemoryOrder", 3, "XY ");
 144
 145       istatus = nc_put_att_text(ncid, msftx_id, "description", 29, "Map scale factor on mass grid");
 146       istatus = nc_put_att_text(ncid, msfty_id, "description", 29, "Map scale factor on mass grid");
 147       istatus = nc_put_att_text(ncid, msfux_id, "description", 26, "Map scale factor on u-grid");
 148       istatus = nc_put_att_text(ncid, msfuy_id, "description", 26, "Map scale factor on u-grid");
 149       istatus = nc_put_att_text(ncid, msfvx_id, "description", 26, "Map scale factor on v-grid");
 150       istatus = nc_put_att_text(ncid, msfvy_id, "description", 26, "Map scale factor on v-grid");
 151       istatus = nc_put_att_text(ncid, msfvy_inv_id, "description", 26, "Map scale factor on v-grid");
 152
 153       istatus = nc_put_att_text(ncid, msftx_id, "units", 0, "");
 154       istatus = nc_put_att_text(ncid, msfty_id, "units", 0, "");
 155       istatus = nc_put_att_text(ncid, msfux_id, "units", 0, "");
 156       istatus = nc_put_att_text(ncid, msfuy_id, "units", 0, "");
 157       istatus = nc_put_att_text(ncid, msfvx_id, "units", 0, "");
 158       istatus = nc_put_att_text(ncid, msfvy_id, "units", 0, "");
 159       istatus = nc_put_att_text(ncid, msfvy_inv_id, "units", 0, "");
 160
 161       istatus = nc_put_att_text(ncid, msftx_id, "stagger", 0, "");
 162       istatus = nc_put_att_text(ncid, msfty_id, "stagger", 0, "");
 163       istatus = nc_put_att_text(ncid, msfux_id, "stagger", 1, "X");
 164       istatus = nc_put_att_text(ncid, msfuy_id, "stagger", 1, "X");
 165       istatus = nc_put_att_text(ncid, msfvx_id, "stagger", 1, "Y");
 166       istatus = nc_put_att_text(ncid, msfvy_id, "stagger", 1, "Y");
 167       istatus = nc_put_att_text(ncid, msfvy_inv_id, "stagger", 1, "Y");
 168
 169       istatus = nc_put_att_text(ncid, msftx_id, "coordianates", 10, "XLONG XLAT");
 170       istatus = nc_put_att_text(ncid, msfty_id, "coordianates", 10, "XLONG XLAT");
 171       istatus = nc_put_att_text(ncid, msfux_id, "coordianates", 14, "XLONG_U XLAT_U");
 172       istatus = nc_put_att_text(ncid, msfuy_id, "coordianates", 14, "XLONG_U XLAT_U");
 173       istatus = nc_put_att_text(ncid, msfvx_id, "coordianates", 14, "XLONG_V XLAT_V");
 174       istatus = nc_put_att_text(ncid, msfvy_id, "coordianates", 14, "XLONG_V XLAT_V");
 175       istatus = nc_put_att_text(ncid, msfvy_inv_id, "coordianates", 14, "XLONG_V XLAT_V");
 176
 177       istatus = nc_enddef(ncid);
 178
 179       istatus = nc_put_var_float(ncid, msftx_id, msft);
 180       istatus |= nc_put_var_float(ncid, msfty_id, msft);
 181
 182       istatus |= nc_put_var_float(ncid, msfux_id, msfu);
 183       istatus |= nc_put_var_float(ncid, msfuy_id, msfu);
 184
 185       istatus |= nc_put_var_float(ncid, msfvx_id, msfv);
 186       istatus |= nc_put_var_float(ncid, msfvy_id, msfv);
 187
 188       for(j=0; j<sn_stag_len; j++) {
 189          for(i=0; i<we_len; i++) {
 190             if (msfv[i+we_len*j] != 0.)
 191                msfv[i+we_len*j] = 1.0/msfv[i+we_len*j];
 192             else
 193                msfv[i+we_len*j] = 0.0;
 194          }
 195       }
 196       istatus |= nc_put_var_float(ncid, msfvy_inv_id, msfv);
 197
 198       if (istatus != NC_NOERR) {
 199          fprintf(stderr,"Error: Could not write array to %s!\n",argv[1]);
 200          free(msft);
 201          free(msfu);
 202          free(msfv);
 203          return 1;
 204       }
 205
 206       /* Close file */
 207       free(msft);
 208       free(msfu);
 209       free(msfv);
 210       istatus = nc_close(ncid);
 211    }
 212    else {
 213       fprintf(stderr,"Error: Could not open %s as a NetCDF file!\n",argv[1]);
 214       return 1;
 215    }
 216
 217    return 0;
 218 }