/*********************************************************************** * GNU Lesser General Public License * * This file is part of the GFDL Flexible Modeling System (FMS). * * FMS is free software: you can redistribute it and/or modify it under * the terms of the GNU Lesser General Public License as published by * the Free Software Foundation, either version 3 of the License, or (at * your option) any later version. * * FMS is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License * for more details. * * You should have received a copy of the GNU Lesser General Public * License along with FMS. If not, see <http://www.gnu.org/licenses/>. **********************************************************************/ #include <stdlib.h> #include <stdio.h> #include <math.h> #include <string.h> #include "read_mosaic.h" #include "constant.h" #include "mosaic_util.h" #ifdef use_netCDF #include <netcdf.h> #endif /********************************************************************* void netcdf_error( int status ) status is the returning value of netcdf call. this routine will handle the error when status is not NC_NOERR. ********************************************************************/ void handle_netcdf_error(const char *msg, int status ) { char errmsg[512]; sprintf( errmsg, "%s: %s", msg, (char *)nc_strerror(status) ); error_handler(errmsg); } /* handle_netcdf_error */ /*************************************************************************** void get_file_dir(const char *file, char *dir) get the directory where file is located. The dir will be the complate path before the last "/". If no "/" exist in file, the path will be current ".". ***************************************************************************/ void get_file_dir(const char *file, char *dir) { int len; char *strptr = NULL; /* get the diretory */ strptr = strrchr(file, '/'); if(strptr) { len = strptr - file; strncpy(dir, file, len); } else { len = 1; strcpy(dir, "."); } dir[len] = 0; } /* get_file_dir */ int field_exist(const char* file, const char *name) { int ncid, varid, status; char msg[512]; int existed=0; #ifdef use_netCDF status = nc_open(file, NC_NOWRITE, &ncid); if(status != NC_NOERR) { sprintf(msg, "field_exist: in opening file %s", file); handle_netcdf_error(msg, status); } status = nc_inq_varid(ncid, name, &varid); if(status == NC_NOERR){ existed = 1; } status = nc_close(ncid); if(status != NC_NOERR) { sprintf(msg, "field_exist: in closing file %s.", file); handle_netcdf_error(msg, status); } #else /* ndef use_netCDF */ error_handler("read_mosaic: Add flag -Duse_netCDF when compiling"); #endif /* use_netcdf */ return existed; } /* field_exist */ int get_dimlen(const char* file, const char *name) { int ncid, dimid, status, len; size_t size; char msg[512]; len = 0; #ifdef use_netCDF status = nc_open(file, NC_NOWRITE, &ncid); if(status != NC_NOERR) { sprintf(msg, "in opening file %s", file); handle_netcdf_error(msg, status); } status = nc_inq_dimid(ncid, name, &dimid); if(status != NC_NOERR) { sprintf(msg, "in getting dimid of %s from file %s.", name, file); handle_netcdf_error(msg, status); } status = nc_inq_dimlen(ncid, dimid, &size); if(status != NC_NOERR) { sprintf(msg, "in getting dimension size of %s from file %s.", name, file); handle_netcdf_error(msg, status); } status = nc_close(ncid); if(status != NC_NOERR) { sprintf(msg, "in closing file %s.", file); handle_netcdf_error(msg, status); } len = size; if(status != NC_NOERR) { sprintf(msg, "in closing file %s", file); handle_netcdf_error(msg, status); } #else error_handler("read_mosaic: Add flag -Duse_netCDF when compiling"); #endif return len; } /* get_dimlen */ /******************************************************************************* void get_string_data(const char *file, const char *name, char *data) get string data of field with "name" from "file". ******************************************************************************/ void get_string_data(const char *file, const char *name, char *data) { int ncid, varid, status; char msg[512]; #ifdef use_netCDF status = nc_open(file, NC_NOWRITE, &ncid); if(status != NC_NOERR) { sprintf(msg, "in opening file %s", file); handle_netcdf_error(msg, status); } status = nc_inq_varid(ncid, name, &varid); if(status != NC_NOERR) { sprintf(msg, "in getting varid of %s from file %s.", name, file); handle_netcdf_error(msg, status); } status = nc_get_var_text(ncid, varid, data); if(status != NC_NOERR) { sprintf(msg, "in getting data of %s from file %s.", name, file); handle_netcdf_error(msg, status); } status = nc_close(ncid); if(status != NC_NOERR) { sprintf(msg, "in closing file %s.", file); handle_netcdf_error(msg, status); } #else error_handler("read_mosaic: Add flag -Duse_netCDF when compiling"); #endif } /* get_string_data */ /******************************************************************************* void get_string_data_level(const char *file, const char *name, const size_t *start, const size_t *nread, char *data) get string data of field with "name" from "file". ******************************************************************************/ void get_string_data_level(const char *file, const char *name, char *data, const unsigned int *level) { int ncid, varid, status, i; size_t start[4], nread[4]; char msg[512]; #ifdef use_netCDF status = nc_open(file, NC_NOWRITE, &ncid); if(status != NC_NOERR) { sprintf(msg, "in opening file %s", file); handle_netcdf_error(msg, status); } status = nc_inq_varid(ncid, name, &varid); if(status != NC_NOERR) { sprintf(msg, "in getting varid of %s from file %s.", name, file); handle_netcdf_error(msg, status); } for(i=0; i<4; i++) { start[i] = 0; nread[i] = 1; } start[0] = *level; nread[1] = STRING; status = nc_get_vara_text(ncid, varid, start, nread, data); if(status != NC_NOERR) { sprintf(msg, "in getting data of %s from file %s.", name, file); handle_netcdf_error(msg, status); } status = nc_close(ncid); if(status != NC_NOERR) { sprintf(msg, "in closing file %s.", file); handle_netcdf_error(msg, status); } #else error_handler("read_mosaic: Add flag -Duse_netCDF when compiling"); #endif } /* get_string_data_level */ /******************************************************************************* void get_var_data(const char *file, const char *name, double *data) get var data of field with "name" from "file". ******************************************************************************/ void get_var_data(const char *file, const char *name, void *data) { int ncid, varid, status; nc_type vartype; char msg[512]; #ifdef use_netCDF status = nc_open(file, NC_NOWRITE, &ncid); if(status != NC_NOERR) { sprintf(msg, "in opening file %s", file); handle_netcdf_error(msg, status); } status = nc_inq_varid(ncid, name, &varid); if(status != NC_NOERR) { sprintf(msg, "in getting varid of %s from file %s.", name, file); handle_netcdf_error(msg, status); } status = nc_inq_vartype(ncid, varid, &vartype); if(status != NC_NOERR) { sprintf(msg, "get_var_data: in getting vartype of of %s in file %s ", name, file); handle_netcdf_error(msg, status); } switch (vartype) { case NC_DOUBLE:case NC_FLOAT: #ifdef OVERLOAD_R4 status = nc_get_var_float(ncid, varid, data); #else status = nc_get_var_double(ncid, varid, data); #endif break; case NC_INT: status = nc_get_var_int(ncid, varid, data); break; default: sprintf(msg, "get_var_data: field %s in file %s has an invalid type, " "the type should be NC_DOUBLE, NC_FLOAT or NC_INT", name, file); error_handler(msg); } if(status != NC_NOERR) { sprintf(msg, "in getting data of %s from file %s.", name, file); handle_netcdf_error(msg, status); } status = nc_close(ncid); if(status != NC_NOERR) { sprintf(msg, "in closing file %s.", file); handle_netcdf_error(msg, status); } #else error_handler("read_mosaic: Add flag -Duse_netCDF when compiling"); #endif } /* get_var_data */ /******************************************************************************* void get_var_data(const char *file, const char *name, double *data) get var data of field with "name" from "file". ******************************************************************************/ void get_var_data_region(const char *file, const char *name, const size_t *start, const size_t *nread, void *data) { int ncid, varid, status; nc_type vartype; char msg[512]; #ifdef use_netCDF status = nc_open(file, NC_NOWRITE, &ncid); if(status != NC_NOERR) { sprintf(msg, "get_var_data_region: in opening file %s", file); handle_netcdf_error(msg, status); } status = nc_inq_varid(ncid, name, &varid); if(status != NC_NOERR) { sprintf(msg, "in getting varid of %s from file %s.", name, file); handle_netcdf_error(msg, status); } status = nc_inq_vartype(ncid, varid, &vartype); if(status != NC_NOERR) { sprintf(msg, "get_var_data_region: in getting vartype of of %s in file %s ", name, file); handle_netcdf_error(msg, status); } switch (vartype) { case NC_DOUBLE:case NC_FLOAT: #ifdef OVERLOAD_R4 status = nc_get_vara_float(ncid, varid, start, nread, data); #else status = nc_get_vara_double(ncid, varid, start, nread, data); #endif break; case NC_INT: status = nc_get_vara_int(ncid, varid, start, nread, data); break; default: sprintf(msg, "get_var_data_region: field %s in file %s has an invalid type, " "the type should be NC_DOUBLE, NC_FLOAT or NC_INT", name, file); error_handler(msg); } if(status != NC_NOERR) { sprintf(msg, "get_var_data_region: in getting data of %s from file %s.", name, file); handle_netcdf_error(msg, status); } status = nc_close(ncid); if(status != NC_NOERR) { sprintf(msg, "get_var_data_region: in closing file %s.", file); handle_netcdf_error(msg, status); } #else error_handler("read_mosaic: Add flag -Duse_netCDF when compiling"); #endif } /* get_var_data_region */ /****************************************************************************** void get_var_text_att(const char *file, const char *name, const char *attname, char *att) get text attribute of field 'name' from 'file ******************************************************************************/ void get_var_text_att(const char *file, const char *name, const char *attname, char *att) { int ncid, varid, status; char msg[512]; #ifdef use_netCDF status = nc_open(file, NC_NOWRITE, &ncid); if(status != NC_NOERR) { sprintf(msg, "in opening file %s", file); handle_netcdf_error(msg, status); } status = nc_inq_varid(ncid, name, &varid); if(status != NC_NOERR) { sprintf(msg, "in getting varid of %s from file %s.", name, file); handle_netcdf_error(msg, status); } status = nc_get_att_text(ncid, varid, attname, att); if(status != NC_NOERR) { sprintf(msg, "in getting attribute %s of %s from file %s.", attname, name, file); handle_netcdf_error(msg, status); } status = nc_close(ncid); if(status != NC_NOERR) { sprintf(msg, "in closing file %s.", file); handle_netcdf_error(msg, status); } #else error_handler("read_mosaic: Add flag -Duse_netCDF when compiling"); #endif } /* get_var_text_att */ /*********************************************************************** return number of overlapping cells. ***********************************************************************/ #ifndef __AIX int read_mosaic_xgrid_size_( const char *xgrid_file ) { return read_mosaic_xgrid_size(xgrid_file); } #endif int read_mosaic_xgrid_size( const char *xgrid_file ) { int ncells; ncells = get_dimlen(xgrid_file, "ncells"); return ncells; } #ifdef OVERLOAD_R4 float get_global_area(void) { float garea; #else double get_global_area(void) { double garea; #endif garea = 4*M_PI*RADIUS*RADIUS; return garea; } #ifndef __AIX #ifdef OVERLOAD_R4 float get_global_area_(void) { float garea; #else double get_global_area_(void) { double garea; #endif garea = 4*M_PI*RADIUS*RADIUS; return garea; } #endif /****************************************************************************/ #ifndef __AIX #ifdef OVERLOAD_R4 void read_mosaic_xgrid_order1_(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, float *area ) #else void read_mosaic_xgrid_order1_(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, double *area ) #endif { read_mosaic_xgrid_order1(xgrid_file, i1, j1, i2, j2, area); } #endif #ifdef OVERLOAD_R4 void read_mosaic_xgrid_order1(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, float *area ) #else void read_mosaic_xgrid_order1(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, double *area ) #endif { int ncells, n; int *tile1_cell, *tile2_cell; #ifdef OVERLOAD_R4 float garea; #else double garea; #endif ncells = get_dimlen(xgrid_file, "ncells"); tile1_cell = (int *)malloc(ncells*2*sizeof(int)); tile2_cell = (int *)malloc(ncells*2*sizeof(int)); get_var_data(xgrid_file, "tile1_cell", tile1_cell); get_var_data(xgrid_file, "tile2_cell", tile2_cell); get_var_data(xgrid_file, "xgrid_area", area); garea = 4*M_PI*RADIUS*RADIUS; for(n=0; n<ncells; n++) { i1[n] = tile1_cell[n*2] - 1; j1[n] = tile1_cell[n*2+1] - 1; i2[n] = tile2_cell[n*2] - 1; j2[n] = tile2_cell[n*2+1] - 1; area[n] /= garea; /* rescale the exchange grid area to unit earth area */ } free(tile1_cell); free(tile2_cell); } /* read_mosaic_xgrid_order1 */ #ifndef __AIX #ifdef OVERLOAD_R4 void read_mosaic_xgrid_order1_region_(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, float *area, int *isc, int *iec ) #else void read_mosaic_xgrid_order1_region_(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, double *area, int *isc, int *iec ) #endif { read_mosaic_xgrid_order1_region(xgrid_file, i1, j1, i2, j2, area, isc, iec); } #endif #ifdef OVERLOAD_R4 void read_mosaic_xgrid_order1_region(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, float *area, int *isc, int *iec ) #else void read_mosaic_xgrid_order1_region(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, double *area, int *isc, int *iec ) #endif { int ncells, n, i; int *tile1_cell, *tile2_cell; size_t start[4], nread[4]; #ifdef OVERLOAD_R4 float garea; #else double garea; #endif ncells = *iec-*isc+1; tile1_cell = (int *)malloc(ncells*2*sizeof(int)); tile2_cell = (int *)malloc(ncells*2*sizeof(int)); for(i=0; i<4; i++) { start[i] = 0; nread[i] = 1; } start[0] = *isc; nread[0] = ncells; nread[1] = 2; get_var_data_region(xgrid_file, "tile1_cell", start, nread, tile1_cell); get_var_data_region(xgrid_file, "tile2_cell", start, nread, tile2_cell); nread[1] = 1; get_var_data_region(xgrid_file, "xgrid_area", start, nread, area); garea = 4*M_PI*RADIUS*RADIUS; for(n=0; n<ncells; n++) { i1[n] = tile1_cell[n*2] - 1; j1[n] = tile1_cell[n*2+1] - 1; i2[n] = tile2_cell[n*2] - 1; j2[n] = tile2_cell[n*2+1] - 1; area[n] /= garea; /* rescale the exchange grid area to unit earth area */ } free(tile1_cell); free(tile2_cell); } /* read_mosaic_xgrid_order1 */ /* NOTE: di, dj is for tile1, */ /****************************************************************************/ #ifndef __AIX #ifdef OVERLOAD_R4 void read_mosaic_xgrid_order2_(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, float *area, float *di, float *dj ) #else void read_mosaic_xgrid_order2_(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, double *area, double *di, double *dj ) #endif { read_mosaic_xgrid_order2(xgrid_file, i1, j1, i2, j2, area, di, dj); } #endif #ifdef OVERLOAD_R4 void read_mosaic_xgrid_order2(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, float *area, float *di, float *dj ) #else void read_mosaic_xgrid_order2(const char *xgrid_file, int *i1, int *j1, int *i2, int *j2, double *area, double *di, double *dj ) #endif { int ncells, n; int *tile1_cell, *tile2_cell; double *tile1_distance; #ifdef OVERLOAD_R4 float garea; #else double garea; #endif ncells = get_dimlen(xgrid_file, "ncells"); tile1_cell = (int *)malloc(ncells*2*sizeof(int )); tile2_cell = (int *)malloc(ncells*2*sizeof(int )); tile1_distance = (double *)malloc(ncells*2*sizeof(double)); get_var_data(xgrid_file, "tile1_cell", tile1_cell); get_var_data(xgrid_file, "tile2_cell", tile2_cell); get_var_data(xgrid_file, "xgrid_area", area); get_var_data(xgrid_file, "tile1_distance", tile1_distance); garea = 4*M_PI*RADIUS*RADIUS; for(n=0; n<ncells; n++) { i1[n] = tile1_cell[n*2] - 1; j1[n] = tile1_cell[n*2+1] - 1; i2[n] = tile2_cell[n*2] - 1; j2[n] = tile2_cell[n*2+1] - 1; di[n] = tile1_distance[n*2]; dj[n] = tile1_distance[n*2+1]; area[n] /= garea; /* rescale the exchange grid area to unit earth area */ } free(tile1_cell); free(tile2_cell); free(tile1_distance); } /* read_mosaic_xgrid_order2 */ /****************************************************************************** int read_mosaic_ntiles(const char *mosaic_file) return number tiles in mosaic_file ******************************************************************************/ #ifndef __AIX int read_mosaic_ntiles_(const char *mosaic_file) { return read_mosaic_ntiles(mosaic_file); } #endif int read_mosaic_ntiles(const char *mosaic_file) { int ntiles; ntiles = get_dimlen(mosaic_file, "ntiles"); return ntiles; } /* read_mosaic_ntiles */ /****************************************************************************** int read_mosaic_ncontacts(const char *mosaic_file) return number of contacts in mosaic_file ******************************************************************************/ #ifndef __AIX int read_mosaic_ncontacts_(const char *mosaic_file) { return read_mosaic_ncontacts(mosaic_file); } #endif int read_mosaic_ncontacts(const char *mosaic_file) { int ncontacts; if(field_exist(mosaic_file, "contacts") ) ncontacts = get_dimlen(mosaic_file, "ncontact"); else ncontacts = 0; return ncontacts; } /* read_mosaic_ncontacts */ /***************************************************************************** void read_mosaic_grid_sizes(const char *mosaic_file, int *nx, int *ny) read mosaic grid size of each tile, currently we are assuming the refinement is 2. We assume the grid files are located at the same directory as mosaic_file. *****************************************************************************/ #ifndef __AIX void read_mosaic_grid_sizes_(const char *mosaic_file, int *nx, int *ny) { read_mosaic_grid_sizes(mosaic_file, nx, ny); } #endif void read_mosaic_grid_sizes(const char *mosaic_file, int *nx, int *ny) { unsigned int ntiles, n; char gridfile[STRING], tilefile[2*STRING]; char dir[STRING]; const int x_refine = 2, y_refine = 2; get_file_dir(mosaic_file, dir); ntiles = get_dimlen(mosaic_file, "ntiles"); for(n = 0; n < ntiles; n++) { get_string_data_level(mosaic_file, "gridfiles", gridfile, &n); sprintf(tilefile, "%s/%s", dir, gridfile); nx[n] = get_dimlen(tilefile, "nx"); ny[n] = get_dimlen(tilefile, "ny"); if(nx[n]%x_refine != 0) error_handler("Error from read_mosaic_grid_sizes: nx is not divided by x_refine"); if(ny[n]%y_refine != 0) error_handler("Error from read_mosaic_grid_sizes: ny is not divided by y_refine"); nx[n] /= x_refine; ny[n] /= y_refine; } } /* read_mosaic_grid_sizes */ /****************************************************************************** void read_mosaic_contact(const char *mosaic_file) read mosaic contact information ******************************************************************************/ #ifndef __AIX void read_mosaic_contact_(const char *mosaic_file, int *tile1, int *tile2, int *istart1, int *iend1, int *jstart1, int *jend1, int *istart2, int *iend2, int *jstart2, int *jend2) { read_mosaic_contact(mosaic_file, tile1, tile2, istart1, iend1, jstart1, jend1, istart2, iend2, jstart2, jend2); } #endif /* transfer the index from supergrid to model grid return 0 if istart = iend otherwise return 1. */ int transfer_to_model_index(int istart_in, int iend_in, int *istart_out, int *iend_out, int refine_ratio) { int type; if( istart_in == iend_in ) { type = 0; istart_out[0] = (istart_in+1)/refine_ratio-1; iend_out[0] = istart_out[0]; } else { type = 1; if( iend_in > istart_in ) { istart_out[0] = istart_in - 1; iend_out[0] = iend_in - refine_ratio; } else { istart_out[0] = istart_in - refine_ratio; iend_out[0] = iend_in - 1; } if( istart_out[0]%refine_ratio || iend_out[0]%refine_ratio) error_handler("Error from read_mosaic: mismatch between refine_ratio and istart_in/iend_in"); istart_out[0] /= refine_ratio; iend_out[0] /= refine_ratio; } return type; } void read_mosaic_contact(const char *mosaic_file, int *tile1, int *tile2, int *istart1, int *iend1, int *jstart1, int *jend1, int *istart2, int *iend2, int *jstart2, int *jend2) { char contacts[STRING]; char **gridtiles; #define MAXVAR 40 char pstring[MAXVAR][STRING]; unsigned int nstr, ntiles, ncontacts, n, m, l, found; const int x_refine = 2, y_refine = 2; int i1_type, j1_type, i2_type, j2_type; ntiles = get_dimlen(mosaic_file, "ntiles"); gridtiles = (char **)malloc(ntiles*sizeof(char *)); for(n=0; n<ntiles; n++) { gridtiles[n] = (char *)malloc(STRING*sizeof(char)); get_string_data_level(mosaic_file, "gridtiles", gridtiles[n], &n); } ncontacts = get_dimlen(mosaic_file, "ncontact"); for(n = 0; n < ncontacts; n++) { get_string_data_level(mosaic_file, "contacts", contacts, &n); /* parse the string contacts to get tile number */ tokenize( contacts, ":", STRING, MAXVAR, (char *)pstring, &nstr); if(nstr != 4) error_handler("Error from read_mosaic: number of elements " "in contact seperated by :/:: should be 4"); found = 0; for(m=0; m<ntiles; m++) { if(strcmp(gridtiles[m], pstring[1]) == 0) { /*found the tile name */ found = 1; tile1[n] = m+1; break; } } if(!found) error_handler("error from read_mosaic: the first tile name specified " "in contact is not found in tile list"); found = 0; for(m=0; m<ntiles; m++) { if(strcmp(gridtiles[m], pstring[3]) == 0) { /*found the tile name */ found = 1; tile2[n] = m+1; break; } } if(!found) error_handler("error from read_mosaic: the second tile name specified " "in contact is not found in tile list"); get_string_data_level(mosaic_file, "contact_index", contacts, &n); /* parse the string to get contact index */ tokenize( contacts, ":,", STRING, MAXVAR, (char *)pstring, &nstr); if(nstr != 8) error_handler("Error from read_mosaic: number of elements " "in contact_index seperated by :/, should be 8"); /* make sure the string is only composed of numbers */ for(m=0; m<nstr; m++){ for(l=0; l<strlen(pstring[m]); l++){ if(pstring[m][l] > '9' || pstring[m][l] < '0' ) { error_handler("Error from read_mosaic: some of the character in " "contact_indices except token is not digit number"); } } } istart1[n] = atoi(pstring[0]); iend1[n] = atoi(pstring[1]); jstart1[n] = atoi(pstring[2]); jend1[n] = atoi(pstring[3]); istart2[n] = atoi(pstring[4]); iend2[n] = atoi(pstring[5]); jstart2[n] = atoi(pstring[6]); jend2[n] = atoi(pstring[7]); i1_type = transfer_to_model_index(istart1[n], iend1[n], istart1+n, iend1+n, x_refine); j1_type = transfer_to_model_index(jstart1[n], jend1[n], jstart1+n, jend1+n, y_refine); i2_type = transfer_to_model_index(istart2[n], iend2[n], istart2+n, iend2+n, x_refine); j2_type = transfer_to_model_index(jstart2[n], jend2[n], jstart2+n, jend2+n, y_refine); if( i1_type == 0 && j1_type == 0 ) error_handler("Error from read_mosaic_contact:istart1==iend1 and jstart1==jend1"); if( i2_type == 0 && j2_type == 0 ) error_handler("Error from read_mosaic_contact:istart2==iend2 and jstart2==jend2"); if( i1_type + j1_type != i2_type + j2_type ) error_handler("Error from read_mosaic_contact: It is not a line or overlap contact"); } for(m=0; m<ntiles; m++) { free(gridtiles[m]); } free(gridtiles); } /* read_mosaic_contact */ /****************************************************************************** void read_mosaic_grid_data(const char *mosaic_file, const char *name, int nx, int ny, double *data, int level, int ioff, int joff) read mosaic grid information onto model grid. We assume the refinement is 2 right now. We may remove this restriction in the future. nx and ny are model grid size. level is the tile number. ioff and joff to indicate grid location. ioff =0 and joff = 0 for C-cell. ioff=0 and joff=1 for E-cell, ioff=1 and joff=0 for N-cell, ioff=1 and joff=1 for T-cell ******************************************************************************/ void read_mosaic_grid_data(const char *mosaic_file, const char *name, int nx, int ny, double *data, unsigned int level, int ioff, int joff) { char tilefile[STRING], gridfile[STRING], dir[STRING]; double *tmp; int ni, nj, nxp, nyp, i, j; get_file_dir(mosaic_file, dir); get_string_data_level(mosaic_file, "gridfiles", gridfile, &level); sprintf(tilefile, "%s/%s", dir, gridfile); ni = get_dimlen(tilefile, "nx"); nj = get_dimlen(tilefile, "ny"); if( ni != nx*2 || nj != ny*2) error_handler("supergrid size should be double of the model grid size"); tmp = (double *)malloc((ni+1)*(nj+1)*sizeof(double)); get_var_data( tilefile, name, tmp); nxp = nx + 1 - ioff; nyp = ny + 1 - joff; 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]; free(tmp); } /* read_mosaic_grid_data */