#include #include #include #include #include "constant.h" #include "mosaic_util.h" #include "tool_util.h" #include "interp.h" #include "mpp.h" #include "mpp_domain.h" #include "mpp_io.h" #include "conserve_interp.h" #include "bilinear_interp.h" #include "fregrid_util.h" #define D2R (M_PI/180.) #define R2D (180./M_PI) const double SMALL = 1.0e-4; double distant(double a, double b, double met1, double met2); double bp_lam(double x, double y, double bpeq, double rp); double bp_phi(double x, double y, double bpsp, double bpnp); double lon_in_range(double lon, double lon_strt); void vtx_insert(double *x, double *y, int *n, int n_ins); void vtx_delete(double *x, double *y, int *n, int n_del); int lon_fix(double *x, double *y, int n_in, double tlon); /*************************************************************************** void get_file_path(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_path(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_path */ int get_int_entry(char *line, int* value) { char* pch; int num; pch = strtok(line, ", "); num = 0; while( pch != NULL) { value[num++] = atoi(pch); pch = strtok(NULL, ", "); } return num; }; int get_double_entry(char *line, double *value) { char* pch; int num; pch = strtok(line, ", "); num = 0; while( pch != NULL) { value[num++] = atof(pch); pch = strtok(NULL, ", "); } return num; }; /********************************************************************* double spherical_dist(double x1, double y1, double x2, double y2) return distance between spherical grid on the earth *********************************************************************/ double spherical_dist(double x1, double y1, double x2, double y2) { double dist = 0.0; double h1, h2; if(x1 == x2) { h1 = RADIUS; h2 = RADIUS; dist = distant(y1,y2,h1,h2); } else if(y1 == y2) { h1 = RADIUS * cos(y1*D2R); h2 = RADIUS * cos(y2*D2R); dist = distant(x1,x2,h1,h2); } else mpp_error("tool_till: This is not rectangular grid"); return dist; }; /* spherical_dist */ /********************************************************************* void double bipolar_dist(double x1, double y1, double x2, double y2) return distance of bipolar grids *********************************************************************/ double bipolar_dist(double x1, double y1, double x2, double y2, double bpeq, double bpsp, double bpnp, double rp ) { double dist, x[2],y[2], bp_lon[2], bp_lat[2], metric[2]; double h1[2], h2[2], chic; int n; x[0] = x1; x[1] = x2; y[0] = y1; y[1] = y2; /*--- get the bipolar grid and metric term ----------------------------*/ for(n=0; n<2; n++){ bp_lon[n] = bp_lam(x[n],y[n],bpeq, rp); /* longitude (degrees) in bipolar grid system */ bp_lat[n] = bp_phi(x[n],y[n],bpsp, bpnp); /* latitude (degrees) in bipolar grid system */ h1[n] = RADIUS*cos(bp_lat[n]*D2R); h2[n] = RADIUS; metric[n] = 1.0; if (fabs(y[n]-90.0) < SMALL || fabs(bp_lon[n]*D2R) >= SMALL || fabs(bp_lat[n]*D2R) >= SMALL) { chic = acos(cos(bp_lon[n]*D2R)*cos(bp_lat[n]*D2R)); /* eqn. 6 */ metric[n] = rp*(1/pow(cos(chic/2),2))/(1+(pow(rp,2))*(pow(tan(chic/2),2)));/* eq 3 */ } } /*--- then calculate the distance -------------------------------------*/ if(x1 == x2) dist = distant(bp_lon[0],bp_lon[1],metric[0]*h1[0],metric[1]*h1[1]); else if(y1 == y2) dist = distant(bp_lat[0],bp_lat[1],metric[0]*h2[0],metric[1]*h2[1]); else mpp_error("tool_util: This tripolar grid not transformed from rectangular grid"); return dist; }; /* bipolar_dist */ /********************************************************************* double distant(double a, double b, double met1, double met2) return distant on the earth *********************************************************************/ double distant(double a, double b, double met1, double met2) { return fabs(a-b)*D2R*(met1+met2)/2. ; }; /* distant */ /********************************************************************* double spherical_area(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4 ) rectangular grid box area ********************************************************************/ double spherical_area(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4 ) { double area, dx, lat1, lat2, x[4],y[4]; int i, ip; x[0] = x1; y[0] = y1; x[1] = x2; y[1] = y2; x[2] = x3; y[2] = y3; x[3] = x4; y[3] = y4; area = 0.0; for(i=0; i<4; i++) { ip = i+1; if(ip ==4) ip = 0; dx = (x[ip] - x[i])*D2R; lat1 = y[ip]*D2R; lat2 = y[i]*D2R; if(dx==0.0) continue; if(dx > M_PI) dx = dx - 2.0*M_PI; if(dx < -M_PI) dx = dx + 2.0*M_PI; if (lat1 == lat2) /* cheap area calculation along latitude */ area = area - dx*sin(lat1); else area = area - dx*(sin(lat1)+sin(lat2))/2; /* TRAPEZOID_RULE */ } area = area * RADIUS * RADIUS; return area; }; /* spherical_area */ /********************************************************************* double bipolar_area(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4 ) bipolar grid area ********************************************************************/ double bipolar_area(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4 ) { double area, dx, lat1, lat2, x[8],y[8]; int i, ip, n; x[0] = x1; y[0] = y1; x[1] = x2; y[1] = y2; x[2] = x3; y[2] = y3; x[3] = x4; y[3] = y4; /*--- first fix the longitude at the pole -----------------------------*/ n = lon_fix(x, y, 4, 180.); /*--- calculate the area ---------------------------------------------- */ area = 0.0; for(i=0; i M_PI) dx = dx - 2.0*M_PI; if(dx < -M_PI) dx = dx + 2.0*M_PI; if (lat1 == lat2) /* cheap area calculation along latitude */ area = area - dx*sin(lat1); else area = area - dx*(sin(lat1)+sin(lat2))/2; /* TRAPEZOID_RULE */ } area = area * RADIUS * RADIUS; return area; }; /* bipolar_area */ /********************************************************************* double lat_dist(double x1, double x2) distance (in degrees) between points on lat. circle ********************************************************************/ double lat_dist(double x1, double x2) { return min(fmod(x1-x2+720,360.),fmod(x2-x1+720,360.)); }; /********************************************************************* double bp_lam(double x, double y, double bpeq) find bipolar grid longitude given geo. coordinates ********************************************************************/ double bp_lam(double x, double y, double bpeq, double rp) { double bp_lam; /* bp_lam = ((90-y)/(90-lat_join))*90 */ /* invert eqn. 5 with phic=0 to place point at specified geo. lat */ bp_lam = 2.*atan(tan((0.5*M_PI-y*D2R)/2)/rp)*R2D; if (lat_dist(x,bpeq)<90.) bp_lam = -bp_lam; return bp_lam; }; /* bp_lam */ /********************************************************************* double bp_phi(double x, double y, double bpsp, double bpnp) find bipolar grid latitude given geo. coordinates ********************************************************************/ double bp_phi(double x, double y, double bpsp, double bpnp) { if (lat_dist(x,bpsp)<90.) return (-90+lat_dist(x,bpsp)); else return ( 90-lat_dist(x,bpnp)); }; /* bp_phi */ /********************************************************************* void tp_trans(double& lon, double& lat, double lon_ref) calculate tripolar grid ********************************************************************/ void tp_trans(double *lon, double *lat, double lon_ref, double lon_start, double lam0, double bpeq, double bpsp, double bpnp, double rp ) { double lamc, phic, lams, chic, phis; lamc = bp_lam(*lon, *lat, bpeq, rp )*D2R; phic = bp_phi(*lon, *lat, bpsp, bpnp)*D2R; if (fabs(*lat-90.) < SMALL) { if (phic > 0) *lon=lon_in_range(lon_start,lon_ref); else *lon=lon_start+180.; chic = acos(cos(lamc)*cos(phic)); /* eqn. 6 */ phis = M_PI*0.5-2*atan(rp*tan(chic/2)); /* eqn. 5 */ *lat = phis*R2D; return; } if (fabs(lamc) < SMALL && fabs(phic) < SMALL) { *lat=90.; *lon=lon_ref; } else { lams = fmod(lam0+M_PI+M_PI/2-atan2(sin(lamc),tan(phic)),2*M_PI); /* eqn. 5 */ chic = acos(cos(lamc)*cos(phic)); /* eqn. 6 */ phis = M_PI*0.5-2*atan(rp*tan(chic/2)); /* eqn. 5 */ *lon = lams*R2D; *lon = lon_in_range(*lon,lon_ref); *lat = phis*R2D; } }; /* tp_trans */ /********************************************************************* double Lon_in_range(double lon, double lon_strt) Returns lon_strt <= longitude <= lon_strt+360 ********************************************************************/ double lon_in_range(double lon, double lon_strt) { double lon_in_range, lon_end; lon_in_range = lon; lon_end = lon_strt+360.; if (fabs(lon_in_range - lon_strt) < SMALL) lon_in_range = lon_strt; else if (fabs(lon_in_range - lon_end) < SMALL) lon_in_range = lon_strt; else { while(1) { if (lon_in_range < lon_strt) lon_in_range = lon_in_range + 360.; else if (lon_in_range > lon_end) lon_in_range = lon_in_range - 360.; else break; } } return lon_in_range; }; /* lon_in_range */ /********************************************************************* int lon_fix(double *x, double *y, int n_in, double tlon) fix longitude at pole. ********************************************************************/ int lon_fix(double *x, double *y, int n_in, double tlon) { int i, ip, im, n_out; double x_sum, dx; n_out = n_in; i = 0; while( i < n_out) { if(fabs(y[i]) >= 90.-SMALL) { im = i - 1; if(im < 0) im = im + n_out; ip = i + 1; if(ip >= n_out) ip = ip - n_out; /*--- all pole points must be paired ---------------------------- */ if(y[im] == y[i] && y[ip] == y[i] ) { vtx_delete(x,y, &n_out, i); i = i - 1; } else if(y[im] != y[i] && y[ip] != y[i] ) { vtx_insert(x,y,&n_out,i); i = i + 1; } } i = i + 1; } /*--- first of pole pair has longitude of previous vertex ------------- --- second of pole pair has longitude of subsequent vertex ---------- */ for(i=0;i= 90.-SMALL) { im= i - 1; if(im < 0) im = im + n_out; ip = i + 1; if(ip >= n_out) ip = ip - n_out; if(y[im] != y[i]) x[i] = x[im]; if(y[ip] != y[i]) x[i] = x[ip]; } } if(n_out == 0) return 0; x_sum = x[1]; for(i=1;i< n_out;i++){ dx = x[i] - x[i-1]; if(dx < -180) dx = dx + 360; else if (dx > 180) dx = dx - 360; x[i] = x[i-1] + dx; x_sum = x_sum + x[i]; } dx = x_sum/(n_out) - tlon; if (dx < -180.) for(i=0;i 180.) for(i=0;i=n_ins; i--){ x[i+1] = x[i]; y[i+1] = y[i]; } (*n)++; }; /* vtx_insert */ /*---------------------------------------------------------------------- void vect_cross(e, p1, p2) Perform cross products of 3D vectors: e = P1 X P2 -------------------------------------------------------------------*/ /******************************************************************************** void compute_grid_bound(int nb, const couble *bnds, const int *npts, int *grid_size, const char *center_cell) compute the 1-D grid location. ********************************************************************************/ double* compute_grid_bound(int nb, const double *bnds, const int *npts, int *grid_size, const char *center) { int refine, i, n, np; double *grid=NULL, *tmp=NULL; double *grid1=NULL, *grid2=NULL; if(!strcmp(center, "none") ) refine = 1; else if(!strcmp(center, "t_cell") || !strcmp(center, "c_cell") ) refine = 2; else mpp_error("tool_util: center should be 'none', 'c_cell' or 't_cell' "); grid1 = (double *)malloc(nb*sizeof(double)); grid1[0] = 1; n = 0; for(i=1; i 1) mpp_error( "make_hgrid: cubic grid generation must be run one processor, contact developer"); } /* Currently we restrict nx can be divided by ndivx and ny can be divided by ndivy */ if(ntilex >0 && ntilex != ntiles) mpp_error("make_hgrid: number of entry specified through --ndivx does not equal ntiles"); if(ntiley >0 && ntiley != ntiles) mpp_error("make_hgrid: number of entry specified through --ndivy does not equal ntiles"); for(n=0; n1) sprintf(outfile, "%s.tile%d.nc", gridname, l); else sprintf(outfile, "%s.nc", gridname); fid = mpp_open(outfile, MPP_WRITE); /* define dimenison */ ni = nx/ndivx[n]; nj = ny/ndivy[n]; nip = ni + 1; njp = nj + 1; dimlist[0] = mpp_def_dim(fid, "string", STRINGLEN); dimlist[1] = mpp_def_dim(fid, "nx", ni); dimlist[2] = mpp_def_dim(fid, "ny", nj); dimlist[3] = mpp_def_dim(fid, "nxp", nip); dimlist[4] = mpp_def_dim(fid, "nyp", njp); /* define variable */ if( strcmp(north_pole_tile, "none") == 0) /* no north pole, then no projection */ id_tile = mpp_def_var(fid, "tile", MPP_CHAR, 1, dimlist, 4, "standard_name", "grid_tile_spec", "geometry", geometry, "discretization", discretization, "conformal", conformal ); else if( strcmp(projection, "none") == 0) id_tile = mpp_def_var(fid, "tile", MPP_CHAR, 1, dimlist, 5, "standard_name", "grid_tile_spec", "geometry", geometry, "north_pole", north_pole_tile, "discretization", discretization, "conformal", conformal ); else id_tile = mpp_def_var(fid, "tile", MPP_CHAR, 1, dimlist, 6, "standard_name", "grid_tile_spec", "geometry", geometry, "north_pole", north_pole_tile, "projection", projection, "discretization", discretization, "conformal", conformal ); dims[0] = dimlist[4]; dims[1] = dimlist[3]; id_x = mpp_def_var(fid, "x", MPP_DOUBLE, 2, dims, 2, "standard_name", "geographic_longitude", "units", "degree_east"); id_y = mpp_def_var(fid, "y", MPP_DOUBLE, 2, dims, 2, "standard_name", "geographic_latitude", "units", "degree_north"); dims[0] = dimlist[4]; dims[1] = dimlist[1]; id_dx = mpp_def_var(fid, "dx", MPP_DOUBLE, 2, dims, 2, "standard_name", "grid_edge_x_distance", "units", "meters"); dims[0] = dimlist[2]; dims[1] = dimlist[3]; id_dy = mpp_def_var(fid, "dy", MPP_DOUBLE, 2, dims, 2, "standard_name", "grid_edge_y_distance", "units", "meters"); dims[0] = dimlist[2]; dims[1] = dimlist[1]; id_area = mpp_def_var(fid, "area", MPP_DOUBLE, 2, dims, 2, "standard_name", "grid_cell_area", "units", "m2" ); dims[0] = dimlist[4]; dims[1] = dimlist[3]; id_angle_dx = mpp_def_var(fid, "angle_dx", MPP_DOUBLE, 2, dims, 2, "standard_name", "grid_vertex_x_angle_WRT_geographic_east", "units", "degrees_east"); if(strcmp(conformal, "true") != 0) id_angle_dy = mpp_def_var(fid, "angle_dy", MPP_DOUBLE, 2, dims, 2, "standard_name", "grid_vertex_y_angle_WRT_geographic_north", "units", "degrees_north"); if( strcmp(north_pole_arcx, "none") == 0) id_arcx = mpp_def_var(fid, "arcx", MPP_CHAR, 1, dimlist, 1, "standard_name", "grid_edge_x_arc_type" ); else id_arcx = mpp_def_var(fid, "arcx", MPP_CHAR, 1, dimlist, 2, "standard_name", "grid_edge_x_arc_type", "north_pole", north_pole_arcx ); mpp_def_global_att(fid, "grid_version", GRID_VERSION); mpp_def_global_att(fid, "code_version", TAGNAME); mpp_def_global_att(fid, "history", history); mpp_end_def(fid); for(m=0; m<4; m++) { start[m] = 0; nwrite[m] = 0; } nwrite[0] = strlen(tilename); mpp_put_var_value_block(fid, id_tile, start, nwrite, tilename ); tmp = get_subregion(nxp, x+n*nxp*nyp, i*ni, (i+1)*ni, j*nj, (j+1)*nj); mpp_put_var_value(fid, id_x, tmp); free(tmp); tmp = get_subregion(nxp, y+n*nxp*nyp, i*ni, (i+1)*ni, j*nj, (j+1)*nj); mpp_put_var_value(fid, id_y, tmp); free(tmp); gdata = (double *)malloc(nx*nyp*sizeof(double)); mpp_global_field_double(domain, nxc, nyc+1, dx+n*nx*nyp, gdata); tmp = get_subregion( nx, gdata, i*ni, (i+1)*ni-1, j*nj, (j+1)*nj); mpp_put_var_value(fid, id_dx, tmp); free(tmp); free(gdata); gdata = (double *)malloc(nxp*ny*sizeof(double)); mpp_global_field_double(domain, nxc+1, nyc, dy+n*nxp*ny, gdata); tmp = get_subregion( nxp, gdata, i*ni, (i+1)*ni, j*nj, (j+1)*nj-1); mpp_put_var_value(fid, id_dy, tmp); free(tmp); free(gdata); gdata = (double *)malloc(nx*ny*sizeof(double)); mpp_global_field_double(domain, nxc, nyc, area+n*nx*ny, gdata); tmp = get_subregion( nx, gdata, i*ni, (i+1)*ni-1, j*nj, (j+1)*nj-1); mpp_put_var_value(fid, id_area, tmp); free(tmp); free(gdata); gdata = (double *)malloc(nxp*nyp*sizeof(double)); mpp_global_field_double(domain, nxc+1, nyc+1, angle_dx+n*nxp*nyp, gdata); tmp = get_subregion( nxp, gdata, i*ni, (i+1)*ni, j*nj, (j+1)*nj); mpp_put_var_value(fid, id_angle_dx, tmp); free(tmp); free(gdata); if(strcmp(conformal, "true") != 0) { gdata = (double *)malloc(nxp*nyp*sizeof(double)); mpp_global_field_double(domain, nxc+1, nyc+1, angle_dy+n*nxp*nyp, gdata); tmp = get_subregion( nxp, gdata, i*ni, (i+1)*ni, j*nj, (j+1)*nj); mpp_put_var_value(fid, id_angle_dy, tmp); free(tmp); free(gdata); } nwrite[0] = strlen(arcx); mpp_put_var_value_block(fid, id_arcx, start, nwrite, arcx ); mpp_close(fid); } } } } free(x); free(y); free(dx); free(dy); free(area); free(angle_dx); if(strcmp(conformal, "true") != 0) free(angle_dy); } int gs_fregrid(char *history, char *mosaic_in, char *mosaic_out, char *dir_in, char *dir_out, char **input_file, int nfiles, char **output_file, int nfiles_out, char *remap_file, char **scalar_name, int nscalar, char **u_name, int nvector, char **v_name, int nvector2, char *interp_method, char *test_case, double test_param, unsigned int opcode, int grid_type, unsigned int finer_step, int fill_missing, int nlon, int nlat, int check_conserve, int y_at_center, double lonbegin, double lonend, double latbegin, double latend, int kbegin, int kend, int lbegin, int lend) { int ntiles_in = 0; /* number of tiles in input mosaic */ int ntiles_out = 0; /* number of tiles in output mosaic */ int option_index, c, i, n, m, l; char txt[STRING]; Grid_config *grid_in = NULL; /* store input grid */ Grid_config *grid_out = NULL; /* store output grid */ Field_config *scalar_in = NULL; /* store input scalar data */ Field_config *scalar_out = NULL; /* store output scalar data */ Field_config *u_in = NULL; /* store input vector u-component */ Field_config *v_in = NULL; /* store input vector v-component */ Field_config *u_out = NULL; /* store input vector u-component */ Field_config *v_out = NULL; /* store input vector v-component */ File_config *file_in = NULL; /* store input file information */ File_config *file_out = NULL; /* store output file information */ File_config *file2_in = NULL; /* store input file information */ File_config *file2_out = NULL; /* store output file information */ Bound_config *bound_T = NULL; /* store halo update information for T-cell*/ Interp_config *interp = NULL; /* store remapping information */ int save_weight_only = 0; int fid; /* check the arguments */ if( !mosaic_in ) mpp_error("fregrid: input_mosaic is not specified"); if( !mosaic_out ) { if(nlon == 0 || nlat ==0 ) mpp_error("fregrid: when output_mosaic is not specified, nlon and nlat should be specified"); if(lonend <= lonbegin) mpp_error("fregrid: when output_mosaic is not specified, lonEnd should be larger than lonBegin"); if(latend <= latbegin) mpp_error("fregrid: when output_mosaic is not specified, latEnd should be larger than latBegin"); } else { if(nlon !=0 || nlat != 0) mpp_error("fregrid: when output_mosaic is specified, nlon and nlat should not be specified"); } if( nfiles == 0) { if(nvector > 0 || nscalar > 0 || nvector2 > 0) mpp_error("fregrid: when --input_file is not specified, --scalar_field, --u_field and --v_field should also not be specified"); if(!remap_file) mpp_error("fregrid: when --input_file is not specified, remap_file must be specified to save weight information"); save_weight_only = 1; if(mpp_pe()==mpp_root_pe())printf("NOTE: No input file specified in this run, no data file will be regridded " "and only weight information is calculated.\n"); } else if( nfiles == 1 || nfiles ==2) { if( nvector != nvector2 ) mpp_error("fregrid: number of fields specified in u_field must be the same as specified in v_field"); if( nscalar+nvector==0 ) mpp_error("fregrid: both scalar_field and vector_field are not specified"); /* when nvector =2 and nscalar=0, nfiles can be 2 otherwise nfiles must be 1 */ if( nscalar && nfiles != 1 ) mpp_error("fregrid: when scalar_field is specified, number of files must be 1"); if( nfiles_out == 0 ) { for(i=0; i 0) { opcode |= VECTOR; if(grid_type == AGRID) opcode |= AGRID; else if(grid_type == BGRID) opcode |= BGRID; } /* get the mosaic information of input and output mosaic*/ fid = mpp_open(mosaic_in, MPP_READ); ntiles_in = mpp_get_dimlen(fid, "ntiles"); mpp_close(fid); if(mosaic_out) { fid = mpp_open(mosaic_out, MPP_READ); ntiles_out = mpp_get_dimlen(fid, "ntiles"); mpp_close(fid); } else ntiles_out = 1; if(!strcmp(interp_method, "conserve_order1") ) { if(mpp_pe() == mpp_root_pe())printf("****fregrid: first order conservative scheme will be used for regridding.\n"); opcode |= CONSERVE_ORDER1; } else if(!strcmp(interp_method, "conserve_order2") ) { if(mpp_pe() == mpp_root_pe())printf("****fregrid: second order conservative scheme will be used for regridding.\n"); opcode |= CONSERVE_ORDER2; } else if(!strcmp(interp_method, "bilinear") ) { if(mpp_pe() == mpp_root_pe())printf("****fregrid: bilinear remapping scheme will be used for regridding.\n"); opcode |= BILINEAR; } else mpp_error("fregrid: interp_method must be 'conserve_order1', 'conserve_order2' or 'bilinear'"); if(test_case) { if(nfiles != 1) mpp_error("fregrid: when test_case is specified, nfiles should be 1"); sprintf(output_file[0], "%s.%s.output", test_case, interp_method); } if(check_conserve) opcode |= CHECK_CONSERVE; if( opcode & BILINEAR ) { int ncontact; ncontact = read_mosaic_ncontacts(mosaic_in); if( nlon == 0 || nlat == 0) mpp_error("fregrid: when interp_method is bilinear, nlon and nlat should be specified"); if(ntiles_in != 6) mpp_error("fregrid: when interp_method is bilinear, the input mosaic should be 6 tile cubic grid"); if(ncontact !=12) mpp_error("fregrid: when interp_method is bilinear, the input mosaic should be 12 contact cubic grid"); if(mpp_npes() > 1) mpp_error("fregrid: parallel is not implemented for bilinear remapping"); } else y_at_center = 1; /* memory allocation for data structure */ grid_in = (Grid_config *)malloc(ntiles_in *sizeof(Grid_config)); grid_out = (Grid_config *)malloc(ntiles_out*sizeof(Grid_config)); bound_T = (Bound_config *)malloc(ntiles_in *sizeof(Bound_config)); interp = (Interp_config *)malloc(ntiles_out*sizeof(Interp_config)); get_input_grid( ntiles_in, grid_in, bound_T, mosaic_in, opcode ); if(mosaic_out) get_output_grid_from_mosaic( ntiles_out, grid_out, mosaic_out, opcode ); else get_output_grid_by_size(ntiles_out, grid_out, lonbegin, lonend, latbegin, latend, nlon, nlat, finer_step, y_at_center, opcode); if(remap_file) set_remap_file(ntiles_out, mosaic_out, remap_file, interp, &opcode, save_weight_only); /* preparing for the interpolation, if remapping information exist, read it from remap_file, otherwise create the remapping information and write it to remap_file */ if( opcode & BILINEAR ) /* bilinear interpolation from cubic to lalon */ setup_bilinear_interp(ntiles_in, grid_in, ntiles_out, grid_out, interp, opcode ); else setup_conserve_interp(ntiles_in, grid_in, ntiles_out, grid_out, interp, opcode); if(save_weight_only) { if(mpp_pe() == mpp_root_pe() ) { printf("NOTE: Successfully running fregrid and the following files which store weight information are generated.\n"); for(n=0; n 0) { scalar_in = (Field_config *)malloc(ntiles_in *sizeof(Field_config)); scalar_out = (Field_config *)malloc(ntiles_out *sizeof(Field_config)); } if(nvector > 0) { u_in = (Field_config *)malloc(ntiles_in *sizeof(Field_config)); u_out = (Field_config *)malloc(ntiles_out *sizeof(Field_config)); v_in = (Field_config *)malloc(ntiles_in *sizeof(Field_config)); v_out = (Field_config *)malloc(ntiles_out *sizeof(Field_config)); } set_mosaic_data_file(ntiles_in, mosaic_in, dir_in, file_in, input_file[0]); set_mosaic_data_file(ntiles_out, mosaic_out, dir_out, file_out, output_file[0]); if(nfiles == 2) { set_mosaic_data_file(ntiles_in, mosaic_in, dir_in, file2_in, input_file[1]); set_mosaic_data_file(ntiles_out, mosaic_out, dir_out, file2_out, output_file[1]); } for(n=0; n 0) get_field_missing(ntiles_in, scalar_in); if(nvector > 0) { get_field_missing(ntiles_in, u_in); get_field_missing(ntiles_in, v_in); } /* set time step to 1, only test scalar field now, nz need to be 1 */ if(test_case) { if(nscalar != 1 || nvector != 0) mpp_error("fregrid: when test_case is specified, nscalar must be 1 and nvector must be 0"); if(scalar_in->var->nz != 1) mpp_error("fregrid: when test_case is specified, number of vertical level must be 1"); file_in->nt = 1; file_out->nt = 1; } /* Then doing the regridding */ for(m=0; mnt; m++) { int memsize, level; write_output_time(ntiles_out, file_out, m); if(nfiles > 1) write_output_time(ntiles_out, file2_out, m); /* first interp scalar variable */ for(l=0; lvar[l].has_taxis && m>0) continue; level = m + scalar_in->var[l].lstart; if(test_case) get_test_input_data(test_case, test_param, ntiles_in, scalar_in, grid_in, bound_T, opcode); else get_input_data(ntiles_in, scalar_in, grid_in, bound_T, l, level, opcode); allocate_field_data(ntiles_out, scalar_out, grid_out, l); if( opcode & BILINEAR ) do_scalar_bilinear_interp(interp, l, ntiles_in, grid_in, grid_out, scalar_in, scalar_out, finer_step, fill_missing); else do_scalar_conserve_interp(interp, l, ntiles_in, grid_in, ntiles_out, grid_out, scalar_in, scalar_out, opcode); write_field_data(ntiles_out, scalar_out, grid_out, l, m); if(opcode & CONSERVE_ORDER2) { for(n=0; n0) continue; level = m + u_in->var[l].lstart; get_input_data(ntiles_in, u_in, grid_in, bound_T, l, level, opcode); get_input_data(ntiles_in, v_in, grid_in, bound_T, l, level, opcode); allocate_field_data(ntiles_out, u_out, grid_out, l); allocate_field_data(ntiles_out, v_out, grid_out, l); if( opcode & BILINEAR ) do_vector_bilinear_interp(interp, l, ntiles_in, grid_in, ntiles_out, grid_out, u_in, v_in, u_out, v_out, finer_step, fill_missing); else do_vector_conserve_interp(interp, l, ntiles_in, grid_in, ntiles_out, grid_out, u_in, v_in, u_out, v_out, opcode); write_field_data(ntiles_out, u_out, grid_out, l, m); write_field_data(ntiles_out, v_out, grid_out, l, m); for(n=0; n 1 ) { mpp_close(file2_out[n].fid); printf("****%s\n", file2_out[n].name); } } } return 0; } int gs_make_topog(char *history, char *mosaic_file, char *topog_type, int x_refine, int y_refine, double basin_depth, char *topog_file, double bottom_depth, double min_depth, double scale_factor, int num_filter_pass, double gauss_amp, double gauss_scale, double slope_x, double slope_y, double bowl_south, double bowl_north, double bowl_west, double bowl_east, int fill_first_row, int filter_topog, int round_shallow, int fill_shallow, int deepen_shallow, int smooth_topo_allow_deepening, char *output_file) { char *topog_field = NULL; int option_index, i, c; if(x_refine != 2 || y_refine != 2 ) mpp_error("Error from make_topog: x_refine and y_refine should be 2, contact developer"); if(mpp_pe() == mpp_root_pe()) printf("NOTE from make_topog ==> x_refine is %d, y_refine is %d\n", x_refine, y_refine); if (strcmp(topog_type,"rectangular_basin") == 0) { if(mpp_pe() == mpp_root_pe()) printf("NOTE from make_topog ==> the basin depth is %f\n",basin_depth); } else if (strcmp(topog_type,"gaussian") == 0) { if(mpp_pe() == mpp_root_pe()){ printf("NOTE from make_topog ==> bottom_depth is: %f\n", bottom_depth ); printf("NOTE from make_topog ==> min_depth is: %f\n", min_depth ); printf("NOTE from make_topog ==> gauss_amp is: %f\n", gauss_amp ); printf("NOTE from make_topog ==> gauss_scale is: %f\n", gauss_scale ); printf("NOTE from make_topog ==> slope_x is: %f\n", slope_x ); printf("NOTE from make_topog ==> slope_y is: %f\n", slope_y ); } } else if(strcmp(topog_type,"bowl") == 0) { if(mpp_pe() == mpp_root_pe()){ printf("NOTE from make_topog ==> bottom_depth is: %f\n",bottom_depth); printf("NOTE from make_topog ==> min_depth is: %f\n",min_depth); printf("NOTE from make_topog ==> bowl_south is: %f\n",bowl_south); printf("NOTE from make_topog ==> bowl_north is: %f\n",bowl_north); printf("NOTE from make_topog ==> bowl_west is: %f\n",bowl_west); printf("NOTE from make_topog ==> bowl_east is: %f\n",bowl_east); } } else if(strcmp(topog_type,"idealized") == 0) { if(mpp_pe() == mpp_root_pe()){ printf("NOTE from make_topog ==> bottom_depth is: %f\n",bottom_depth); printf("NOTE from make_topog ==> min_depth is: %f\n",min_depth); } } else if(strcmp(topog_type,"realistic") == 0) { if(!topog_file || !topog_field) mpp_error("Error from make_topog: when topog_type is realistic, topog_file and topog_field must be specified."); if(mpp_pe() == mpp_root_pe()){ printf("NOTE from make_topog ==> bottom_depth is: %f\n",bottom_depth); printf("NOTE from make_topog ==> min_depth is: %f\n",min_depth); printf("NOTE from make_topog ==> topog_file is: %s\n", topog_file); printf("NOTE from make_topog ==> topog_field is: %s\n", topog_field); printf("NOTE from make_topog ==> scale_factor is: %f\n", scale_factor); printf("NOTE from make_topog ==> num_filter_pass is: %d\n", num_filter_pass); if(fill_first_row) printf("NOTE from make_topog ==>make first row of ocean model all land points.\n"); if(filter_topog) printf("NOTE from make_topog ==>will apply filter to topography.\n"); if(round_shallow) printf("NOTE from make_topog ==>Make cells land if depth is less than 1/2 " "mimumim depth, otherwise make ocean.\n"); if(fill_shallow) printf("NOTE from make_topog ==>Make cells less than minimum depth land.\n"); if(deepen_shallow) printf("NOTE from make_topog ==>Make cells less than minimum depth equal to minimum depth.\n"); if(smooth_topo_allow_deepening) printf("NOTE from make_topog ==>allow filter to deepen cells.\n"); } } else { mpp_error("make_topog: topog_type should be rectangular_basin, gaussian, bowl, idealized or realistic"); } if(mpp_pe() == mpp_root_pe()) { printf("**************************************************\n"); printf("Begin to generate topography \n"); } { #define STRING 255 int m_fid, g_fid, vid; int ntiles, fid, dim_ntiles, n, dims[2]; size_t start[4], nread[4], nwrite[4]; int *nx, *ny, *nxp, *nyp; int *id_depth; double *depth, *x, *y; char **tile_files; char history[512], dimx_name[128], dimy_name[128], depth_name[128]; char gridfile[256], griddir[256]; /* history will be write out as global attribute in output file to specify the command line arguments */ /* grid should be located in the same directory of mosaic file */ get_file_path(mosaic_file, griddir); /* get mosaic dimension */ m_fid = mpp_open(mosaic_file, MPP_READ); ntiles = mpp_get_dimlen( m_fid, "ntiles"); tile_files = (char **)malloc(ntiles*sizeof(double *)); id_depth = (int *)malloc(ntiles*sizeof(int)); /* loop through each tile to get tile information and set up meta data for output file */ fid = mpp_open(output_file, MPP_WRITE); mpp_def_global_att(fid, "grid_version", GRID_VERSION); mpp_def_global_att(fid, "code_version", TAGNAME); mpp_def_global_att(fid, "history", history); dim_ntiles = mpp_def_dim(fid, "ntiles", ntiles); nx = (int *)malloc(ntiles*sizeof(int)); ny = (int *)malloc(ntiles*sizeof(int)); nxp = (int *)malloc(ntiles*sizeof(int)); nyp = (int *)malloc(ntiles*sizeof(int)); for( n = 0; n < ntiles; n++ ) { tile_files[n] = (char *)malloc(STRING*sizeof(double)); start[0] = n; start[1] = 0; nread[0] = 1; nread[1] = STRING; vid = mpp_get_varid(m_fid, "gridfiles"); mpp_get_var_value_block(m_fid, vid, start, nread, gridfile); sprintf(tile_files[n], "%s/%s", griddir, gridfile); g_fid = mpp_open(tile_files[n], MPP_READ); nx[n] = mpp_get_dimlen(g_fid, "nx"); ny[n] = mpp_get_dimlen(g_fid, "ny"); if( nx[n]%x_refine != 0 ) mpp_error("make_topog: supergrid x-size can not be divided by x_refine"); if( ny[n]%y_refine != 0 ) mpp_error("make_topog: supergrid y-size can not be divided by y_refine"); nx[n] /= x_refine; ny[n] /= y_refine; nxp[n] = nx[n] + 1; nyp[n] = ny[n] + 1; if(ntiles == 1) { strcpy(dimx_name, "nx"); strcpy(dimy_name, "ny"); strcpy(depth_name, "depth"); } else { sprintf(dimx_name, "nx_tile%d", n+1); sprintf(dimy_name, "ny_tile%d", n+1); sprintf(depth_name, "depth_tile%d", n+1); } dims[1] = mpp_def_dim(fid, dimx_name, nx[n]); dims[0] = mpp_def_dim(fid, dimy_name, ny[n]); id_depth[n] = mpp_def_var(fid, depth_name, NC_DOUBLE, 2, dims, 2, "standard_name", "topographic depth at T-cell centers", "units", "meters"); mpp_close(g_fid); } mpp_close(m_fid); mpp_end_def(fid); /* Generate topography and write out to the output_file */ for(n=0; nMAXCONTACT) mpp_error("make_solo_mosaic: number of contacts is more than MAXCONTACT"); for(l=0; l0) dim_ncontact = mpp_def_dim(fid, "ncontact", ncontact); dim_string = mpp_def_dim(fid, "string", STRING); /* define variable */ id_mosaic = mpp_def_var(fid, "mosaic", MPP_CHAR, 1, &dim_string, 4, "standard_name", "grid_mosaic_spec", "children", "gridtiles", "contact_regions", "contacts", "grid_descriptor", grid_descriptor); dim[0] = dim_ntiles; dim[1] = dim_string; id_griddir = mpp_def_var(fid, "gridlocation", MPP_CHAR, 1, &dim[1], 1, "standard_name", "grid_file_location"); id_gridfiles = mpp_def_var(fid, "gridfiles", MPP_CHAR, 2, dim, 0); id_gridtiles = mpp_def_var(fid, "gridtiles", MPP_CHAR, 2, dim, 0); if(ncontact>0) { dim[0] = dim_ncontact; dim[1] = dim_string; id_contacts = mpp_def_var(fid, "contacts", MPP_CHAR, 2, dim, 5, "standard_name", "grid_contact_spec", "contact_type", "boundary", "alignment", "true", "contact_index", "contact_index", "orientation", "orient"); id_contact_index = mpp_def_var(fid, "contact_index", MPP_CHAR, 2, dim, 1, "standard_name", "starting_ending_point_index_of_contact"); } mpp_def_global_att(fid, "grid_version", GRID_VERSION); mpp_def_global_att(fid, "code_version", TAGNAME); mpp_def_global_att(fid, "history", history); mpp_end_def(fid); /* write out data */ for(i=0; i<4; i++) { start[i] = 0; nwrite[i] = 1; } nwrite[0] = strlen(mosaic_name); mpp_put_var_value_block(fid, id_mosaic, start, nwrite, mosaic_name); nwrite[0] = strlen(dir); mpp_put_var_value_block(fid, id_griddir, start, nwrite, dir); nwrite[0] = 1; for(n=0; nnx = nx; river_data->ny = ny; river_data->xt = (double *)malloc(nx*sizeof(double)); river_data->yt = (double *)malloc(ny*sizeof(double)); river_data->xb = (double *)malloc(nxp*sizeof(double)); river_data->yb = (double *)malloc(nyp*sizeof(double)); river_data->xb_r = (double *)malloc(nxp*sizeof(double)); river_data->yb_r = (double *)malloc(nyp*sizeof(double)); river_data->subA = (double *)malloc(nx*ny*sizeof(double)); mpp_get_var_att(fid, vid, "missing_value", &(river_data->subA_missing) ); vid = mpp_get_varid(fid, tocell_name); mpp_get_var_att(fid, vid, "missing_value", &dbl_missing ); river_data->tocell_missing = dbl_missing; vid = mpp_get_varid(fid, travel_name); mpp_get_var_att(fid, vid, "missing_value", &dbl_missing ); river_data->travel_missing = dbl_missing; vid = mpp_get_varid(fid, basin_name); mpp_get_var_att(fid, vid, "missing_value", &dbl_missing ); river_data->basin_missing = dbl_missing; vid = mpp_get_varid(fid, cellarea_name); mpp_get_var_att(fid, vid, "missing_value", &(river_data->cellarea_missing) ); vid = mpp_get_varid(fid, celllength_name); mpp_get_var_att(fid, vid, "missing_value", &(river_data->celllength_missing) ); vid = mpp_get_varid(fid, subA_name); mpp_get_var_value(fid, vid, river_data->subA); vid = mpp_get_varid(fid, xaxis_name); mpp_get_var_value(fid, vid, river_data->xt); vid = mpp_get_varid(fid, yaxis_name); mpp_get_var_value(fid, vid, river_data->yt); mpp_close(fid); for(i=1; ixb[i] = 0.5*(river_data->xt[i-1]+river_data->xt[i]); for(j=1; jyb[j] = 0.5*(river_data->yt[j-1]+river_data->yt[j]); river_data->xb[0] = 2*river_data->xt[0] - river_data->xb[1]; river_data->yb[0] = 2*river_data->yt[0] - river_data->yb[1]; river_data->xb[nx] = 2*river_data->xt[nx-1] - river_data->xb[nx-1]; river_data->yb[ny] = 2*river_data->yt[ny-1] - river_data->yb[ny-1]; /* make sure the xb is in the range [0,360] and yb is in the range [-90,90] */ if(fabs(river_data->xb[0]) > EPSLN) mpp_error("river_regrid: The starting longitude of grid bound is not 0"); if(fabs(river_data->xb[nx]-360) > EPSLN) mpp_error("river_regrid: The ending longitude of grid bound is not 360"); if(fabs(river_data->yb[0]+90) > EPSLN) mpp_error("river_regrid: The starting latitude of grid bound is not -90"); if(fabs(river_data->yb[ny]-90) > EPSLN) mpp_error("river_regrid: The ending latitude of grid bound is not 90"); river_data->xb[0] = 0.; river_data->xb[nx] = 360.; river_data->yb[0] = -90.; river_data->yb[ny] = 90.; for(j=0; jyb_r[j] = river_data->yb[j]*D2R; for(i=0; ixb_r[i] = river_data->xb[i]*D2R; } /*---------------------------------------------------------------------- read the grid of detination mosaic. void get_mosaic_grid(char *file) where file is the coupler mosaic file. --------------------------------------------------------------------*/ void get_mosaic_grid(const char *coupler_mosaic, const char *land_mosaic, int ntiles, river_type *river_data, unsigned int *opcode) { int n_xgrid_files, nx, ny, nxp, nyp, ni, nj, nip, njp; int n, m, i, j, ii, jj, nxp2, nyp2; size_t start[4], nread[4]; char dir[STRING], gridfile[STRING], tilefile[STRING]; char tilename[STRING], land_mosaic_name[STRING]; char **xgrid_file; double *x, *y; double **pxt, **pyt; char *pfile; int m_fid, m_vid, g_fid, g_vid; /* coupler_mosaic, land_mosaic, and exchange grid file should be located in the same directory */ get_file_path(coupler_mosaic, dir); m_fid = mpp_open(coupler_mosaic, MPP_READ); m_vid = mpp_get_varid(m_fid, "lnd_mosaic"); mpp_get_var_value(m_fid, m_vid, land_mosaic_name); /* get the exchange grid file name */ n_xgrid_files = mpp_get_dimlen(m_fid, "nfile_aXl"); xgrid_file = (char **)malloc(n_xgrid_files*sizeof(char *)); m_vid = mpp_get_varid(m_fid, "aXl_file"); for(n=0; n 1 || river_data[n].landfrac[m] < 0) mpp_error("river_regrid: land_frac should be between 0 or 1"); } free(area); } /* n = 0, ntiles */ mpp_close(m_fid); /* currently we are assuming all the times have the same grid size */ for(n=1; n0) { read_mosaic_contact(land_mosaic, tile1, tile2, istart1, iend1, jstart1, jend1, istart2, iend2, jstart2, jend2); if(ncontact <= 2) { /* x-cyclic of y-cyclic */ if(ntiles !=1) mpp_error("river_regrid: number of tiles must be 1 for single contact"); for(n=0; nnx; ny = river_out->ny; nxp2 = nx + 2; nyp2 = ny + 2; subA_missing = river_in->subA_missing; tocell_missing = river_in->tocell_missing; travel_missing = river_in->travel_missing; basin_missing = river_in->basin_missing; cellarea_missing = river_in->cellarea_missing; celllength_missing = river_in->celllength_missing; for(n=0; nnx; ny_in = river_in->ny; xb_in = river_in->xb_r; yb_in = river_in->yb_r; missing = river_out->subA_missing; for(n=0; n=ur_y) && (yv1[1]>=ur_y) && (yv1[2]>=ur_y) && (yv1[3]>=ur_y) ) continue; for(ii=0; iisubA[jj*nx_in+ii] == missing) continue; ll_x = xb_in[ii]; ur_x = xb_in[ii+1]; /* adjust xv1 to make sure it is in the range as ll_x and ur_x */ adjust_lon(xv1, 0.5*(ll_x + ur_x) ); if ( (xv1[0]<=ll_x) && (xv1[1]<=ll_x) && (xv1[2]<=ll_x) && (xv1[3]<=ll_x) ) continue; if ( (xv1[0]>=ur_x) && (xv1[1]>=ur_x) && (xv1[2]>=ur_x) && (xv1[3]>=ur_x) ) continue; if ( (n_out = clip ( xv1, yv1, 4, ll_x, ll_y, ur_x, ur_y, xv2, yv2 )) > 0 ) { double xarea; xarea = poly_area (xv2, yv2, n_out ); if( xarea/river_out[n].area[j*nx_out+i] < MIN_AREA_RATIO ) continue; if(river_in->subA[jj*nx_in+ii] > river_out[n].subA[j1*nxp2+i1]) river_out[n].subA[j1*nxp2+i1] = river_in->subA[jj*nx_in+ii]; } } } } } /* fill the halo of subA */ psubA = (double **)malloc(ntiles*sizeof(double *)); for(n=0; nnx, river_out->ny, opcode); free(psubA); };/* calc_max_subA */ /*------------------------------------------------------------------------------ void update_halo_double(int ntiles, double *data, int nx, int ny, unsigned int opcode) We assume all the tiles have the same size. -----------------------------------------------------------------------------*/ void update_halo_double(int ntiles, double **data, int nx, int ny, unsigned int opcode) { int nxp1, nyp1, nxp2, i, j, n; int te, tw, ts, tn; nxp1 = nx + 1; nyp1 = ny + 1; nxp2 = nx + 2; if(opcode & X_CYCLIC) { for(j=1; j M_PI) dx = dx - 2*M_PI; x_sum += (x[i] = x[i-1] + dx); } dx = (x_sum/npts)-tlon; if (dx < -M_PI) for (i=0;i M_PI) for (i=0;inx; ny = river_data->ny; nxp2 = nx + 2; subA_missing = river_data->subA_missing; ptocell = (int **)malloc(ntiles*sizeof(int *)); pdir = (int **)malloc(ntiles*sizeof(int *)); for(n=0; n suba_cutoff ) { /* tocell will be land cell with larger subA value, instead of neighboring ocean cell. */ for(joff=0; joff subA_me) { if(tval == -999 ) { if(subA > tval) { iget = ioff; jget = joff; tval = subA; } } else if(tval < LARGE_VALUE) { if(subA > tval && subA < LARGE_VALUE) { iget = ioff; jget = joff; tval = subA; } } else { iget = ioff; jget = joff; tval = subA; } } } } } else { for(joff=0; joff= tval) { iget = ioff; jget = joff; tval = subA; } } } } } if(iget >= 0 && jget >= 0) { if(tval > river_data[n].subA[j*nxp2+i]) { ptocell[n][j*nxp2+i] = out_flow[jget*ncells+iget]; pdir [n][j*nxp2+i] = out_dir [jget*ncells+iget]; if(pow(2,pdir[n][j*nxp2+i]) != ptocell[n][j*nxp2+i] ) mpp_error("river_regrid: pow(2,dir) should equal to tocell"); } else ptocell[n][j*nxp2+i] = 0; } else ptocell[n][j*nxp2+i] = 0; } } } update_halo_int(ntiles, ptocell, nx, ny, opcode); update_halo_int(ntiles, pdir, nx, ny, opcode); free(ptocell); free(pdir); };/* calc_tocell */ /*------------------------------------------------------------------------------ void calc_river_data() calculate travel and basin according to tocell, as well as celllength, cellarea. For each basin, one and only one river point will have tocell = 0, this point will have travel = 1. ----------------------------------------------------------------------------*/ void calc_river_data(int ntiles, river_type* river_data, unsigned int opcode ) { int nx, ny, nxp2, nyp2, n, i, j, ii, jj, nxp, nyp, im1, jm1, ioff, joff; int cur_basin, cur_travel; int not_done, dir, i_dest, j_dest; int basin_missing, travel_missing, tocell_missing; int maxtravel, maxbasin, travelnow; double subA_missing; int **pbasin, **ptravel, **pdir; double **psubA; double xv[4], yv[4]; const int di[] = {1,1,0,-1,-1,-1,0,1}; const int dj[] = {0,-1,-1,-1,0,1,1,1}; nx = river_data->nx; ny = river_data->ny; nxp = nx+1; nyp = ny+1; nxp2 = nx+2; nyp2 = ny+2; basin_missing = river_data->basin_missing; subA_missing = river_data->subA_missing; tocell_missing = river_data->tocell_missing; cur_basin = 0; /* set up pointer */ pbasin = (int **)malloc(ntiles*sizeof(int *)); ptravel = (int **)malloc(ntiles*sizeof(int *)); pdir = (int **)malloc(ntiles*sizeof(int *)); psubA = (double **)malloc(ntiles*sizeof(double *)); for(n=0; n 0) { ii = i+1; jj = j+1; xv[0] = river_data[n].xb_r[j*nxp+i]; xv[1] = river_data[n].xb_r[j*nxp+ii]; xv[2] = river_data[n].xb_r[jj*nxp+ii]; xv[3] = river_data[n].xb_r[jj*nxp+i]; yv[0] = river_data[n].yb_r[j*nxp+i]; yv[1] = river_data[n].yb_r[j*nxp+ii]; yv[2] = river_data[n].yb_r[jj*nxp+ii]; yv[3] = river_data[n].yb_r[jj*nxp+i]; river_data[n].cellarea[j*nx+i] = river_data[n].area[j*nx+i]*river_data[n].landfrac[j*nx+i]; dir = pdir[n][jj*nxp2+ii]; if( dir >= 0) { i_dest = ii + di[dir]; j_dest = jj + dj[dir]; river_data[n].celllength[j*nx+i] = distance(river_data[n].xt[jj*nxp2+ii], river_data[n].yt[jj*nxp2+ii], river_data[n].xt[j_dest*nxp2+i_dest], river_data[n].yt[j_dest*nxp2+i_dest] ); } else river_data[n].celllength[j*nx+i] = 0; } } } /* define the basinid and travel for the coast point */ for(n=0; n 0){ dir = pdir[n][j*nxp2+i]; i_dest = i + di[dir]; j_dest = j + dj[dir]; if(river_data[n].tocell[j_dest*nxp2+i_dest] == tocell_missing) { pbasin[n] [j*nxp2+i] = ++cur_basin; ptravel[n][j*nxp2+i] = 1; river_data[n].last_point[(j-1)*nx+i-1] = 1; } } } } update_halo_int(ntiles, pbasin, nx, ny, opcode); update_halo_int(ntiles, ptravel, nx, ny, opcode); /* then define the travel and basin for all other points */ cur_travel = 0; not_done = 1; while(not_done) { not_done = 0; for(n=0; n= 0 && pbasin[n][j*nxp2+i] == basin_missing ) { i_dest = i + di[dir]; j_dest = j + dj[dir]; if( ptravel[n][j_dest*nxp2+i_dest] == cur_travel ) { not_done = 1; /* still have points need to be updateed */ if(pbasin[n][j_dest*nxp2+i_dest] == basin_missing) { mpp_error("river_grid: the tocell should have valid basin value"); } pbasin[n][j*nxp2+i] = pbasin[n][j_dest*nxp2+i_dest]; ptravel[n][j*nxp2+i] = cur_travel+1; } } } } cur_travel++; update_halo_int(ntiles, pbasin, nx, ny, opcode); update_halo_int(ntiles, ptravel, nx, ny, opcode); } /* figure out maximum travel and maximum basin*/ maxtravel = -1; maxbasin = -1; basin_missing = river_data->basin_missing; for(n=0; n=0; travelnow--) { for(n=0; n=0 ) { if( opcode & CUBIC_GRID ) { if(ii == 0) { /* west */ if(n%2==0) dir = (dir+2)%8; /* tile 1 3 5 */ } else if(ii == nxp) { /*east */ if(n%2==1) dir = (dir+2)%8; /* tile 2 4 6 */ } else if(jj == 0) { /* south */ if(n%2==1) dir = (dir+6)%8; /* tile 2 4 6 */ } else if(jj == nyp) { /*north */ if(n%2==0) dir = (dir+6)%8; /* tile 1 3 5 */ } } i_dest = ii + di[dir]; j_dest = jj + dj[dir]; if(i_dest == i && j_dest == j) psubA[n][j*nxp2+i] += psubA[n][jj*nxp2+ii]; } } } } update_halo_double(ntiles, psubA, nx, ny, opcode); } free(pbasin); free(ptravel); free(psubA); free(pdir); }; /* calc_travel */ /*------------------------------------------------------------------------------ void check_river_data( ) check to make sure all the river points have been assigned travel and basin value and all the ocean points will have missing value. ----------------------------------------------------------------------------*/ void check_river_data(int ntiles, river_type *river_data ) { int maxtravel = -1; int maxbasin = -1; int nx, ny, nxp2, nyp2; int n, i, j, im1, jm1, ioff, joff, ii, jj; int tocell, travel, basin; int tocell_missing, travel_missing, basin_missing; int ncoast_full_land, nsink; double subA, subA_missing; int *ncoast; /* print out maximum travel and number of rivers */ maxtravel = -1; maxbasin = -1; nx = river_data->nx; ny = river_data->ny; nxp2 = nx + 2; nyp2 = ny + 2; subA_missing = river_data->subA_missing; tocell_missing = river_data->tocell_missing; basin_missing = river_data->basin_missing; travel_missing = river_data->travel_missing; for(n=0; n NOTE from river_regrid: maximum travel is %d and maximum basin is %d.\n", maxtravel, maxbasin); ncoast_full_land = 0; nsink = 0; for(n=0; n 0) printf("Warning from river_regrid: there are %d coast points is a full land cell\n", ncoast_full_land); else printf("NOTE from river_regrid: there are no coast points is a full land cell\n"); if(nsink > 0) printf("Warning from river_regrid: there are %d sink points is a full land cell\n", nsink); else printf("NOTE from river_regrid: there are no sink points is a full land cell\n"); /* check river travel to make sure there is one and only one point with travel = 0. */ ncoast = (int *)malloc(maxbasin*sizeof(int)); for(n=0; n1) { printf("river with basin = %d has more than one point with travel = 0.\n", n+1); mpp_error("river_regrid: some river has more than one point with travel = 0"); } } free(ncoast); }; /* check_river_data */ /*------------------------------------------------------------------------------ void sort_basin(int ntiles, river_type* river_data); sorting the basin according to subA The river with larger subA at coast point will get smaller basinid. The basinid will be reorganized. -----------------------------------------------------------------------------*/ void sort_basin(int ntiles, river_type* river_data) { int nx, ny, nxp2, nyp2, i, j, n; int maxbasin, basin_missing, basin; int *rank, *indx; double subA_missing, *maxsuba; river_type *river_tmp; nx = river_data->nx; ny = river_data->ny; nxp2 = nx + 2; nyp2 = ny + 2; river_tmp = (river_type *)malloc(ntiles*sizeof(river_type)); /* calculate maximum basin */ maxbasin = -1; basin_missing = river_data->basin_missing; for(n=0; nsubA_missing; maxsuba = (double *)malloc(maxbasin*sizeof(double)); indx = (int *)malloc(maxbasin*sizeof(int )); rank = (int *)malloc(maxbasin*sizeof(int )); for(n=0; n maxbasin || basin < 1) mpp_error("river_regrid: basin should be between 1 and maxbasin"); maxsuba[basin-1] = river_data[n].subA[j*nxp2+i]; } } /* make sure maxsuba is assigned properly */ for(n=0; n maxbasin || basin < 1) mpp_error("river_regrid: basin should be a positive integer no larger than maxbasin"); river_data[n].basin[j*nxp2+i] = maxbasin - rank[basin-1]; } } /* release the memory */ for(n=0; n1) sprintf(river_data[n].filename, "%s.tile%d.nc", output_file, n+1); else sprintf(river_data[n].filename, "%s.nc", output_file); nx = river_data[n].nx; ny = river_data[n].ny; nxp2 = nx + 2; fid = mpp_open(river_data[n].filename, MPP_WRITE); mpp_def_global_att(fid, "version", version); mpp_def_global_att(fid, "code_version", TAGNAME); mpp_def_global_att(fid, "history", history); dimid[1] = mpp_def_dim(fid, gridx_name, nx); dimid[0] = mpp_def_dim(fid, gridy_name, ny); id_xaxis = mpp_def_var(fid, gridx_name, MPP_DOUBLE, 1, &(dimid[1]), 0); id_yaxis = mpp_def_var(fid, gridy_name, MPP_DOUBLE, 1, dimid, 0); mpp_copy_var_att(fid_in, id_xaxis_in, fid, id_xaxis); mpp_copy_var_att(fid_in, id_yaxis_in, fid, id_yaxis); id_subA = mpp_def_var(fid, subA_name, MPP_DOUBLE, 2, dimid , 0); mpp_copy_var_att(fid_in, id_subA_in, fid, id_subA); id_tocell = mpp_def_var(fid, tocell_name, MPP_INT, 2, dimid, 0); mpp_copy_var_att(fid_in, id_tocell_in, fid, id_tocell); id_travel = mpp_def_var(fid, travel_name, MPP_INT, 2, dimid , 0); mpp_copy_var_att(fid_in, id_travel_in, fid, id_travel); id_basin = mpp_def_var(fid, basin_name, MPP_INT, 2, dimid , 0); mpp_copy_var_att(fid_in, id_basin_in, fid, id_basin); id_cellarea = mpp_def_var(fid, cellarea_name, MPP_DOUBLE, 2, dimid , 0); mpp_copy_var_att(fid_in, id_cellarea_in, fid, id_cellarea); id_celllength = mpp_def_var(fid, celllength_name, MPP_DOUBLE, 2, dimid , 0); mpp_copy_var_att(fid_in, id_celllength_in, fid, id_celllength); id_landfrac = mpp_def_var(fid, landfrac_name, MPP_DOUBLE, 2, dimid, 2, "long_name", "land fraction", "units", "none"); id_x = mpp_def_var(fid, x_name, MPP_DOUBLE, 2, dimid, 2, "long_name", "Geographic longitude", "units", "degree_east"); id_y = mpp_def_var(fid, y_name, MPP_DOUBLE, 2, dimid, 2, "long_name", "Geographic latitude", "units", "degree_north"); mpp_end_def(fid); xt = (double *)malloc(nx*sizeof(double)); yt = (double *)malloc(ny*sizeof(double)); /* for lat-lon grid, actual lon, lat will be written out, for cubic grid, index will be written out */ if(ntiles == 1) { for(i=0; i 180.) dx = dx - 360.; if(dx < -180.) dx = dx + 360.; if(lon1 == lon2) dist = fabs(lat2 - lat1) * D2R; else if(lat1 == lat2) dist = fabs(dx) * D2R * cos(lat1*D2R); else { /* diagonal distance */ s1 = fabs(dx) * D2R * cos(lat1*D2R); s2 = fabs(lat2 - lat1) * D2R; dist = sqrt(s1*s1+s2*s2); } dist *= RADIUS; return dist; } /*------------------------------------------------------------------------------ void qsort_index(double array[], int start, int end, int rank[]) sort array in increasing order and get the rank of each member in the original array the array size of array and rank will be size. ----------------------------------------------------------------------------*/ #define swap(a,b,t) ((t)=(a),(a)=(b),(b)=(t)) void qsort_index(double array[], int start, int end, int rank[]) { double pivot; int tmp_int; double tmp_double; if (end > start) { int l = start + 1; int r = end; int pivot_index = (start+(end-start)/2); swap(array[start], array[pivot_index], tmp_double); /*** choose arbitrary pivot ***/ swap(rank[start], rank[pivot_index], tmp_int); pivot = array[start]; while(l < r) { if (array[l] <= pivot) { l++; } else { while(l < r && array[r] >= pivot) /*** skip superfluous swaps ***/ { --r; } swap(array[l], array[r], tmp_double); swap(rank[l], rank[r], tmp_int); } } if(l != end || array[end] > pivot ) { l--; swap(array[start], array[l], tmp_double); swap(rank[start], rank[l], tmp_int); qsort_index(array, start, l, rank); qsort_index(array, r, end, rank); } else { /* the first element is the largest one */ swap(array[start], array[l], tmp_double); swap(rank[start], rank[l], tmp_int); qsort_index(array, start, l-1, rank); } } } #define PI M_PI int gs_transfer_to_mosaic(char *old_file, char *mosaic_dir) { char *pch=NULL, history[512], entry[1280]; char tilefile[MAXTILE][STRING], tiletype[MAXTILE][SHORTSTRING]; char tile_name[MAXTILE][STRING]; int ntiles=0, nfiles=0, ncontact=0; double periodx=0, periody=0; int contact_tile1[MAXCONTACT], contact_tile2[MAXCONTACT]; int contact_tile1_istart[MAXCONTACT], contact_tile1_iend[MAXCONTACT]; int contact_tile1_jstart[MAXCONTACT], contact_tile1_jend[MAXCONTACT]; int contact_tile2_istart[MAXCONTACT], contact_tile2_iend[MAXCONTACT]; int contact_tile2_jstart[MAXCONTACT], contact_tile2_jend[MAXCONTACT]; char mosaic_name[128] = "atmos_mosaic"; char grid_descriptor[128] = ""; int c, i, n, m, l, errflg, check=0; char atmos_name[128]="atmos", land_name[128]="land", ocean_name[128]="ocean"; char agrid_file[128], lgrid_file[128], ogrid_file[128]; int is_coupled_grid = 0, is_ocean_only =1; int interp_order=1; double *x, *y, *dx, *dy, *area, *angle_dx, *angle_dy; int nx, ny, nxp, nyp; const char mosaic_version[] = "0.2"; if(!old_file) mpp_error("Usage:\n \t \t transfer_to_mosaic --input_file input_file.nc"); if(mpp_field_exist(old_file, "AREA_ATMxOCN") ) is_coupled_grid = 1; if(mpp_field_exist(old_file, "AREA_ATM") ) is_ocean_only = 0; if(mpp_field_exist(old_file, "DI_ATMxOCN") ) interp_order= 2; // -----------------ocean_hgrid.nc---------------------------------[ // Ocean { int nx_C, ny_C, nx_T, ny_T, nvertex; int fid_old, vid; double *x_C, *y_C, *x_T, *y_T; double *x_vert_T, *y_vert_T; double *ds_00_02_C, *ds_00_20_C, *ds_01_11_C, *ds_01_21_C, *ds_02_22_C, *ds_10_11_C, *ds_10_12_C, *ds_11_21_C, *ds_11_12_C, *ds_20_22_C, *ds_01_11_T, *ds_01_21_T, *ds_02_22_T, *ds_10_11_T, *ds_10_12_T, *ds_11_12_T, *ds_11_21_T, *ds_20_22_T, *ds_01_21_E, *ds_10_12_N, *angle_C; int i,j,k, i1, i2, j1, j2, k0, k1, k2, k3; int ji, ji0, ji1, ji2, ji3; int j1i1, j1i2, j2i1, j2i2; double unknown = 0.0; /* double unknown = -1.0E+10; */ fid_old = mpp_open(old_file, MPP_READ); nx_C = mpp_get_dimlen(fid_old, "grid_x_C"); ny_C = mpp_get_dimlen(fid_old, "grid_y_C"); nx_T = mpp_get_dimlen(fid_old, "grid_x_T"); ny_T = mpp_get_dimlen(fid_old, "grid_y_T"); nvertex = mpp_get_dimlen(fid_old, "vertex"); printf("\nReading file: %s \n", old_file); /* Reading Ocean Grid Variables */ x_C = (double *) malloc(ny_C * nx_C * sizeof(double)); y_C = (double *) malloc(ny_C * nx_C * sizeof(double)); x_T = (double *) malloc(ny_T * nx_T * sizeof(double)); y_T = (double *) malloc(ny_T * nx_T * sizeof(double)); x_vert_T = (double *) malloc(nvertex*ny_T*nx_T*sizeof(double)); y_vert_T = (double *) malloc(nvertex*ny_T*nx_T*sizeof(double)); ds_01_11_C = (double *) malloc(ny_C * nx_C * sizeof(double)); ds_10_11_C = (double *) malloc(ny_C * nx_C * sizeof(double)); ds_11_21_C = (double *) malloc(ny_C * nx_C * sizeof(double)); ds_11_12_C = (double *) malloc(ny_C * nx_C * sizeof(double)); ds_01_11_T = (double *) malloc(ny_T * nx_T * sizeof(double)); ds_10_11_T = (double *) malloc(ny_T * nx_T * sizeof(double)); ds_11_12_T = (double *) malloc(ny_T * nx_T * sizeof(double)); ds_11_21_T = (double *) malloc(ny_T * nx_T * sizeof(double)); ds_02_22_T = (double *) malloc(ny_T * nx_T * sizeof(double)); ds_20_22_T = (double *) malloc(ny_T * nx_T * sizeof(double)); angle_C = (double *) malloc(ny_C * nx_C * sizeof(double)); vid = mpp_get_varid(fid_old, "x_C"); mpp_get_var_value(fid_old, vid, x_C); vid = mpp_get_varid(fid_old, "y_C"); mpp_get_var_value(fid_old, vid, y_C); vid = mpp_get_varid(fid_old, "x_T"); mpp_get_var_value(fid_old, vid, x_T); vid = mpp_get_varid(fid_old, "y_T"); mpp_get_var_value(fid_old, vid, y_T); vid = mpp_get_varid(fid_old, "x_vert_T"); mpp_get_var_value(fid_old, vid, x_vert_T); vid = mpp_get_varid(fid_old, "y_vert_T"); mpp_get_var_value(fid_old, vid, y_vert_T); vid = mpp_get_varid(fid_old, "ds_01_11_C"); mpp_get_var_value(fid_old, vid, ds_01_11_C); vid = mpp_get_varid(fid_old, "ds_10_11_C"); mpp_get_var_value(fid_old, vid, ds_10_11_C); vid = mpp_get_varid(fid_old, "ds_11_21_C"); mpp_get_var_value(fid_old, vid, ds_11_21_C); vid = mpp_get_varid(fid_old, "ds_11_12_C"); mpp_get_var_value(fid_old, vid, ds_11_12_C); vid = mpp_get_varid(fid_old, "ds_01_11_T"); mpp_get_var_value(fid_old, vid, ds_01_11_T); vid = mpp_get_varid(fid_old, "ds_10_11_T"); mpp_get_var_value(fid_old, vid, ds_10_11_T); vid = mpp_get_varid(fid_old, "ds_11_21_T"); mpp_get_var_value(fid_old, vid, ds_11_21_T); vid = mpp_get_varid(fid_old, "ds_11_12_T"); mpp_get_var_value(fid_old, vid, ds_11_12_T); vid = mpp_get_varid(fid_old, "ds_02_22_T"); mpp_get_var_value(fid_old, vid, ds_02_22_T); vid = mpp_get_varid(fid_old, "ds_20_22_T"); mpp_get_var_value(fid_old, vid, ds_20_22_T); vid = mpp_get_varid(fid_old, "angle_C"); mpp_get_var_value(fid_old, vid, angle_C); mpp_close(fid_old); nx = 2*nx_T; nxp = nx +1; ny = 2*ny_T; nyp = ny +1; x = (double *) malloc(nxp*nyp * sizeof(double)); y = (double *) malloc(nxp*nyp* sizeof(double)); dx = (double *) malloc(nx*nyp * sizeof(double)); dy = (double *) malloc(nxp*ny * sizeof(double)); area = (double *) malloc(nx*ny * sizeof(double)); angle_dx = (double *) malloc(nxp*nyp * sizeof(double)); for(j = 0; j0) dx[(2*j+2)*nx + (2*i)] = ds_11_21_C[j*nx_T + (i-1)]; dx[(2*j+2)*nx + (2*i+1)] = ds_01_11_C[j*nx_T + i]; dy[(2*j )*nxp + (2*i+1)] = ds_10_11_T[j*nx_T + i]; if(j>0) dy[(2*j )*nxp + (2*i+2)] = ds_11_12_C[(j-1)*nx_T + i]; dy[(2*j+1)*nxp + (2*i+1)] = ds_11_12_T[j*nx_T + i]; dy[(2*j+1)*nxp + (2*i+2)] = ds_10_11_C[j*nx_T + i]; angle_dx[(2*j+1)*nxp + 2*i+1] = unknown; angle_dx[(2*j+1)*nxp + 2*i+2] = unknown; angle_dx[(2*j+2)*nxp + 2*i+1] = unknown; angle_dx[(2*j+2)*nxp + 2*i+2] = angle_C[j*nx_C + i ]; if (j==0) { x[i1] = 0.5*(x_vert_T[ji0]+x_vert_T[ji1]); x[i2] = x_vert_T[ji1]; y[i1] = 0.5*(y_vert_T[ji0]+y_vert_T[ji1]); y[i2] = y_vert_T[ji1]; dx[ (2*i+0)] = unknown; dx[ (2*i+1)] = unknown; /* dy[(2*j )*nxp + (2*i+2)] = unknown; */ dy[(2*j )*nxp + (2*i+2)] = ds_20_22_T[j*nx_T + i] - ds_10_11_C[j*nx_T + i]; angle_dx[2*i+1] = unknown; angle_dx[2*i+2] = unknown; } if (i==0) { x[j1*nxp] = 0.5*(x_vert_T[ji0]+x_vert_T[ji3]); x[j2*nxp] = x_vert_T[ji3]; y[j1*nxp] = 0.5*(y_vert_T[ji0]+y_vert_T[ji3]); y[j2*nxp] = y_vert_T[ji3]; /* dx[(2*j+2)*nx + (2*i)] = unknown;*/ dx[(2*j+2)*nx + (2*i)] = ds_02_22_T[j*nx_T + i] - ds_01_11_C[j*nx_T + i]; dy[(2*j )*nxp ] = unknown; dy[(2*j+1)*nxp ] = unknown; angle_dx[(2*j+1)*nxp] = unknown; angle_dx[(2*j+2)*nxp] = unknown; } if (i==0 && j==0) { x[0] = x_vert_T[ji0]; y[0] = y_vert_T[ji0]; angle_dx[0] = unknown; } } } //for(j = 0; j0.0) printf("%d %d %f\n ",i,j,angle_dx[j*nxp+i]); } } //for(j = 0; j85.) printf("%f ",y_T[(ny_T-1)*nx_T+i]); } //printf("\n\n"); //for(j = 0; j85.) printf("%f ",y[(nyp-1)*nxp+i]); } /* write ocean-grid data */ /*-------ocean_hgrid.nc-------*/ { int fid, id_tile, id_x, id_y, id_dx, id_dy, id_area, id_angle_dx, id_angle_dy, id_arcx; int dimlist[5], dims[2], i, j, l, ni, nj, nip, njp; char tile_spec_version[] = "0.2", geometry[] = "spherical", discretization[] = "logically_rectangular", conformal[]="true"; char outfile[128] = "", tilename[128]=""; size_t start[4], nwrite[4]; sprintf(outfile, "%s_hgrid.nc",ocean_name); printf("Writing %s\n", outfile); printf("\t nx_C=%d, ny_C=%d, x_T=%d, ny_T=%d \n", nx_C, ny_C,nx_T, ny_T); sprintf(tilename, "tile%d", 1); fid = mpp_open(outfile, MPP_WRITE); dimlist[0] = mpp_def_dim(fid, "string", STRINGLEN); dimlist[1] = mpp_def_dim(fid, "nxp", nxp); dimlist[2] = mpp_def_dim(fid, "nyp", nyp); dimlist[3] = mpp_def_dim(fid, "nx", nx); dimlist[4] = mpp_def_dim(fid, "ny", ny); id_tile = mpp_def_var(fid, "tile", MPP_CHAR, 1, dimlist, 5, "standard_name", "grid_tile_spec", "tile_spec_version", tile_spec_version, "geometry", geometry, "discretization", discretization, "conformal", conformal ); dims[0] = dimlist[2]; dims[1] = dimlist[1]; id_x = mpp_def_var(fid, "x", MPP_DOUBLE, 2, dims, 2, "standard_name", "geographic_longitude", "units", "degree_east"); id_y = mpp_def_var(fid, "y", MPP_DOUBLE, 2, dims, 2, "standard_name", "geographic_latitude", "units", "degree_north"); dims[0] = dimlist[2]; dims[1] = dimlist[3]; id_dx = mpp_def_var(fid, "dx", MPP_DOUBLE, 2, dims, 2, "standard_name", "grid_edge_x_distance", "units", "meters"); dims[0] = dimlist[4]; dims[1] = dimlist[1]; id_dy = mpp_def_var(fid, "dy", MPP_DOUBLE, 2, dims, 2, "standard_name", "grid_edge_y_distance", "units", "meters"); dims[0] = dimlist[2]; dims[1] = dimlist[1]; id_angle_dx= mpp_def_var(fid, "angle_dx", MPP_DOUBLE, 2, dims, 2, "standard_name", "grid_vertex_x_angle_WRT_geographic_east", "units", "degrees_east"); dims[0] = dimlist[4]; dims[1] = dimlist[3]; id_area= mpp_def_var(fid, "area", MPP_DOUBLE, 2, dims, 2, "standard_name", "grid_cell_area", "units", "m2"); mpp_end_def(fid); for(i=0; i<4; i++) { start[i] = 0; nwrite[i] = 0; } nwrite[0] = strlen(tilename); mpp_put_var_value_block(fid, id_tile, start, nwrite, tilename ); mpp_put_var_value(fid, id_x, x); mpp_put_var_value(fid, id_y, y); mpp_put_var_value(fid, id_dx, dx); mpp_put_var_value(fid, id_dy, dy); mpp_put_var_value(fid, id_angle_dx, angle_dx); mpp_close(fid); } } // ----------------------------------------------------------------] // -----------------atmos_hgrid.nc---------------------------------[ // Atmosphere int nxba, nyba, nxta, nyta, nxap, nyap; if(is_coupled_grid) { double *x_atmos, *y_atmos; double *xba, *yba, *xta, *yta; int i,j, fid_old, vid; fid_old = mpp_open(old_file, MPP_READ); nxba = mpp_get_dimlen(fid_old, "xba"); nyba = mpp_get_dimlen(fid_old, "yba"); nxta = mpp_get_dimlen(fid_old, "xta"); nyta = mpp_get_dimlen(fid_old, "yta"); xba = (double *) malloc(nxba * sizeof(double)); yba = (double *) malloc(nyba * sizeof(double)); xta = (double *) malloc(nxta * sizeof(double)); yta = (double *) malloc(nyta * sizeof(double)); vid = mpp_get_varid(fid_old, "xba"); mpp_get_var_value(fid_old, vid, xba); vid = mpp_get_varid(fid_old, "yba"); mpp_get_var_value(fid_old, vid, yba); vid = mpp_get_varid(fid_old, "xta"); mpp_get_var_value(fid_old, vid, xta); vid = mpp_get_varid(fid_old, "yta"); mpp_get_var_value(fid_old, vid, yta); mpp_close(fid_old); nxap = nxba + nxta; nyap = nyba + nyta; x_atmos = (double *) malloc(nxap*nyap * sizeof(double)); y_atmos = (double *) malloc(nxap*nyap * sizeof(double)); for(i = 0; i0) dim_ncontact = mpp_def_dim(fid, "ncontact", ncontact); dim_string = mpp_def_dim(fid, "string", STRING); id_mosaic = mpp_def_var(fid, "mosaic", MPP_CHAR, 1, &dim_string, 5, "standard_name", "grid_mosaic_spec", "mosaic_spec_version", mosaic_version, "children", "gridtiles", "contact_regions", "contacts", "grid_descriptor", grid_descriptor); dim[0] = dim_ntiles; dim[1] = dim_string; id_griddir = mpp_def_var(fid, "gridlocation", MPP_CHAR, 1, &dim[1], 1, "standard_name", "grid_file_location"); id_gridfiles = mpp_def_var(fid, "gridfiles", MPP_CHAR, 2, dim, 0); id_gridtiles = mpp_def_var(fid, "gridtiles", MPP_CHAR, 2, dim, 0); if(ncontact>0) { dim[0] = dim_ncontact; dim[1] = dim_string; id_contacts = mpp_def_var(fid, "contacts", MPP_CHAR, 2, dim, 5, "standard_name", "grid_contact_spec", "contact_type", "boundary", "alignment", "true", "contact_index", "contact_index", "orientation", "orient"); id_contact_index = mpp_def_var(fid, "contact_index", MPP_CHAR, 2, dim, 1, "standard_name", "starting_ending_point_index_of_contact"); } //mpp_def_global(fid, "history", history); mpp_end_def(fid); /* write out data */ for(i=0; i<4; i++) { start[i] = 0; nwrite[i] = 1; } nwrite[0] = strlen(mosaic_name); mpp_put_var_value_block(fid, id_mosaic, start, nwrite, mosaic_name); nwrite[0] = strlen(mosaic_dir); mpp_put_var_value_block(fid, id_griddir, start, nwrite, mosaic_dir); nwrite[0] = 1; for(n=0; n0) dim_ncontact = mpp_def_dim(fid, "ncontact", ncontact); dim_string = mpp_def_dim(fid, "string", STRING); id_mosaic = mpp_def_var(fid, "mosaic", MPP_CHAR, 1, &dim_string, 5, "standard_name", "grid_mosaic_spec", "mosaic_spec_version", mosaic_version, "children", "gridtiles", "contact_regions", "contacts", "grid_descriptor", grid_descriptor); dim[0] = dim_ntiles; dim[1] = dim_string; id_griddir = mpp_def_var(fid, "gridlocation", MPP_CHAR, 1, &dim[1], 1, "standard_name", "grid_file_location"); id_gridfiles = mpp_def_var(fid, "gridfiles", MPP_CHAR, 2, dim, 0); id_gridtiles = mpp_def_var(fid, "gridtiles", MPP_CHAR, 2, dim, 0); if(ncontact>0) { dim[0] = dim_ncontact; dim[1] = dim_string; id_contacts = mpp_def_var(fid, "contacts", MPP_CHAR, 2, dim, 5, "standard_name", "grid_contact_spec", "contact_type", "boundary", "alignment", "true", "contact_index", "contact_index", "orientation", "orient"); id_contact_index = mpp_def_var(fid, "contact_index", MPP_CHAR, 2, dim, 1, "standard_name", "starting_ending_point_index_of_contact"); } //mpp_def_global(fid, "history", history); mpp_end_def(fid); /* write out data */ for(i=0; i<4; i++) { start[i] = 0; nwrite[i] = 1; } nwrite[0] = strlen(mosaic_name); mpp_put_var_value_block(fid, id_mosaic, start, nwrite, mosaic_name); nwrite[0] = strlen(mosaic_dir); mpp_put_var_value_block(fid, id_griddir, start, nwrite, mosaic_dir); nwrite[0] = 1; for(n=0; n0) dim_ncontact = mpp_def_dim(fid, "ncontact", ncontact); dim_string = mpp_def_dim(fid, "string", STRING); id_mosaic = mpp_def_var(fid, "mosaic", MPP_CHAR, 1, &dim_string, 5, "standard_name", "grid_mosaic_spec", "mosaic_spec_version", mosaic_version, "children", "gridtiles", "contact_regions", "contacts", "grid_descriptor", grid_descriptor); dim[0] = dim_ntiles; dim[1] = dim_string; id_griddir = mpp_def_var(fid, "gridlocation", MPP_CHAR, 1, &dim[1], 1, "standard_name", "grid_file_location"); id_gridfiles = mpp_def_var(fid, "gridfiles", MPP_CHAR, 2, dim, 0); id_gridtiles = mpp_def_var(fid, "gridtiles", MPP_CHAR, 2, dim, 0); if(ncontact>0) { dim[0] = dim_ncontact; dim[1] = dim_string; id_contacts = mpp_def_var(fid, "contacts", MPP_CHAR, 2, dim, 5, "standard_name", "grid_contact_spec", "contact_type", "boundary", "alignment", "true", "contact_index", "contact_index", "orientation", "orient"); id_contact_index = mpp_def_var(fid, "contact_index", MPP_CHAR, 2, dim, 1, "standard_name", "starting_ending_point_index_of_contact"); } //mpp_def_global(fid, "history", history); mpp_end_def(fid); /* write out data */ for(i=0; i<4; i++) { start[i] = 0; nwrite[i] = 1; } nwrite[0] = strlen(mosaic_name); mpp_put_var_value_block(fid, id_mosaic, start, nwrite, mosaic_name); nwrite[0] = strlen(mosaic_dir); mpp_put_var_value_block(fid, id_griddir, start, nwrite, mosaic_dir); nwrite[0] = 1; for(n=0; n0) { size_t start[4], nwrite[4]; int fid, dim_string, dim_ncells, dim_two, dims[4]; int id_xgrid_area, id_contact, n; int id_tile1_cell, id_tile2_cell, id_tile1_dist, id_tile2_dist; char contact[STRING]; for(i=0; i<4; i++) { start[i] = 0; nwrite[i] = 1; } sprintf(axl_file, "%s_mosaicX%s_mosaic.nc", amosaic_name, lmosaic_name); printf("Writing %s_mosaicX%s_mosaic.nc\n", amosaic_name, lmosaic_name); sprintf(contact, "%s_mosaic:%s::%s_mosaic:%s", amosaic_name, atile_name, lmosaic_name, ltile_name); fid = mpp_open(axl_file, MPP_WRITE); //mpp_def_global(fid, "history", history); dim_string = mpp_def_dim(fid, "string", STRING); dim_ncells = mpp_def_dim(fid, "ncells", naxl); dim_two = mpp_def_dim(fid, "two", 2); if(interp_order == 2) { id_contact = mpp_def_var(fid, "contact", MPP_CHAR, 1, &dim_string, 8, "standard_name", "grid_contact_spec", "contact_spec_version", version, "contact_type", "exchange", "parent1_cell", "tile1_cell", "parent2_cell", "tile2_cell", "xgrid_area_field", "xgrid_area", "distant_to_parent1_centroid", "tile1_distance", "distant_to_parent2_centroid", "tile2_distance"); } else { id_contact = mpp_def_var(fid, "contact", MPP_CHAR, 1, &dim_string, 6, "standard_name", "grid_contact_spec", "contact_spec_version", version, "contact_type", "exchange", "parent1_cell", "tile1_cell", "parent2_cell", "tile2_cell", "xgrid_area_field", "xgrid_area"); } dims[0] = dim_ncells; dims[1] = dim_two; id_tile1_cell = mpp_def_var(fid, "tile1_cell", MPP_INT, 2, dims, 1, "standard_name", "parent_cell_indices_in_mosaic1"); id_tile2_cell = mpp_def_var(fid, "tile2_cell", MPP_INT, 2, dims, 1, "standard_name", "parent_cell_indices_in_mosaic2"); id_xgrid_area = mpp_def_var(fid, "xgrid_area", MPP_DOUBLE, 1, &dim_ncells, 2, "standard_name", "exchange_grid_area", "units", "m2"); if(interp_order == 2) { id_tile1_dist = mpp_def_var(fid, "tile1_distance", MPP_DOUBLE, 2, dims, 1, "standard_name", "distance_from_parent1_cell_centroid"); id_tile2_dist = mpp_def_var(fid, "tile2_distance", MPP_DOUBLE, 2, dims, 1, "standard_name", "distance_from_parent2_cell_centroid"); } mpp_end_def(fid); for(i=0; i<4; i++) {start[i] = 0; nwrite[i] = 1;} nwrite[0] = strlen(contact); mpp_put_var_value_block(fid, id_contact, start, nwrite, contact); nwrite[0] = naxl; atmxlnd_area = (double *)malloc(naxl*sizeof(double)); atmxlnd_ia = (int *)malloc(naxl*sizeof(int )); atmxlnd_ja = (int *)malloc(naxl*sizeof(int )); atmxlnd_il = (int *)malloc(naxl*sizeof(int )); atmxlnd_jl = (int *)malloc(naxl*sizeof(int )); vid = mpp_get_varid(fid_old, "AREA_ATMxLND"); mpp_get_var_value(fid_old, vid, atmxlnd_area); vid = mpp_get_varid(fid_old, "I_ATM_ATMxLND"); mpp_get_var_value(fid_old, vid, atmxlnd_ia); vid = mpp_get_varid(fid_old, "J_ATM_ATMxLND"); mpp_get_var_value(fid_old, vid, atmxlnd_ja); vid = mpp_get_varid(fid_old, "I_LND_ATMxLND"); mpp_get_var_value(fid_old, vid, atmxlnd_il); vid = mpp_get_varid(fid_old, "J_LND_ATMxLND"); mpp_get_var_value(fid_old, vid, atmxlnd_jl); for(i=0;i0) { size_t start[4], nwrite[4]; int fid, dim_string, dim_ncells, dim_two, dims[4]; int id_xgrid_area, id_contact, n; int id_tile1_cell, id_tile2_cell, id_tile1_dist, id_tile2_dist; char contact[STRING]; for(i=0; i<4; i++) { start[i] = 0; nwrite[i] = 1; } sprintf(axo_file, "%s_mosaicX%s_mosaic.nc", amosaic_name, omosaic_name); printf("Writing %s_mosaicX%s_mosaic.nc\n", amosaic_name, omosaic_name); sprintf(contact, "%s_mosaic:%s::%s_mosaic:%s", amosaic_name, atile_name, omosaic_name, otile_name); fid = mpp_open(axo_file, MPP_WRITE); //mpp_def_global(fid, "history", history); dim_string = mpp_def_dim(fid, "string", STRING); dim_ncells = mpp_def_dim(fid, "ncells", naxo); dim_two = mpp_def_dim(fid, "two", 2); if(interp_order == 2) { id_contact = mpp_def_var(fid, "contact", MPP_CHAR, 1, &dim_string, 8, "standard_name", "grid_contact_spec", "contact_spec_version", version, "contact_type", "exchange", "parent1_cell", "tile1_cell", "parent2_cell", "tile2_cell", "xgrid_area_field", "xgrid_area", "distant_to_parent1_centroid", "tile1_distance", "distant_to_parent2_centroid", "tile2_distance"); } else { id_contact = mpp_def_var(fid, "contact", MPP_CHAR, 1, &dim_string, 6, "standard_name", "grid_contact_spec", "contact_spec_version", version, "contact_type", "exchange", "parent1_cell", "tile1_cell", "parent2_cell", "tile2_cell", "xgrid_area_field", "xgrid_area"); } dims[0] = dim_ncells; dims[1] = dim_two; id_tile1_cell = mpp_def_var(fid, "tile1_cell", MPP_INT, 2, dims, 1, "standard_name", "parent_cell_indices_in_mosaic1"); id_tile2_cell = mpp_def_var(fid, "tile2_cell", MPP_INT, 2, dims, 1, "standard_name", "parent_cell_indices_in_mosaic2"); id_xgrid_area = mpp_def_var(fid, "xgrid_area", MPP_DOUBLE, 1, &dim_ncells, 2, "standard_name", "exchange_grid_area", "units", "m2"); if(interp_order == 2) { id_tile1_dist = mpp_def_var(fid, "tile1_distance", MPP_DOUBLE, 2, dims, 1, "standard_name", "distance_from_parent1_cell_centroid"); id_tile2_dist = mpp_def_var(fid, "tile2_distance", MPP_DOUBLE, 2, dims, 1, "standard_name", "distance_from_parent2_cell_centroid"); } mpp_end_def(fid); for(i=0; i<4; i++) { start[i] = 0; nwrite[i] = 1;} nwrite[0] = strlen(contact); mpp_put_var_value_block(fid, id_contact, start, nwrite, contact); nwrite[0] = naxo; atmxocn_area = (double *)malloc(naxo*sizeof(double)); atmxocn_ia = (int *)malloc(naxo*sizeof(int )); atmxocn_ja = (int *)malloc(naxo*sizeof(int )); atmxocn_io = (int *)malloc(naxo*sizeof(int )); atmxocn_jo = (int *)malloc(naxo*sizeof(int )); vid = mpp_get_varid(fid_old, "AREA_ATMxOCN"); mpp_get_var_value(fid_old, vid, atmxocn_area); vid = mpp_get_varid(fid_old, "I_ATM_ATMxOCN"); mpp_get_var_value(fid_old, vid, atmxocn_ia); vid = mpp_get_varid(fid_old, "J_ATM_ATMxOCN"); mpp_get_var_value(fid_old, vid, atmxocn_ja); vid = mpp_get_varid(fid_old, "I_OCN_ATMxOCN"); mpp_get_var_value(fid_old, vid, atmxocn_io); vid = mpp_get_varid(fid_old, "J_OCN_ATMxOCN"); mpp_get_var_value(fid_old, vid, atmxocn_jo); for(i=0;i0) { size_t start[4], nwrite[4]; int fid, dim_string, dim_ncells, dim_two, dims[4]; int id_xgrid_area, id_contact, n; int id_tile1_cell, id_tile2_cell, id_tile1_dist, id_tile2_dist; char contact[STRING]; for(i=0; i<4; i++) { start[i] = 0; nwrite[i] = 1; } sprintf(lxo_file, "%s_mosaicX%s_mosaic.nc", lmosaic_name, omosaic_name); printf("Writing %s_mosaicX%s_mosaic.nc\n", lmosaic_name, omosaic_name); sprintf(contact, "%s_mosaic:%s::%s_mosaic:%s", lmosaic_name, ltile_name, omosaic_name, otile_name); fid = mpp_open(lxo_file, MPP_WRITE); //mpp_def_global(fid, "history", history); dim_string = mpp_def_dim(fid, "string", STRING); dim_ncells = mpp_def_dim(fid, "ncells", nlxo); dim_two = mpp_def_dim(fid, "two", 2); if(interp_order == 2) { id_contact = mpp_def_var(fid, "contact", MPP_CHAR, 1, &dim_string, 8, "standard_name", "grid_contact_spec", "contact_spec_version", version, "contact_type", "exchange", "parent1_cell", "tile1_cell", "parent2_cell", "tile2_cell", "xgrid_area_field", "xgrid_area", "distant_to_parent1_centroid", "tile1_distance", "distant_to_parent2_centroid", "tile2_distance"); } else { id_contact = mpp_def_var(fid, "contact", MPP_CHAR, 1, &dim_string, 6, "standard_name", "grid_contact_spec", "contact_spec_version", version, "contact_type", "exchange", "parent1_cell", "tile1_cell", "parent2_cell", "tile2_cell", "xgrid_area_field", "xgrid_area"); } dims[0] = dim_ncells; dims[1] = dim_two; id_tile1_cell = mpp_def_var(fid, "tile1_cell", MPP_INT, 2, dims, 1, "standard_name", "parent_cell_indices_in_mosaic1"); id_tile2_cell = mpp_def_var(fid, "tile2_cell", MPP_INT, 2, dims, 1, "standard_name", "parent_cell_indices_in_mosaic2"); id_xgrid_area = mpp_def_var(fid, "xgrid_area", MPP_DOUBLE, 1, &dim_ncells, 2, "standard_name", "exchange_grid_area", "units", "m2"); if(interp_order == 2) { id_tile1_dist = mpp_def_var(fid, "tile1_distance", MPP_DOUBLE, 2, dims, 1, "standard_name", "distance_from_parent1_cell_centroid"); id_tile2_dist = mpp_def_var(fid, "tile2_distance", MPP_DOUBLE, 2, dims, 1, "standard_name", "distance_from_parent2_cell_centroid"); } mpp_end_def(fid); for(i=0; i<4; i++) { start[i] = 0; nwrite[i] = 1;} nwrite[0] = strlen(contact); mpp_put_var_value_block(fid, id_contact, start, nwrite, contact); nwrite[0] = nlxo; lndxocn_area = (double *)malloc(nlxo*sizeof(double)); lndxocn_il = (int *)malloc(nlxo*sizeof(int )); lndxocn_jl = (int *)malloc(nlxo*sizeof(int )); lndxocn_io = (int *)malloc(nlxo*sizeof(int )); lndxocn_jo = (int *)malloc(nlxo*sizeof(int )); vid = mpp_get_varid(fid_old, "AREA_LNDxOCN"); mpp_get_var_value(fid_old, vid, lndxocn_area); vid = mpp_get_varid(fid_old, "I_LND_LNDxOCN"); mpp_get_var_value(fid_old, vid, lndxocn_il); vid = mpp_get_varid(fid_old, "J_LND_LNDxOCN"); mpp_get_var_value(fid_old, vid, lndxocn_jl); vid = mpp_get_varid(fid_old, "I_OCN_LNDxOCN"); mpp_get_var_value(fid_old, vid, lndxocn_io); vid = mpp_get_varid(fid_old, "J_OCN_LNDxOCN"); mpp_get_var_value(fid_old, vid, lndxocn_jo); for(i=0;i min_lat) min_atm_lat = min_lat; } if(tmpy[0]*D2R > min_atm_lat + TINY_VALUE) { /* extend one point in south direction*/ ocn_south_ext = 1; } } nxo[n] /= x_refine; nyo[n] /= y_refine; nyo[n] += ocn_south_ext; xocn[n] = (double *)malloc((nxo[n]+1)*(nyo[n]+1)*sizeof(double)); yocn[n] = (double *)malloc((nxo[n]+1)*(nyo[n]+1)*sizeof(double)); area_ocn[n] = (double *)malloc((nxo[n] )*(nyo[n] )*sizeof(double)); for(j = 0; j < nyo[n]+1; j++) for(i = 0; i < nxo[n]+1; i++) { xocn[n][(j+ocn_south_ext)*(nxo[n]+1)+i] = tmpx[(j*y_refine)*(nxo[n]*x_refine+1)+i*x_refine] * D2R; yocn[n][(j+ocn_south_ext)*(nxo[n]+1)+i] = tmpy[(j*y_refine)*(nxo[n]*x_refine+1)+i*x_refine] * D2R; } if(ocn_south_ext==1) { for(i=0; isea_level) omask[n][(j+ocn_south_ext)*nx+i] = 1; } free(depth); } mpp_close(t_fid); } /* Either omosaic is different from both lmosaic and amosaic, or all the three mosaic are the same */ if(strcmp(omosaic_name, amosaic_name)) { /* omosaic is different from amosaic */ if(!strcmp(omosaic_name, lmosaic_name)) mpp_error("make_coupler_mosaic: omosaic is the same as lmosaic, " "but different from amosaic."); same_mosaic = 0; } else { /* omosaic is same as amosaic */ if(strcmp(omosaic_name, lmosaic_name)) mpp_error("make_coupler_mosaic: omosaic is the same as amosaic, " "but different from lmosaic."); same_mosaic = 1; } /*************************************************************************************** First generate the exchange grid between atmos mosaic and land/ocean mosaic ***************************************************************************************/ nfile_axo = 0; nfile_axl = 0; nfile_lxo = 0; { int no, nl, na, n; size_t **naxl, **naxo; int ***atmxlnd_ia, ***atmxlnd_ja, ***atmxlnd_il, ***atmxlnd_jl; int ***atmxocn_ia, ***atmxocn_ja, ***atmxocn_io, ***atmxocn_jo; double ***atmxlnd_area, ***atmxlnd_dia, ***atmxlnd_dja, ***atmxlnd_dil, ***atmxlnd_djl; double ***atmxocn_area, ***atmxocn_dia, ***atmxocn_dja, ***atmxocn_dio, ***atmxocn_djo; double ***atmxocn_clon, ***atmxocn_clat, ***atmxlnd_clon, ***atmxlnd_clat; double min_area; naxl = (size_t ** )malloc(ntile_atm*sizeof(size_t *)); naxo = (size_t ** )malloc(ntile_atm*sizeof(size_t *)); atmxlnd_area = (double ***)malloc(ntile_atm*sizeof(double **)); atmxlnd_ia = (int ***)malloc(ntile_atm*sizeof(int **)); atmxlnd_ja = (int ***)malloc(ntile_atm*sizeof(int **)); atmxlnd_il = (int ***)malloc(ntile_atm*sizeof(int **)); atmxlnd_jl = (int ***)malloc(ntile_atm*sizeof(int **)); atmxocn_area = (double ***)malloc(ntile_atm*sizeof(double **)); atmxocn_ia = (int ***)malloc(ntile_atm*sizeof(int **)); atmxocn_ja = (int ***)malloc(ntile_atm*sizeof(int **)); atmxocn_io = (int ***)malloc(ntile_atm*sizeof(int **)); atmxocn_jo = (int ***)malloc(ntile_atm*sizeof(int **)); if(interp_order == 2 ) { atmxlnd_dia = (double ***)malloc(ntile_atm*sizeof(double **)); atmxlnd_dja = (double ***)malloc(ntile_atm*sizeof(double **)); atmxlnd_dil = (double ***)malloc(ntile_atm*sizeof(double **)); atmxlnd_djl = (double ***)malloc(ntile_atm*sizeof(double **)); atmxocn_dia = (double ***)malloc(ntile_atm*sizeof(double **)); atmxocn_dja = (double ***)malloc(ntile_atm*sizeof(double **)); atmxocn_dio = (double ***)malloc(ntile_atm*sizeof(double **)); atmxocn_djo = (double ***)malloc(ntile_atm*sizeof(double **)); atmxlnd_clon = (double ***)malloc(ntile_atm*sizeof(double **)); atmxlnd_clat = (double ***)malloc(ntile_atm*sizeof(double **)); atmxocn_clon = (double ***)malloc(ntile_atm*sizeof(double **)); atmxocn_clat = (double ***)malloc(ntile_atm*sizeof(double **)); } for(na=0; na= ya_max || yl_max <= ya_min ) continue; nl_in = fix_lon(xl, yl, 4, xa_avg); xl_min = minval_double(nl_in, xl); xl_max = maxval_double(nl_in, xl); /* xl should in the same range as xa after lon_fix, so no need to consider cyclic condition */ if(xa_min >= xl_max || xa_max <= xl_min ) continue; if ( (n_out = clip_2dx2d( xa, ya, na_in, xl, yl, nl_in, x_out, y_out )) > 0 ) { xarea = poly_area(x_out, y_out, n_out); min_area = min(area_lnd[nl][jl*nxl[nl]+il], area_atm[na][la]); if( xarea/min_area > AREA_RATIO_THRESH ) { /* remember the exchange grid vertices */ for(n=0; nMX) mpp_error("make_coupler_mosaic: count is greater than MX, increase MX"); } } } } /* calculate atmos/ocean x-cells */ for(no=0; no 0.5) { /* over sea/ice */ /* xo should in the same range as xa after lon_fix, so no need to consider cyclic condition */ if(xa_min >= xo_max || xa_max <= xo_min || yo_min >= ya_max || yo_max <= ya_min ) continue; if ( (n_out = clip_2dx2d( xa, ya, na_in, xo, yo, no_in, x_out, y_out )) > 0) { xarea = poly_area(x_out, y_out, n_out ); min_area = min(area_ocn[no][jo*nxo[no]+io], area_atm[na][la]); if(xarea/min_area > AREA_RATIO_THRESH) { atmxocn_area[na][no][naxo[na][no]] = xarea; atmxocn_io[na][no][naxo[na][no]] = io; atmxocn_jo[na][no][naxo[na][no]] = jo; atmxocn_ia[na][no][naxo[na][no]] = ia; atmxocn_ja[na][no][naxo[na][no]] = ja; if(interp_order == 2) { atmxocn_clon[na][no][naxo[na][no]] = poly_ctrlon ( x_out, y_out, n_out, xa_avg); atmxocn_clat[na][no][naxo[na][no]] = poly_ctrlat ( x_out, y_out, n_out ); } ++(naxo[na][no]); if(naxo[na][no] > MAXXGRID) mpp_error("naxo is greater than MAXXGRID, increase MAXXGRID"); } } } else { /* over land */ /* find the overlap of atmxlnd and ocean cell */ for(l=0; l= xo_max || axl_xmax[l] <= xo_min || axl_ymin[l] >= ya_max || axl_ymax[l] <= ya_min ) continue; if((n_out = clip_2dx2d( atmxlnd_x[l], atmxlnd_y[l], num_v[l], xo, yo, no_in, x_out, y_out )) > 0) { xarea = poly_area(x_out, y_out, n_out ); min_area = min(area_lnd[axl_t[l]][axl_j[l]*nxl[axl_t[l]]+axl_i[l]], area_atm[na][la]); if(xarea/min_area > AREA_RATIO_THRESH) { axl_area[l] += xarea; if(interp_order == 2) { axl_clon[l] += poly_ctrlon ( x_out, y_out, n_out, xa_avg); axl_clat[l] += poly_ctrlat ( x_out, y_out, n_out); } } } } } } } /* get the exchange grid between land and atmos. */ for(l=0; l AREA_RATIO_THRESH) { atmxlnd_area[na][nl][naxl[na][nl]] = axl_area[l]; atmxlnd_ia [na][nl][naxl[na][nl]] = ia; atmxlnd_ja [na][nl][naxl[na][nl]] = ja; atmxlnd_il [na][nl][naxl[na][nl]] = axl_i[l]; atmxlnd_jl [na][nl][naxl[na][nl]] = axl_j[l]; if(interp_order == 2) { atmxlnd_clon[na][nl][naxl[na][nl]] = axl_clon[l]; atmxlnd_clat[na][nl][naxl[na][nl]] = axl_clat[l]; } ++(naxl[na][nl]); if(naxl[na][nl] > MAXXGRID) mpp_error("naxl is greater than MAXXGRID, increase MAXXGRID"); } } }/* end of la loop */ mpp_delete_domain2d(&Dom); } /* end of na loop */ /* calculate the centroid of model grid, as well as land_mask and ocean_mask */ { double **l_area, **o_area; int nl, no, ll, lo; l_area = (double **)malloc(ntile_lnd*sizeof(double *)); o_area = (double **)malloc(ntile_ocn*sizeof(double *)); for(nl =0; nl 0) { double *g_area; int *g_il, *g_jl; int ii; g_il = (int *)malloc(nxgrid*sizeof(int )); g_jl = (int *)malloc(nxgrid*sizeof(int )); g_area = (double *)malloc(nxgrid*sizeof(double)); mpp_gather_field_int (naxl[na][nl], atmxlnd_il[na][nl], g_il); mpp_gather_field_int (naxl[na][nl], atmxlnd_jl[na][nl], g_jl); mpp_gather_field_double(naxl[na][nl], atmxlnd_area[na][nl], g_area); for(i=0; i 0) { double *g_area; int *g_io, *g_jo; int ii; g_io = (int *)malloc(nxgrid*sizeof(int )); g_jo = (int *)malloc(nxgrid*sizeof(int )); g_area = (double *)malloc(nxgrid*sizeof(double)); mpp_gather_field_int (naxo[na][no], atmxocn_io[na][no], g_io); mpp_gather_field_int (naxo[na][no], atmxocn_jo[na][no], g_jo); mpp_gather_field_double(naxo[na][no], atmxocn_area[na][no], g_area); for(i=0; i 0) { double *g_area, *g_clon, *g_clat; int *g_ia, *g_ja, *g_il, *g_jl; int ii; g_ia = (int *)malloc(nxgrid*sizeof(int )); g_ja = (int *)malloc(nxgrid*sizeof(int )); g_il = (int *)malloc(nxgrid*sizeof(int )); g_jl = (int *)malloc(nxgrid*sizeof(int )); g_area = (double *)malloc(nxgrid*sizeof(double)); g_clon = (double *)malloc(nxgrid*sizeof(double)); g_clat = (double *)malloc(nxgrid*sizeof(double)); mpp_gather_field_int (naxl[na][nl], atmxlnd_ia[na][nl], g_ia); mpp_gather_field_int (naxl[na][nl], atmxlnd_ja[na][nl], g_ja); mpp_gather_field_int (naxl[na][nl], atmxlnd_il[na][nl], g_il); mpp_gather_field_int (naxl[na][nl], atmxlnd_jl[na][nl], g_jl); mpp_gather_field_double(naxl[na][nl], atmxlnd_area[na][nl], g_area); mpp_gather_field_double(naxl[na][nl], atmxlnd_clon[na][nl], g_clon); mpp_gather_field_double(naxl[na][nl], atmxlnd_clat[na][nl], g_clat); for(i=0; i 0) { double *g_area, *g_clon, *g_clat; int *g_ia, *g_ja, *g_io, *g_jo; int ii; g_ia = (int *)malloc(nxgrid*sizeof(int )); g_ja = (int *)malloc(nxgrid*sizeof(int )); g_io = (int *)malloc(nxgrid*sizeof(int )); g_jo = (int *)malloc(nxgrid*sizeof(int )); g_area = (double *)malloc(nxgrid*sizeof(double)); g_clon = (double *)malloc(nxgrid*sizeof(double)); g_clat = (double *)malloc(nxgrid*sizeof(double)); mpp_gather_field_int (naxo[na][no], atmxocn_ia[na][no], g_ia); mpp_gather_field_int (naxo[na][no], atmxocn_ja[na][no], g_ja); mpp_gather_field_int (naxo[na][no], atmxocn_io[na][no], g_io); mpp_gather_field_int (naxo[na][no], atmxocn_jo[na][no], g_jo); mpp_gather_field_double(naxo[na][no], atmxocn_area[na][no], g_area); mpp_gather_field_double(naxo[na][no], atmxocn_clon[na][no], g_clon); mpp_gather_field_double(naxo[na][no], atmxocn_clat[na][no], g_clat); for(i=0; i 0) { a_clon[la] /= a_area[la]; a_clat[la] /= a_area[la]; } } /* substract atmos centroid to get the centroid distance between atmos grid and exchange grid. */ for(nl=0; nl 0) { l_clon[nl][ll] /= l_area[nl][ll]; l_clat[nl][ll] /= l_area[nl][ll]; } } for(na=0; na 0) { o_clon[no][lo] /= o_area[no][lo]; o_clat[no][lo] /= o_area[no][lo]; } } for(na=0; na TOLORENCE ) { printf("at ocean point (%d,%d), omask = %f, ocn_frac = %f, diff = %f\n", io, jo, omask[no][i], ocn_frac, omask[no][i] - ocn_frac); mpp_error("make_coupler_mosaic: omask is not equal ocn_frac"); } mask[jo*nxo[no]+io] = ocn_frac; } if(ntile_ocn > 1) sprintf(ocn_mask_file, "ocean_mask_tile%d.nc", no+1); else strcpy(ocn_mask_file, "ocean_mask.nc"); fid = mpp_open(ocn_mask_file, MPP_WRITE); mpp_def_global_att(fid, "grid_version", GRID_VERSION); mpp_def_global_att(fid, "code_version", TAGNAME); mpp_def_global_att(fid, "history", history); dims[1] = mpp_def_dim(fid, "nx", nxo[no]); dims[0] = mpp_def_dim(fid, "ny", ny); id_mask = mpp_def_var(fid, "mask", MPP_DOUBLE, 2, dims, 2, "standard_name", "ocean fraction at T-cell centers", "units", "none"); mpp_end_def(fid); mpp_put_var_value(fid, id_mask, mask); mpp_close(fid); free(mask); } } /* calculate land_frac and write out land_frac */ { int il, jl; int id_mask, fid, dims[2]; char lnd_mask_file[STRING]; double *mask; for(nl=0; nl 1) sprintf(lnd_mask_file, "land_mask_tile%d.nc", nl+1); else strcpy(lnd_mask_file, "land_mask.nc"); fid = mpp_open(lnd_mask_file, MPP_WRITE); mpp_def_global_att(fid, "grid_version", GRID_VERSION); mpp_def_global_att(fid, "code_version", TAGNAME); mpp_def_global_att(fid, "history", history); dims[1] = mpp_def_dim(fid, "nx", nxl[nl]); dims[0] = mpp_def_dim(fid, "ny", nyl[nl]); id_mask = mpp_def_var(fid, "mask", MPP_DOUBLE, 2, dims, 2, "standard_name", "land fraction at T-cell centers", "units", "none"); mpp_end_def(fid); mpp_put_var_value(fid, id_mask, mask); free(mask); mpp_close(fid); } } for(nl=0; nl0) { size_t start[4], nwrite[4]; int *gdata_int; double *gdata_dbl; int fid, dim_string, dim_ncells, dim_two, dims[4]; int id_xgrid_area, id_contact, n; int id_tile1_cell, id_tile2_cell, id_tile1_dist, id_tile2_dist; char contact[STRING]; for(i=0; i<4; i++) { start[i] = 0; nwrite[i] = 1; } if(same_mosaic) sprintf(axl_file[nfile_axl], "atm_%s_%sXlnd_%s_%s.nc", amosaic_name, atile_name[na], lmosaic_name, ltile_name[nl]); else sprintf(axl_file[nfile_axl], "%s_%sX%s_%s.nc", amosaic_name, atile_name[na], lmosaic_name, ltile_name[nl]); sprintf(contact, "%s:%s::%s:%s", amosaic_name, atile_name[na], lmosaic_name, ltile_name[nl]); fid = mpp_open(axl_file[nfile_axl], MPP_WRITE); mpp_def_global_att(fid, "grid_version", GRID_VERSION); mpp_def_global_att(fid, "code_version", TAGNAME); mpp_def_global_att(fid, "history", history); dim_string = mpp_def_dim(fid, "string", STRING); dim_ncells = mpp_def_dim(fid, "ncells", nxgrid); dim_two = mpp_def_dim(fid, "two", 2); if(interp_order == 2) { id_contact = mpp_def_var(fid, "contact", MPP_CHAR, 1, &dim_string, 7, "standard_name", "grid_contact_spec", "contact_type", "exchange", "parent1_cell", "tile1_cell", "parent2_cell", "tile2_cell", "xgrid_area_field", "xgrid_area", "distant_to_parent1_centroid", "tile1_distance", "distant_to_parent2_centroid", "tile2_distance"); } else { id_contact = mpp_def_var(fid, "contact", MPP_CHAR, 1, &dim_string, 5, "standard_name", "grid_contact_spec", "contact_type", "exchange", "parent1_cell", "tile1_cell", "parent2_cell", "tile2_cell", "xgrid_area_field", "xgrid_area"); } dims[0] = dim_ncells; dims[1] = dim_two; id_tile1_cell = mpp_def_var(fid, "tile1_cell", MPP_INT, 2, dims, 1, "standard_name", "parent_cell_indices_in_mosaic1"); id_tile2_cell = mpp_def_var(fid, "tile2_cell", MPP_INT, 2, dims, 1, "standard_name", "parent_cell_indices_in_mosaic2"); id_xgrid_area = mpp_def_var(fid, "xgrid_area", MPP_DOUBLE, 1, &dim_ncells, 2, "standard_name", "exchange_grid_area", "units", "m2"); if(interp_order == 2) { id_tile1_dist = mpp_def_var(fid, "tile1_distance", MPP_DOUBLE, 2, dims, 1, "standard_name", "distance_from_parent1_cell_centroid"); id_tile2_dist = mpp_def_var(fid, "tile2_distance", MPP_DOUBLE, 2, dims, 1, "standard_name", "distance_from_parent2_cell_centroid"); } mpp_end_def(fid); /* the index will start from 1, instead of 0 ( fortran index) */ for(i = 0;i < naxl[na][nl]; i++) { ++(atmxlnd_ia[na][nl][i]); ++(atmxlnd_ja[na][nl][i]); ++(atmxlnd_il[na][nl][i]); ++(atmxlnd_jl[na][nl][i]); } nwrite[0] = strlen(contact); mpp_put_var_value_block(fid, id_contact, start, nwrite, contact); nwrite[0] = nxgrid; gdata_int = (int *)malloc(nxgrid*sizeof(int)); gdata_dbl = (double *)malloc(nxgrid*sizeof(double)); mpp_gather_field_double(naxl[na][nl], atmxlnd_area[na][nl], gdata_dbl); if(check) { for(n=0; n0) { size_t start[4], nwrite[4]; int *gdata_int; double *gdata_dbl; int fid, dim_string, dim_ncells, dim_two, dims[4]; int id_xgrid_area, id_contact, n; int id_tile1_cell, id_tile2_cell, id_tile1_dist, id_tile2_dist; char contact[STRING]; for(i=0; i<4; i++) { start[i] = 0; nwrite[i] = 1; } if(same_mosaic) sprintf(axo_file[nfile_axo], "atm_%s_%sXocn_%s_%s.nc", amosaic_name, atile_name[na], omosaic_name, otile_name[no]); else sprintf(axo_file[nfile_axo], "%s_%sX%s_%s.nc", amosaic_name, atile_name[na], omosaic_name, otile_name[no]); sprintf(contact, "%s:%s::%s:%s", amosaic_name, atile_name[na], omosaic_name, otile_name[no]); fid = mpp_open(axo_file[nfile_axo], MPP_WRITE); mpp_def_global_att(fid, "grid_version", GRID_VERSION); mpp_def_global_att(fid, "code_version", TAGNAME); mpp_def_global_att(fid, "history", history); dim_string = mpp_def_dim(fid, "string", STRING); dim_ncells = mpp_def_dim(fid, "ncells", nxgrid); dim_two = mpp_def_dim(fid, "two", 2); if(interp_order == 2) { id_contact = mpp_def_var(fid, "contact", MPP_CHAR, 1, &dim_string, 7, "standard_name", "grid_contact_spec", "contact_type", "exchange", "parent1_cell", "tile1_cell", "parent2_cell", "tile2_cell", "xgrid_area_field", "xgrid_area", "distant_to_parent1_centroid", "tile1_distance", "distant_to_parent2_centroid", "tile2_distance"); } else { id_contact = mpp_def_var(fid, "contact", MPP_CHAR, 1, &dim_string, 5, "standard_name", "grid_contact_spec", "contact_type", "exchange", "parent1_cell", "tile1_cell", "parent2_cell", "tile2_cell", "xgrid_area_field", "xgrid_area" ); } dims[0] = dim_ncells; dims[1] = dim_two; id_tile1_cell = mpp_def_var(fid, "tile1_cell", MPP_INT, 2, dims, 1, "standard_name", "parent_cell_indices_in_mosaic1"); id_tile2_cell = mpp_def_var(fid, "tile2_cell", MPP_INT, 2, dims, 1, "standard_name", "parent_cell_indices_in_mosaic2"); id_xgrid_area = mpp_def_var(fid, "xgrid_area", MPP_DOUBLE, 1, &dim_ncells, 2, "standard_name", "exchange_grid_area", "units", "m2"); if(interp_order == 2) { id_tile1_dist = mpp_def_var(fid, "tile1_distance", MPP_DOUBLE, 2, dims, 1, "standard_name", "distance_from_parent1_cell_centroid"); id_tile2_dist = mpp_def_var(fid, "tile2_distance", MPP_DOUBLE, 2, dims, 1, "standard_name", "distance_from_parent2_cell_centroid"); } mpp_end_def(fid); /* the index will start from 1, instead of 0 ( fortran index) */ for(i = 0;i < naxo[na][no]; i++) { ++(atmxocn_ia[na][no][i]); ++(atmxocn_ja[na][no][i]); ++(atmxocn_io[na][no][i]); atmxocn_jo[na][no][i] += 1-ocn_south_ext; /* possible one artificial j-level is added at south end */ } nwrite[0] = strlen(contact); mpp_put_var_value_block(fid, id_contact, start, nwrite, contact); nwrite[0] = nxgrid; gdata_int = (int *)malloc(nxgrid*sizeof(int)); gdata_dbl = (double *)malloc(nxgrid*sizeof(double)); mpp_gather_field_double(naxo[na][no], atmxocn_area[na][no], gdata_dbl); if(check) { for(n=0; n 0.5) { n0 = jo *(nxo[no]+1) + io; n1 = jo *(nxo[no]+1) + io+1; n2 = (jo+1)*(nxo[no]+1) + io+1; n3 = (jo+1)*(nxo[no]+1) + io; xo[0] = xocn[no][n0]; yo[0] = yocn[no][n0]; xo[1] = xocn[no][n1]; yo[1] = yocn[no][n1]; xo[2] = xocn[no][n2]; yo[2] = yocn[no][n2]; xo[3] = xocn[no][n3]; yo[3] = yocn[no][n3]; yo_min = minval_double(4, yo); yo_max = maxval_double(4, yo); if(yo_min >= yl_max || yo_max <= yl_min ) continue; no_in = fix_lon(xo, yo, 4, xl_avg); xo_min = minval_double(no_in, xo); xo_max = maxval_double(no_in, xo); /* xo should in the same range as xa after lon_fix, so no need to consider cyclic condition */ if(xl_min >= xo_max || xl_max <= xo_min ) continue; if ( (n_out = clip_2dx2d( xl, yl, nl_in, xo, yo, no_in, x_out, y_out )) > 0 ){ xarea = poly_area(x_out, y_out, n_out ); min_area = min(area_ocn[no][jo*nxo[no]+io], area_lnd[nl][ll] ); if(xarea/min_area > AREA_RATIO_THRESH ) { lndxocn_area[nl][no][nlxo[nl][no]] = xarea; lndxocn_io[nl][no][nlxo[nl][no]] = io; lndxocn_jo[nl][no][nlxo[nl][no]] = jo; lndxocn_il[nl][no][nlxo[nl][no]] = il; lndxocn_jl[nl][no][nlxo[nl][no]] = jl; if(interp_order == 2) { lndxocn_clon[nl][no][nlxo[nl][no]] = poly_ctrlon ( x_out, y_out, n_out, xl_avg); lndxocn_clat[nl][no][nlxo[nl][no]] = poly_ctrlat ( x_out, y_out, n_out ); } ++(nlxo[nl][no]); if(nlxo[nl][no] > MAXXGRID) mpp_error("nlxo is greater than MAXXGRID, increase MAXXGRID"); } } } /* end of io, jo loop */ } } /* end of ll( or il, jl) loop */ mpp_delete_domain2d(&Dom); }/* for(nl=0; nl 0) { double *g_area, *g_clon, *g_clat; int *g_il, *g_jl, *g_io, *g_jo; int ii; g_il = (int *)malloc(nxgrid*sizeof(int )); g_jl = (int *)malloc(nxgrid*sizeof(int )); g_io = (int *)malloc(nxgrid*sizeof(int )); g_jo = (int *)malloc(nxgrid*sizeof(int )); g_area = (double *)malloc(nxgrid*sizeof(double)); g_clon = (double *)malloc(nxgrid*sizeof(double)); g_clat = (double *)malloc(nxgrid*sizeof(double)); mpp_gather_field_int (nlxo[nl][no], lndxocn_il[nl][no], g_il); mpp_gather_field_int (nlxo[nl][no], lndxocn_jl[nl][no], g_jl); mpp_gather_field_int (nlxo[nl][no], lndxocn_io[nl][no], g_io); mpp_gather_field_int (nlxo[nl][no], lndxocn_jo[nl][no], g_jo); mpp_gather_field_double(nlxo[nl][no], lndxocn_area[nl][no], g_area); mpp_gather_field_double(nlxo[nl][no], lndxocn_clon[nl][no], g_clon); mpp_gather_field_double(nlxo[nl][no], lndxocn_clat[nl][no], g_clat); for(i=0; i 0) { l_clon[ll] /= l_area[ll]; l_clat[ll] /= l_area[ll]; } } /* substract land centroid to get the centroid distance between land grid and exchange grid. */ for(no=0; no 0) { o_clon[no][lo] /= o_area[no][lo]; o_clat[no][lo] /= o_area[no][lo]; } } for(nl=0; nl0) { size_t start[4], nwrite[4]; int *gdata_int; double *gdata_dbl; char contact[STRING]; int fid, dim_string, dim_ncells, dim_two, dims[4]; int id_contact, id_xgrid_area, n; int id_tile1_cell, id_tile2_cell, id_tile1_dist, id_tile2_dist; for(i=0; i<4; i++) { start[i] = 0; nwrite[i] = 1; } if(same_mosaic) sprintf(lxo_file[nfile_lxo], "lnd_%s_%sXocn_%s_%s.nc", lmosaic_name, ltile_name[nl], omosaic_name, otile_name[no]); else sprintf(lxo_file[nfile_lxo], "%s_%sX%s_%s.nc", lmosaic_name, ltile_name[nl], omosaic_name, otile_name[no]); sprintf(contact, "%s:%s::%s:%s", lmosaic_name, ltile_name[nl], omosaic_name, otile_name[no]); fid = mpp_open(lxo_file[nfile_lxo], MPP_WRITE); mpp_def_global_att(fid, "grid_version", GRID_VERSION); mpp_def_global_att(fid, "code_version", TAGNAME); mpp_def_global_att(fid, "history", history); dim_string = mpp_def_dim(fid, "string", STRING); dim_ncells = mpp_def_dim(fid, "ncells", nxgrid); dim_two = mpp_def_dim(fid, "two", 2); if(interp_order == 2) { id_contact = mpp_def_var(fid, "contact", MPP_CHAR, 1, &dim_string, 7, "standard_name", "grid_contact_spec", "contact_type", "exchange", "parent1_cell", "tile1_cell", "parent2_cell", "tile2_cell", "xgrid_area_field", "xgrid_area", "distant_to_parent1_centroid", "tile1_distance", "distant_to_parent2_centroid", "tile2_distance"); } else { id_contact = mpp_def_var(fid, "contact", MPP_CHAR, 1, &dim_string, 5, "standard_name", "grid_contact_spec", "contact_type", "exchange", "parent1_cell", "tile1_cell", "parent2_cell", "tile2_cell", "xgrid_area_field", "xgrid_area" ); } dims[0] = dim_ncells; dims[1] = dim_two; id_tile1_cell = mpp_def_var(fid, "tile1_cell", MPP_INT, 2, dims, 1, "standard_name", "parent1_cell_indices"); id_tile2_cell = mpp_def_var(fid, "tile2_cell", MPP_INT, 2, dims, 1, "standard_name", "parent2_cell_indices"); id_xgrid_area = mpp_def_var(fid, "xgrid_area", MPP_DOUBLE, 1, &dim_ncells, 2, "standard_name", "exchange_grid_area", "units", "m2"); if(interp_order == 2) { id_tile1_dist = mpp_def_var(fid, "tile1_distance", MPP_DOUBLE, 2, dims, 1, "standard_name", "distance_from_parent1_cell_centroid"); id_tile2_dist = mpp_def_var(fid, "tile2_distance", MPP_DOUBLE, 2, dims, 1, "standard_name", "distance_from_parent2_cell_centroid"); } mpp_end_def(fid); /* the index will start from 1, instead of 0 ( fortran index) */ for(i = 0;i < nlxo[nl][no]; i++) { ++(lndxocn_il[nl][no][i]); ++(lndxocn_jl[nl][no][i]); ++(lndxocn_io[nl][no][i]); lndxocn_jo[nl][no][i] += 1 - ocn_south_ext; /* one artificial j-level may be added in the south end */ } nwrite[0] = strlen(contact); mpp_put_var_value_block(fid, id_contact, start, nwrite, contact); nwrite[0] = nxgrid; gdata_int = (int *)malloc(nxgrid*sizeof(int)); gdata_dbl = (double *)malloc(nxgrid*sizeof(double)); mpp_gather_field_double(nlxo[nl][no], lndxocn_area[nl][no], gdata_dbl); mpp_put_var_value(fid, id_xgrid_area, gdata_dbl); mpp_gather_field_int(nlxo[nl][no], lndxocn_il[nl][no], gdata_int); mpp_put_var_value_block(fid, id_tile1_cell, start, nwrite, gdata_int); mpp_gather_field_int(nlxo[nl][no], lndxocn_io[nl][no], gdata_int); mpp_put_var_value_block(fid, id_tile2_cell, start, nwrite, gdata_int); start[1] = 1; mpp_gather_field_int(nlxo[nl][no], lndxocn_jl[nl][no], gdata_int); mpp_put_var_value_block(fid, id_tile1_cell, start, nwrite, gdata_int); mpp_gather_field_int(nlxo[nl][no], lndxocn_jo[nl][no], gdata_int); mpp_put_var_value_block(fid, id_tile2_cell, start, nwrite, gdata_int); if(interp_order == 2) { start[1] = 0; mpp_gather_field_double(nlxo[nl][no], lndxocn_dil[nl][no], gdata_dbl); mpp_put_var_value_block(fid, id_tile1_dist, start, nwrite, gdata_dbl); mpp_gather_field_double(nlxo[nl][no], lndxocn_dio[nl][no], gdata_dbl); mpp_put_var_value_block(fid, id_tile2_dist, start, nwrite, gdata_dbl); start[1] = 1; mpp_gather_field_double(nlxo[nl][no], lndxocn_djl[nl][no], gdata_dbl); mpp_put_var_value_block(fid, id_tile1_dist, start, nwrite, gdata_dbl); mpp_gather_field_double(nlxo[nl][no], lndxocn_djo[nl][no], gdata_dbl); mpp_put_var_value_block(fid, id_tile2_dist, start, nwrite, gdata_dbl); } mpp_close(fid); free(gdata_int); free(gdata_dbl); ++nfile_lxo; } } /* for(no=0; no