mosaic/read_mosaic.c

   1 /***********************************************************************
   2  *                   GNU Lesser General Public License
   3  *
   4  * This file is part of the GFDL Flexible Modeling System (FMS).
   5  *
   6  * FMS is free software: you can redistribute it and/or modify it under
   7  * the terms of the GNU Lesser General Public License as published by
   8  * the Free Software Foundation, either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * FMS is distributed in the hope that it will be useful, but WITHOUT
  12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14  * for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public
  17  * License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
  18  **********************************************************************/
  19 #include <stdlib.h>
  20 #include <stdio.h>
  21 #include <math.h>
  22 #include <string.h>
  23 #include "read_mosaic.h"
  24 #include "constant.h"
  25 #include "mosaic_util.h"
  26 #include <netcdf.h>
  27
  28 /** \file
  29  *  \ingroup mosaic
  30  *  \brief Support for reading mosaic netcdf grid files.
  31  */
  32
  33 /*********************************************************************
  34     void netcdf_error( int status )
  35     status is the returning value of netcdf call. this routine will
  36     handle the error when status is not NC_NOERR.
  37 ********************************************************************/
  38 void handle_netcdf_error(const char *msg, int status )
  39 {
  40   char errmsg[512];
  41
  42   sprintf( errmsg, "%s: %s", msg, (char *)nc_strerror(status) );
  43   error_handler(errmsg);
  44
  45 } /* handle_netcdf_error */
  46
  47 /***************************************************************************
  48   void get_file_dir(const char *file, char *dir)
  49   get the directory where file is located. The dir will be the complate path
  50   before the last "/". If no "/" exist in file, the path will be current ".".
  51 ***************************************************************************/
  52 void get_file_dir(const char *file, char *dir)
  53 {
  54   int len;
  55   const char *strptr = NULL;
  56
  57   /* get the diretory */
  58
  59   strptr = strrchr(file, '/');
  60   if(strptr) {
  61     len = strptr - file;
  62     strncpy(dir, file, len);
  63   }
  64   else {
  65     len = 1;
  66     strcpy(dir, ".");
  67   }
  68   dir[len] = 0;
  69
  70 } /* get_file_dir */
  71
  72
  73 int field_exist(const char* file, const char *name)
  74 {
  75   int ncid, varid, status;
  76   char msg[512];
  77   int existed=0;
  78
  79 #ifdef use_netCDF
  80
  81   status = nc_open(file, NC_NOWRITE, &ncid);
  82   if(status != NC_NOERR) {
  83     sprintf(msg, "field_exist: in opening file %s", file);
  84     handle_netcdf_error(msg, status);
  85   }
  86
  87   status = nc_inq_varid(ncid, name, &varid);
  88   if(status == NC_NOERR){
  89     existed = 1;
  90   }
  91
  92   status = nc_close(ncid);
  93   if(status != NC_NOERR) {
  94     sprintf(msg, "field_exist: in closing file %s.", file);
  95     handle_netcdf_error(msg, status);
  96   }
  97
  98 #else  /* ndef use_netCDF */
  99   error_handler("read_mosaic: Add flag -Duse_netCDF when compiling");
 100 #endif  /* use_netcdf */
 101
 102   return existed;
 103
 104 } /* field_exist */
 105
 106 int get_dimlen(const char* file, const char *name)
 107 {
 108   int ncid, dimid, status, len;
 109   size_t size;
 110   char msg[512];
 111
 112   len = 0;
 113 #ifdef use_netCDF
 114   status = nc_open(file, NC_NOWRITE, &ncid);
 115   if(status != NC_NOERR) {
 116     sprintf(msg, "in opening file %s", file);
 117     handle_netcdf_error(msg, status);
 118   }
 119
 120   status = nc_inq_dimid(ncid, name, &dimid);
 121   if(status != NC_NOERR) {
 122     sprintf(msg, "in getting dimid of %s from file %s.", name, file);
 123     handle_netcdf_error(msg, status);
 124   }
 125
 126   status = nc_inq_dimlen(ncid, dimid, &size);
 127   if(status != NC_NOERR) {
 128     sprintf(msg, "in getting dimension size of %s from file %s.", name, file);
 129     handle_netcdf_error(msg, status);
 130   }
 131   status = nc_close(ncid);
 132   if(status != NC_NOERR) {
 133     sprintf(msg, "in closing file %s.", file);
 134     handle_netcdf_error(msg, status);
 135   }
 136
 137   len = size;
 138   if(status != NC_NOERR) {
 139     sprintf(msg, "in closing file %s", file);
 140     handle_netcdf_error(msg, status);
 141   }
 142 #else
 143   error_handler("read_mosaic: Add flag -Duse_netCDF when compiling");
 144 #endif
 145
 146   return len;
 147
 148 } /* get_dimlen */
 149
 150 /*******************************************************************************
 151    void get_string_data(const char *file, const char *name, char *data)
 152    get string data of field with "name" from "file".
 153 ******************************************************************************/
 154 void get_string_data(const char *file, const char *name, char *data)
 155 {
 156   int ncid, varid, status;
 157   char msg[512];
 158
 159 #ifdef use_netCDF
 160   status = nc_open(file, NC_NOWRITE, &ncid);
 161   if(status != NC_NOERR) {
 162     sprintf(msg, "in opening file %s", file);
 163     handle_netcdf_error(msg, status);
 164   }
 165   status = nc_inq_varid(ncid, name, &varid);
 166   if(status != NC_NOERR) {
 167     sprintf(msg, "in getting varid of %s from file %s.", name, file);
 168     handle_netcdf_error(msg, status);
 169   }
 170   status = nc_get_var_text(ncid, varid, data);
 171   if(status != NC_NOERR) {
 172     sprintf(msg, "in getting data of %s from file %s.", name, file);
 173     handle_netcdf_error(msg, status);
 174   }
 175   status = nc_close(ncid);
 176   if(status != NC_NOERR) {
 177     sprintf(msg, "in closing file %s.", file);
 178     handle_netcdf_error(msg, status);
 179   }
 180 #else
 181   error_handler("read_mosaic: Add flag -Duse_netCDF when compiling");
 182 #endif
 183
 184 } /* get_string_data */
 185
 186 /*******************************************************************************
 187    void get_string_data_level(const char *file, const char *name, const size_t *start, const size_t *nread, char *data)
 188    get string data of field with "name" from "file".
 189 ******************************************************************************/
 190 void get_string_data_level(const char *file, const char *name, char *data, const unsigned int *level)
 191 {
 192   int ncid, varid, status, i;
 193   size_t start[4], nread[4];
 194   char msg[512];
 195
 196 #ifdef use_netCDF
 197   status = nc_open(file, NC_NOWRITE, &ncid);
 198   if(status != NC_NOERR) {
 199     sprintf(msg, "in opening file %s", file);
 200     handle_netcdf_error(msg, status);
 201   }
 202   status = nc_inq_varid(ncid, name, &varid);
 203   if(status != NC_NOERR) {
 204     sprintf(msg, "in getting varid of %s from file %s.", name, file);
 205     handle_netcdf_error(msg, status);
 206   }
 207   for(i=0; i<4; i++) {
 208     start[i] = 0; nread[i] = 1;
 209   }
 210   start[0] = *level; nread[1] = STRING;
 211   status = nc_get_vara_text(ncid, varid, start, nread, data);
 212   if(status != NC_NOERR) {
 213     sprintf(msg, "in getting data of %s from file %s.", name, file);
 214     handle_netcdf_error(msg, status);
 215   }
 216   status = nc_close(ncid);
 217   if(status != NC_NOERR) {
 218     sprintf(msg, "in closing file %s.", file);
 219     handle_netcdf_error(msg, status);
 220   }
 221 #else
 222   error_handler("read_mosaic: Add flag -Duse_netCDF when compiling");
 223 #endif
 224
 225 } /* get_string_data_level */
 226
 227
 228 /*******************************************************************************
 229    void get_var_data(const char *file, const char *name, double *data)
 230    get var data of field with "name" from "file".
 231 ******************************************************************************/
 232 void get_var_data(const char *file, const char *name, void *data)
 233 {
 234
 235   int ncid, varid, status;
 236   nc_type vartype;
 237   char msg[512];
 238
 239 #ifdef use_netCDF
 240   status = nc_open(file, NC_NOWRITE, &ncid);
 241   if(status != NC_NOERR) {
 242     sprintf(msg, "in opening file %s", file);
 243     handle_netcdf_error(msg, status);
 244   }
 245   status = nc_inq_varid(ncid, name, &varid);
 246   if(status != NC_NOERR) {
 247     sprintf(msg, "in getting varid of %s from file %s.", name, file);
 248     handle_netcdf_error(msg, status);
 249   }
 250
 251   status = nc_inq_vartype(ncid, varid, &vartype);
 252   if(status != NC_NOERR) {
 253     sprintf(msg, "get_var_data: in getting vartype of of %s in file %s ", name, file);
 254     handle_netcdf_error(msg, status);
 255   }
 256
 257   switch (vartype) {
 258     case NC_DOUBLE:case NC_FLOAT:
 259       status = nc_get_var_double(ncid, varid, (double *)data);
 260       break;
 261     case NC_INT:
 262       status = nc_get_var_int(ncid, varid, (int *)data);
 263       break;
 264     default:
 265       sprintf(msg, "get_var_data: field %s in file %s has an invalid type, "
 266               "the type should be NC_DOUBLE, NC_FLOAT or NC_INT", name, file);
 267       error_handler(msg);
 268   }
 269   if(status != NC_NOERR) {
 270     sprintf(msg, "in getting data of %s from file %s.", name, file);
 271     handle_netcdf_error(msg, status);
 272   }
 273   status = nc_close(ncid);
 274   if(status != NC_NOERR) {
 275     sprintf(msg, "in closing file %s.", file);
 276     handle_netcdf_error(msg, status);
 277   }
 278 #else
 279   error_handler("read_mosaic: Add flag -Duse_netCDF when compiling");
 280 #endif
 281
 282 } /* get_var_data */
 283
 284 /*******************************************************************************
 285    void get_var_data(const char *file, const char *name, double *data)
 286    get var data of field with "name" from "file".
 287 ******************************************************************************/
 288 void get_var_data_region(const char *file, const char *name, const size_t *start, const size_t *nread, void *data)
 289 {
 290
 291   int ncid, varid, status;
 292   nc_type vartype;
 293   char msg[512];
 294
 295 #ifdef use_netCDF
 296   status = nc_open(file, NC_NOWRITE, &ncid);
 297   if(status != NC_NOERR) {
 298     sprintf(msg, "get_var_data_region: in opening file %s", file);
 299     handle_netcdf_error(msg, status);
 300   }
 301   status = nc_inq_varid(ncid, name, &varid);
 302   if(status != NC_NOERR) {
 303     sprintf(msg, "in getting varid of %s from file %s.", name, file);
 304     handle_netcdf_error(msg, status);
 305   }
 306
 307   status = nc_inq_vartype(ncid, varid, &vartype);
 308   if(status != NC_NOERR) {
 309     sprintf(msg, "get_var_data_region: in getting vartype of of %s in file %s ", name, file);
 310     handle_netcdf_error(msg, status);
 311   }
 312
 313   switch (vartype) {
 314     case NC_DOUBLE:case NC_FLOAT:
 315       status = nc_get_vara_double(ncid, varid, start, nread, (double *)data);
 316       break;
 317     case NC_INT:
 318       status = nc_get_vara_int(ncid, varid, start, nread, (int *)data);
 319       break;
 320     default:
 321       sprintf(msg, "get_var_data_region: field %s in file %s has an invalid type, "
 322               "the type should be NC_DOUBLE, NC_FLOAT or NC_INT", name, file);
 323       error_handler(msg);
 324   }
 325
 326   if(status != NC_NOERR) {
 327     sprintf(msg, "get_var_data_region: in getting data of %s from file %s.", name, file);
 328     handle_netcdf_error(msg, status);
 329   }
 330   status = nc_close(ncid);
 331   if(status != NC_NOERR) {
 332     sprintf(msg, "get_var_data_region: in closing file %s.", file);
 333     handle_netcdf_error(msg, status);
 334   }
 335 #else
 336   error_handler("read_mosaic: Add flag -Duse_netCDF when compiling");
 337 #endif
 338
 339 } /* get_var_data_region */
 340
 341 /******************************************************************************
 342    void get_var_text_att(const char *file, const char *name, const char *attname, char *att)
 343    get text attribute of field 'name' from 'file
 344 ******************************************************************************/
 345 void get_var_text_att(const char *file, const char *name, const char *attname, char *att)
 346 {
 347   int ncid, varid, status;
 348   char msg[512];
 349
 350 #ifdef use_netCDF
 351   status = nc_open(file, NC_NOWRITE, &ncid);
 352   if(status != NC_NOERR) {
 353     sprintf(msg, "in opening file %s", file);
 354     handle_netcdf_error(msg, status);
 355   }
 356   status = nc_inq_varid(ncid, name, &varid);
 357   if(status != NC_NOERR) {
 358     sprintf(msg, "in getting varid of %s from file %s.", name, file);
 359     handle_netcdf_error(msg, status);
 360   }
 361   status = nc_get_att_text(ncid, varid, attname, att);
 362   if(status != NC_NOERR) {
 363     sprintf(msg, "in getting attribute %s of %s from file %s.", attname, name, file);
 364     handle_netcdf_error(msg, status);
 365   }
 366   status = nc_close(ncid);
 367   if(status != NC_NOERR) {
 368     sprintf(msg, "in closing file %s.", file);
 369     handle_netcdf_error(msg, status);
 370   }
 371 #else
 372   error_handler("read_mosaic: Add flag -Duse_netCDF when compiling");
 373 #endif
 374
 375 } /* get_var_text_att */
 376
 377 /***********************************************************************
 378   return number of overlapping cells.
 379 ***********************************************************************/
 380 int read_mosaic_xgrid_size_( const char *xgrid_file )
 381 {
 382   return read_mosaic_xgrid_size(xgrid_file);
 383 }
 384
 385 int read_mosaic_xgrid_size( const char *xgrid_file )
 386 {
 387   int ncells;
 388
 389   ncells = get_dimlen(xgrid_file, "ncells");
 390   return ncells;
 391 }
 392
 393   double get_global_area(void)
 394   {
 395     double garea;
 396     garea = 4*M_PI*RADIUS*RADIUS;
 397
 398     return garea;
 399   }
 400
 401   double get_global_area_(void)
 402   {
 403     double garea;
 404     garea = 4*M_PI*RADIUS*RADIUS;
 405
 406     return garea;
 407   }
 408
 409
 410     /****************************************************************************/
 411   void read_mosaic_xgrid_order1_(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, double *area )
 412     {
 413       read_mosaic_xgrid_order1(xgrid_file, i1, j1, i2, j2, area);
 414
 415     }
 416
 417   void read_mosaic_xgrid_order1(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, double *area )
 418     {
 419       int    ncells, n;
 420       int    *tile1_cell, *tile2_cell;
 421       double garea;
 422
 423       ncells = get_dimlen(xgrid_file, "ncells");
 424
 425       tile1_cell       = (int *)malloc(ncells*2*sizeof(int));
 426       tile2_cell       = (int *)malloc(ncells*2*sizeof(int));
 427       get_var_data(xgrid_file, "tile1_cell", tile1_cell);
 428       get_var_data(xgrid_file, "tile2_cell", tile2_cell);
 429
 430       get_var_data(xgrid_file, "xgrid_area", area);
 431
 432       garea = 4*M_PI*RADIUS*RADIUS;
 433
 434       for(n=0; n<ncells; n++) {
 435         i1[n] = tile1_cell[n*2] - 1;
 436         j1[n] = tile1_cell[n*2+1] - 1;
 437         i2[n] = tile2_cell[n*2] - 1;
 438         j2[n] = tile2_cell[n*2+1] - 1;
 439         area[n] /= garea; /* rescale the exchange grid area to unit earth area */
 440       }
 441
 442       free(tile1_cell);
 443       free(tile2_cell);
 444
 445     } /* read_mosaic_xgrid_order1 */
 446
 447
 448   void read_mosaic_xgrid_order1_region_(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, double *area, int *isc, int *iec )
 449     {
 450       read_mosaic_xgrid_order1_region(xgrid_file, i1, j1, i2, j2, area, isc, iec);
 451
 452     }
 453
 454   void read_mosaic_xgrid_order1_region(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, double *area, int *isc, int *iec )
 455     {
 456       int    ncells, n, i;
 457       int    *tile1_cell, *tile2_cell;
 458       size_t start[4], nread[4];
 459       double garea;
 460
 461       ncells = *iec-*isc+1;
 462
 463       tile1_cell       = (int *)malloc(ncells*2*sizeof(int));
 464       tile2_cell       = (int *)malloc(ncells*2*sizeof(int));
 465       for(i=0; i<4; i++) {
 466         start[i] = 0; nread[i] = 1;
 467       }
 468
 469       start[0] = *isc;
 470       nread[0] = ncells;
 471       nread[1] = 2;
 472
 473       get_var_data_region(xgrid_file, "tile1_cell", start, nread, tile1_cell);
 474       get_var_data_region(xgrid_file, "tile2_cell", start, nread, tile2_cell);
 475
 476       nread[1] = 1;
 477
 478       get_var_data_region(xgrid_file, "xgrid_area", start, nread, area);
 479
 480       garea = 4*M_PI*RADIUS*RADIUS;
 481
 482       for(n=0; n<ncells; n++) {
 483         i1[n] = tile1_cell[n*2] - 1;
 484         j1[n] = tile1_cell[n*2+1] - 1;
 485         i2[n] = tile2_cell[n*2] - 1;
 486         j2[n] = tile2_cell[n*2+1] - 1;
 487         area[n] /= garea; /* rescale the exchange grid area to unit earth area */
 488       }
 489
 490       free(tile1_cell);
 491       free(tile2_cell);
 492
 493     } /* read_mosaic_xgrid_order1 */
 494
 495     /* NOTE: di, dj is for tile1, */
 496     /****************************************************************************/
 497   void read_mosaic_xgrid_order2_(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, double *area, double *di, double *dj )
 498   {
 499       read_mosaic_xgrid_order2(xgrid_file, i1, j1, i2, j2, area, di, dj);
 500
 501     }
 502   void read_mosaic_xgrid_order2(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, double *area, double *di, double *dj )
 503     {
 504       int    ncells, n;
 505       int    *tile1_cell, *tile2_cell;
 506       double *tile1_distance;
 507       double garea;
 508       ncells = get_dimlen(xgrid_file, "ncells");
 509
 510       tile1_cell       = (int    *)malloc(ncells*2*sizeof(int   ));
 511       tile2_cell       = (int    *)malloc(ncells*2*sizeof(int   ));
 512       tile1_distance   = (double *)malloc(ncells*2*sizeof(double));
 513       get_var_data(xgrid_file, "tile1_cell", tile1_cell);
 514       get_var_data(xgrid_file, "tile2_cell", tile2_cell);
 515       get_var_data(xgrid_file, "xgrid_area", area);
 516       get_var_data(xgrid_file, "tile1_distance", tile1_distance);
 517
 518       garea = 4*M_PI*RADIUS*RADIUS;
 519
 520       for(n=0; n<ncells; n++) {
 521         i1[n] = tile1_cell[n*2] - 1;
 522         j1[n] = tile1_cell[n*2+1] - 1;
 523         i2[n] = tile2_cell[n*2] - 1;
 524         j2[n] = tile2_cell[n*2+1] - 1;
 525         di[n] = tile1_distance[n*2];
 526         dj[n] = tile1_distance[n*2+1];
 527         area[n] /= garea; /* rescale the exchange grid area to unit earth area */
 528       }
 529
 530       free(tile1_cell);
 531       free(tile2_cell);
 532       free(tile1_distance);
 533
 534     } /* read_mosaic_xgrid_order2 */
 535
 536     /******************************************************************************
 537   int read_mosaic_ntiles(const char *mosaic_file)
 538   return number tiles in mosaic_file
 539     ******************************************************************************/
 540     int read_mosaic_ntiles_(const char *mosaic_file)
 541     {
 542       return read_mosaic_ntiles(mosaic_file);
 543     }
 544     int read_mosaic_ntiles(const char *mosaic_file)
 545     {
 546
 547       int ntiles;
 548
 549       ntiles = get_dimlen(mosaic_file, "ntiles");
 550
 551       return ntiles;
 552
 553     } /* read_mosaic_ntiles */
 554
 555     /******************************************************************************
 556   int read_mosaic_ncontacts(const char *mosaic_file)
 557   return number of contacts in mosaic_file
 558     ******************************************************************************/
 559     int read_mosaic_ncontacts_(const char *mosaic_file)
 560     {
 561       return read_mosaic_ncontacts(mosaic_file);
 562     }
 563     int read_mosaic_ncontacts(const char *mosaic_file)
 564     {
 565
 566       int ncontacts;
 567
 568       if(field_exist(mosaic_file, "contacts") )
 569         ncontacts = get_dimlen(mosaic_file, "ncontact");
 570       else
 571         ncontacts = 0;
 572
 573       return ncontacts;
 574
 575     } /* read_mosaic_ncontacts */
 576
 577
 578     /*****************************************************************************
 579   void read_mosaic_grid_sizes(const char *mosaic_file, int *nx, int *ny)
 580   read mosaic grid size of each tile, currently we are assuming the refinement is 2.
 581   We assume the grid files are located at the same directory as mosaic_file.
 582     *****************************************************************************/
 583     void read_mosaic_grid_sizes_(const char *mosaic_file, int *nx, int *ny)
 584     {
 585       read_mosaic_grid_sizes(mosaic_file, nx, ny);
 586     }
 587     void read_mosaic_grid_sizes(const char *mosaic_file, int *nx, int *ny)
 588     {
 589       unsigned int ntiles, n;
 590       char gridfile[STRING], tilefile[2*STRING];
 591       char dir[STRING];
 592       const int x_refine = 2, y_refine = 2;
 593
 594       get_file_dir(mosaic_file, dir);
 595       ntiles = get_dimlen(mosaic_file, "ntiles");
 596       for(n = 0; n < ntiles; n++) {
 597         get_string_data_level(mosaic_file, "gridfiles", gridfile, &n);
 598         sprintf(tilefile, "%s/%s", dir, gridfile);
 599         nx[n] = get_dimlen(tilefile, "nx");
 600         ny[n] = get_dimlen(tilefile, "ny");
 601         if(nx[n]%x_refine != 0) error_handler("Error from read_mosaic_grid_sizes: nx is not divided by x_refine");
 602         if(ny[n]%y_refine != 0) error_handler("Error from read_mosaic_grid_sizes: ny is not divided by y_refine");
 603         nx[n] /= x_refine;
 604         ny[n] /= y_refine;
 605       }
 606
 607     } /* read_mosaic_grid_sizes */
 608
 609
 610     /******************************************************************************
 611   void read_mosaic_contact(const char *mosaic_file)
 612   read mosaic contact information
 613     ******************************************************************************/
 614     void read_mosaic_contact_(const char *mosaic_file, int *tile1, int *tile2, int *istart1, int *iend1,
 615                               int *jstart1, int *jend1, int *istart2, int *iend2, int *jstart2, int *jend2)
 616     {
 617       read_mosaic_contact(mosaic_file, tile1, tile2, istart1, iend1, jstart1, jend1, istart2, iend2, jstart2, jend2);
 618     }
 619
 620     /* transfer the index from supergrid to model grid
 621        return 0 if istart = iend
 622        otherwise return 1.
 623     */
 624
 625     int transfer_to_model_index(int istart_in, int iend_in, int *istart_out, int *iend_out, int refine_ratio)
 626     {
 627
 628       int type;
 629
 630       if( istart_in == iend_in ) {
 631         type = 0;
 632         istart_out[0] = (istart_in+1)/refine_ratio-1;
 633         iend_out[0]  = istart_out[0];
 634       }
 635       else {
 636         type = 1;
 637         if( iend_in > istart_in ) {
 638           istart_out[0] = istart_in - 1;
 639           iend_out[0]   = iend_in - refine_ratio;
 640         }
 641         else {
 642           istart_out[0] = istart_in - refine_ratio;
 643           iend_out[0]   = iend_in - 1;
 644         }
 645
 646         if( istart_out[0]%refine_ratio || iend_out[0]%refine_ratio)
 647           error_handler("Error from read_mosaic: mismatch between refine_ratio and istart_in/iend_in");
 648         istart_out[0] /= refine_ratio;
 649         iend_out[0]   /= refine_ratio;
 650       }
 651
 652       return type;
 653
 654     }
 655
 656
 657     void read_mosaic_contact(const char *mosaic_file, int *tile1, int *tile2, int *istart1, int *iend1,
 658                              int *jstart1, int *jend1, int *istart2, int *iend2, int *jstart2, int *jend2)
 659     {
 660       char contacts[STRING];
 661       char **gridtiles;
 662 #define MAXVAR 40
 663       char pstring[MAXVAR][STRING];
 664       unsigned int nstr, ntiles, ncontacts, n, m, l, found;
 665       const int x_refine = 2, y_refine = 2;
 666       int i1_type, j1_type, i2_type, j2_type;
 667
 668       ntiles = get_dimlen(mosaic_file, "ntiles");
 669       gridtiles = (char **)malloc(ntiles*sizeof(char *));
 670       for(n=0; n<ntiles; n++) {
 671         gridtiles[n] = (char *)malloc(STRING*sizeof(char));
 672         get_string_data_level(mosaic_file, "gridtiles", gridtiles[n], &n);
 673       }
 674
 675       ncontacts = get_dimlen(mosaic_file, "ncontact");
 676       for(n = 0; n < ncontacts; n++) {
 677         get_string_data_level(mosaic_file, "contacts", contacts, &n);
 678         /* parse the string contacts to get tile number */
 679         tokenize( contacts, ":", STRING, MAXVAR, (char *)pstring, &nstr);
 680         if(nstr != 4) error_handler("Error from read_mosaic: number of elements "
 681                                     "in contact seperated by :/:: should be 4");
 682         found = 0;
 683         for(m=0; m<ntiles; m++) {
 684           if(strcmp(gridtiles[m], pstring[1]) == 0) { /*found the tile name */
 685             found = 1;
 686             tile1[n] = m+1;
 687             break;
 688           }
 689         }
 690         if(!found) error_handler("error from read_mosaic: the first tile name specified "
 691                                  "in contact is not found in tile list");
 692         found = 0;
 693         for(m=0; m<ntiles; m++) {
 694           if(strcmp(gridtiles[m], pstring[3]) == 0) { /*found the tile name */
 695             found = 1;
 696             tile2[n] = m+1;
 697             break;
 698           }
 699         }
 700         if(!found) error_handler("error from read_mosaic: the second tile name specified "
 701                                  "in contact is not found in tile list");
 702         get_string_data_level(mosaic_file, "contact_index", contacts, &n);
 703         /* parse the string to get contact index */
 704         tokenize( contacts, ":,", STRING, MAXVAR, (char *)pstring, &nstr);
 705         if(nstr != 8) error_handler("Error from read_mosaic: number of elements "
 706                                     "in contact_index seperated by :/, should be 8");
 707         /* make sure the string is only composed of numbers */
 708         for(m=0; m<nstr; m++){
 709           for(l=0; l<strlen(pstring[m]); l++){
 710             if(pstring[m][l] > '9' ||  pstring[m][l] < '0' ) {
 711               error_handler("Error from read_mosaic: some of the character in "
 712                             "contact_indices except token is not digit number");
 713             }
 714           }
 715         }
 716         istart1[n] = atoi(pstring[0]);
 717         iend1[n]   = atoi(pstring[1]);
 718         jstart1[n] = atoi(pstring[2]);
 719         jend1[n]   = atoi(pstring[3]);
 720         istart2[n] = atoi(pstring[4]);
 721         iend2[n]   = atoi(pstring[5]);
 722         jstart2[n] = atoi(pstring[6]);
 723         jend2[n]   = atoi(pstring[7]);
 724         i1_type = transfer_to_model_index(istart1[n], iend1[n], istart1+n, iend1+n, x_refine);
 725         j1_type = transfer_to_model_index(jstart1[n], jend1[n], jstart1+n, jend1+n, y_refine);
 726         i2_type = transfer_to_model_index(istart2[n], iend2[n], istart2+n, iend2+n, x_refine);
 727         j2_type = transfer_to_model_index(jstart2[n], jend2[n], jstart2+n, jend2+n, y_refine);
 728         if( i1_type == 0 && j1_type == 0 )
 729           error_handler("Error from read_mosaic_contact:istart1==iend1 and jstart1==jend1");
 730         if( i2_type == 0 && j2_type == 0 )
 731           error_handler("Error from read_mosaic_contact:istart2==iend2 and jstart2==jend2");
 732         if( i1_type + j1_type != i2_type + j2_type )
 733           error_handler("Error from read_mosaic_contact: It is not a line or overlap contact");
 734
 735       }
 736
 737       for(m=0; m<ntiles; m++) {
 738         free(gridtiles[m]);
 739       }
 740
 741       free(gridtiles);
 742
 743
 744     } /* read_mosaic_contact */
 745
 746
 747     /******************************************************************************
 748   void read_mosaic_grid_data(const char *mosaic_file, const char *name, int nx, int ny,
 749                              double *data, int level, int ioff, int joff)
 750   read mosaic grid information onto model grid. We assume the refinement is 2 right now.
 751   We may remove this restriction in the future. nx and ny are model grid size. level
 752   is the tile number. ioff and joff to indicate grid location. ioff =0 and joff = 0
 753   for C-cell. ioff=0 and joff=1 for E-cell, ioff=1 and joff=0 for N-cell,
 754   ioff=1 and joff=1 for T-cell
 755     ******************************************************************************/
 756     void read_mosaic_grid_data(const char *mosaic_file, const char *name, int nx, int ny,
 757                                double *data, unsigned int level, int ioff, int joff)
 758     {
 759       char   tilefile[STRING], gridfile[STRING], dir[STRING];
 760       double *tmp;
 761       int    ni, nj, nxp, nyp, i, j;
 762
 763       get_file_dir(mosaic_file, dir);
 764
 765       get_string_data_level(mosaic_file, "gridfiles", gridfile, &level);
 766       sprintf(tilefile, "%s/%s", dir, gridfile);
 767
 768       ni = get_dimlen(tilefile, "nx");
 769       nj = get_dimlen(tilefile, "ny");
 770
 771       if( ni != nx*2 || nj != ny*2) error_handler("supergrid size should be double of the model grid size");
 772       tmp = (double *)malloc((ni+1)*(nj+1)*sizeof(double));
 773       get_var_data( tilefile, name, tmp);
 774       nxp = nx + 1 - ioff;
 775       nyp = ny + 1 - joff;
 776       for(j=0; j<nyp; j++) for(i=0; i<nxp; i++) data[j*nxp+i] = tmp[(2*j+joff)*(ni+1)+2*i+ioff];
 777       free(tmp);
 778
 779     } /* read_mosaic_grid_data */