/* * Modify land grid to match ocean grid at coast and calculate atmos/land, * atmos/ocean, and land/ocean overlaps using the Sutherland-Hodgeman polygon * clipping algorithm (Sutherland, I. E. and G. W. Hodgeman, 1974: Reentrant * polygon clipping, CACM, 17(1), 32-42). * Warning, when the atmos grid is cubic grid, the number of model points should be * even to avoid tiling error. I will come back to solve this issue in the future. * */ #include #include #include #include #include #include #include "constant.h" #include "mpp.h" #include "mpp_io.h" #include "tool_util.h" #include "mpp_domain.h" #include "mosaic_util.h" #include "create_xgrid.h" #define print_grid 0 char *usage[] = { "", " make_coupler_mosaic --atmos_mosaic atmos_mosaic.nc --ocean_mosaic ocean_mosaic.nc ", " --ocean_topog ocean_topog.nc [--land_mosaic land_mosaic.nc] [--wave_mosaic]", " [--sea_level #] [--interp_method #] [--mosaic_name mosaic_name] ", " [--area_ratio_thresh #] [--check ] [--verbose] [--print_memory] ", " ", "make_coupler_mosaic generates three exchange grids for the FMS coupler. The output ", "file includes exchange grid files for fluxes between atmosphere and surface (sea ice ", "and land), exchange grid files for runoff between land and sea ice. There might be more ", "than one exchange grid files between two model solo mosaic because there might be ", "multiple tiles in a solo mosaic. The exchange grid information is between model ", "grid, not between supergrid. We assume the refinement ratio between model grid and ", "supergrid is 2. Currently we only output the exchange grid on T-cell. ", "Besides generating the exchange grid files, make_coupler_mosaic also generates the ", "coupler mosaic file (the file name will be mosaic_name.nc) which contains the atmos, ", "land and ocean mosaic path, ocean mosaic topog path and exchange grid file path. ", "make_coupler_mosaic expects NetCDF format input.", " ", "make_coupler_mosaic takes the following flags:", "", "REQUIRED:", "", "--atmos_mosaic atmos_mosaic.nc specify the atmosphere mosaic information. This file", " contains list of tile files which specify the grid ", " information for each tile. Each grid is required to be ", " logically rectangular grid. The file name can not be 'mosaic.nc' ", "", "--ocean_mosaic ocean_mosaic.nc specify the ocean mosaic information. This file", " contains list of tile files which specify the grid ", " information for each tile. The file name can not be 'mosaic.nc' ", " ", "--ocean_topog ocean_topog.nc specify the topography information for ocean mosaic.", " The field name of the topography is depth_tile# or depth when ", " ntiles = 1, The topography data is positive down.", " ", "OPTIONAL FLAGS", "", "--land_mosaic land_mosaic.nc specify the land mosaic information. This file", " contains list of tile files which specify the grid ", " information for each tile. Each grid is required to be ", " logically rectangular grid. When land_mosaic is not specified,", " atmosphere mosaic will be used to specify land mosaic.", " The file name can not be 'mosaic.nc'.", " ", "--wave_mosaic wave_mosaic.nc specify the wave mosaic information. This file", " contains list of tile files which specify the grid ", " information for each tile. Each grid is required to be ", " logically rectangular grid. When wave_mosaic is specified, ", " exchange grid information between wave mosaic and ocean ", " mosaic will be generated. Otherwise none of both will be created.", " ", "--interp_order # specify the order of conservative interplation. Its value ", " can be 1 ( linear order ) or 2 ( second order ) with default ", " value 2. ", " ", "--sea_level # specify the sea level ( in meters ) and its value will be used", " to determine land/sea mask. When topography of ", " a grid cell is less than sea level, this grid cell will be land,", " otherwise it will be ocean. Default value is 0", " ", "--mosaic_name mosaic_name coupler mosaic name. The output coupler mosaic file will be ", " mosaic_name.nc. default value is 'mosaic'. ", " ", "--area_ratio_thresh # Criteria to decide if an overlap between any two model is an ", " exchange grid or not. When overlap area/model grid area is greater ", " than area_ratio_thresh, that overlap is saved as exchange grid. ", " The default value is 1.e-6 ", " ", "--check check the tiling error", "", "--print_memory debug memory usage when it is set ", " ", "--verbose Set --verbose to print out messages during running. ", "", "A sample call to make_coupler_mosaic that makes exchange grids for atmosphere, land and ocean ", "mosaic (atmosphere and land are coincident) is: ", "", " make_coupler_mosaic --atmos_mosaic atmos_mosaic.nc --ocean_mosaic ocean_mosaic.nc ", " --ocean_topog ocean_topog.nc", "", NULL }; #define D2R (M_PI/180.) #define R2D (180./M_PI) #define MAXXGRIDFILE 100 #define MX 2000 #define TINY_VALUE (1.e-7) #define TOLORENCE (1.e-4) #define MIN_AREA_FRAC (1.e-4) #define GREAT_CIRCLE_CLIP 1 #define LEGACY_CLIP 2 char grid_version[] = "0.2"; char tagname[] = "$Name: tikal $"; /* This file will get the directory that stores the file and the file name (without the dir path) */ void get_file_dir_and_name(char *file, char *filedir, char *filename) { char *fptr=NULL; int siz; fptr = strrchr(file, '/'); if(!fptr) { strcpy(filename, file); strcpy(filedir, "./"); } else { ++fptr; siz = fptr - file; strcpy(filename, fptr); strncpy(filedir, file, siz); } }; void get_global_grid(const char *grid_file, int nx, int ny, int x_refine, int y_refine, double *x, double *y) { double *x_local, *y_local, *tmp; size_t start[4], nread[4]; int nxc, nyc, isc, iec, jsc, jec; int isc2, iec2, jsc2, jec2; int nxc2, nyc2, layout[2]; int fid, vid, i, j; domain2D Dom; /* define a temporary domain to read data on local and then use mpp_global_field to get global field. This may solve the IO issue for high resolution grid */ mpp_define_layout(nx, ny, mpp_npes(), layout); mpp_define_domain2d(nx, ny, layout, 0, 0, &Dom); mpp_get_compute_domain2d(Dom, &isc, &iec, &jsc, &jec); nxc = iec - isc + 1; nyc = jec - jsc + 1; x_local = (double *)malloc((nxc+1)*(nyc+1)*sizeof(double)); y_local = (double *)malloc((nxc+1)*(nyc+1)*sizeof(double)); isc2 = x_refine*isc; iec2 = x_refine*iec + 1; jsc2 = y_refine*jsc; jec2 = y_refine*jec + 1; nxc2 = iec2 - isc2 + 1; nyc2 = jec2 - jsc2 + 1; tmp = (double *)malloc((nxc2+1)*(nyc2+1)*sizeof(double)); start[0] = jsc2; start[1] = isc2; start[2] = 0; start[3] = 0; nread[0] = nyc2+1; nread[1] = nxc2+1; nread[2] = 1; nread[3] = 1; fid = mpp_open(grid_file, MPP_READ); vid = mpp_get_varid(fid, "x"); mpp_get_var_value_block(fid, vid, start, nread, tmp); for(j=0; j<=nyc; j++) for(i=0; i<=nxc; i++) x_local[j*(nxc+1)+i] = tmp[(j*y_refine)*(nxc2+1)+i*x_refine]; mpp_global_field_all_double(Dom, nxc+1, nyc+1, x_local, x); vid = mpp_get_varid(fid, "y"); mpp_get_var_value_block(fid, vid, start, nread, tmp); for(j=0; j<=nyc; j++) for(i=0; i<=nxc; i++) y_local[j*(nxc+1)+i] = tmp[(j*y_refine)*(nxc2+1)+i*x_refine]; mpp_global_field_all_double(Dom, nxc+1, nyc+1, y_local, y); mpp_delete_domain2d(&Dom); mpp_close(fid); free(tmp); free(x_local); free(y_local); } void get_global_data(const char *data_file, const char *fieldname, int nx, int ny, double *data, int ioff, int joff) { double *data_local; size_t start[4], nread[4]; int nxc, nyc, isc, iec, jsc, jec; int fid, vid, layout[2]; domain2D Dom; mpp_define_layout(nx, ny, mpp_npes(), layout); mpp_define_domain2d(nx, ny, layout, 0, 0, &Dom); mpp_get_compute_domain2d(Dom, &isc, &iec, &jsc, &jec); nxc = iec - isc + 1; nyc = jec - jsc + 1; data_local = (double *)malloc((nxc+ioff)*(nyc+joff)*sizeof(double)); start[0] = jsc; start[1] = isc; start[2] = 0; start[3] = 0; nread[0] = nyc+joff; nread[1] = nxc+ioff; nread[2] = 1; nread[3] = 1; fid = mpp_open(data_file, MPP_READ); vid = mpp_get_varid(fid, fieldname); mpp_get_var_value_block(fid, vid, start, nread, data_local); mpp_close(fid); mpp_global_field_all_double(Dom, nxc+ioff, nyc+joff, data_local, data); free(data_local); } void get_grid_global_area(int nx, int ny, const double *x, const double *y, double *area) { double *x_local, *y_local, *area_local; int nxc, nyc, isc, iec, jsc, jec; int isc2, iec2, jsc2, jec2; int nxc2, nyc2, layout[2]; int i, j; domain2D Dom; mpp_define_layout(nx, ny, mpp_npes(), layout); mpp_define_domain2d(nx, ny, layout, 0, 0, &Dom); mpp_get_compute_domain2d(Dom, &isc, &iec, &jsc, &jec); nxc = iec - isc + 1; nyc = jec - jsc + 1; x_local = (double *)malloc((nxc+1)*(nyc+1)*sizeof(double)); y_local = (double *)malloc((nxc+1)*(nyc+1)*sizeof(double)); for(j=0; j<=nyc; j++) for(i=0; i<=nxc; i++) { x_local[j*(nxc+1)+i] = x[(j+jsc)*(nx+1)+i+isc]; y_local[j*(nxc+1)+i] = y[(j+jsc)*(nx+1)+i+isc]; } area_local = (double *)malloc(nxc*nyc*sizeof(double)); get_grid_area(&nxc, &nyc, x_local, y_local, area_local); mpp_global_field_all_double(Dom, nxc, nyc, area_local, area); free(x_local); free(y_local); free(area_local); mpp_delete_domain2d(&Dom); } int get_nest_contact(const int *nx, const int *ny, int ncontacts, const int *tile1, const int *tile2, const int *istart1, const int *iend1, const int *jstart1, const int *jend1, const int *istart2, const int *iend2, const int *jstart2, const int *jend2, int *tile_nest, int *tile_parent, int *is_nest, int *ie_nest, int *js_nest, int *je_nest, int *is_parent, int *ie_parent, int *js_parent, int *je_parent) { int n, nx1_contact, ny1_contact, nx2_contact, ny2_contact; int t1, t2; int nnest=0; for(n=0; n1)mpp_error("make_coupler_mosaic(get_nest_contact): only support one nest region, contact developer"); if(nx2_contact*ny2_contact > nx1_contact*ny1_contact) { if(nx2_contact%nx1_contact || ny2_contact%ny1_contact ) if(mpp_pe()==mpp_root_pe()) mpp_error("make_coupler_mosaic(get_nest_contact):it is not a integer refinement"); is_nest [0] = istart2[n]; ie_nest [0] = iend2 [n]; js_nest [0] = jstart2[n]; je_nest [0] = jend2 [n]; tile_nest [0] = tile2 [n]-1; is_parent [0] = istart1[n]; ie_parent [0] = iend1 [n]; js_parent [0] = jstart1[n]; je_parent [0] = jend1 [n]; tile_parent[0] = tile1 [n]-1; } else { if(nx1_contact%nx2_contact || ny1_contact%ny2_contact ) mpp_error("make_coupler_mosaic(get_nest_contact):it is not a integer refinement"); is_nest [0] = istart1[n]; ie_nest [0] = iend1 [n]; js_nest [0] = jstart1[n]; je_nest [0] = jend1 [n]; tile_nest [0] = tile1 [n]-1; is_parent [0] = istart2[n]; ie_parent [0] = iend2 [n]; js_parent [0] = jstart2[n]; je_parent [0] = jend2 [n]; tile_parent[0] = tile2 [n]-1; } } return nnest; } int main (int argc, char *argv[]) { int c, i, same_mosaic; extern char *optarg; char *omosaic = NULL; char *amosaic = NULL; char *lmosaic = NULL; char *wmosaic = NULL; char *otopog = NULL; char mosaic_name[STRING] = "mosaic", mosaic_file[STRING]; char omosaic_name[STRING], amosaic_name[STRING], lmosaic_name[STRING], wmosaic_name[STRING]; char **otile_name=NULL, **atile_name=NULL, **ltile_name=NULL, **wtile_name=NULL; int x_refine = 2, y_refine = 2; int interp_order = 2; double area_ratio_thresh = 1.0e-6; unsigned int check = 0; unsigned int verbose = 0; int errflg = (argc == 1); char history[512]; int nfile_lxo=0, nfile_axo=0, nfile_axl=0; int nfile_wxo=0; char lxo_file[MAXXGRIDFILE][STRING]; char axo_file[MAXXGRIDFILE][STRING]; char axl_file[MAXXGRIDFILE][STRING]; char wxo_file[MAXXGRIDFILE][STRING]; char amosaic_dir[STRING], amosaic_file[STRING]; char lmosaic_dir[STRING], lmosaic_file[STRING]; char omosaic_dir[STRING], omosaic_file[STRING]; char wmosaic_dir[STRING], wmosaic_file[STRING]; char otopog_dir[STRING], otopog_file[STRING]; int ntile_ocn, ntile_atm, ntile_lnd, ntile_wav; int *nxo = NULL, *nyo = NULL, *nxa = NULL, *nya = NULL, *nxl = NULL, *nyl = NULL; int *nxw = NULL, *nyw = NULL; double **xocn = NULL, **yocn = NULL, **xatm = NULL, **yatm = NULL, **xlnd = NULL, **ylnd = NULL; double **xwav = NULL, **ywav = NULL; double **cart_xatm=NULL, **cart_yatm=NULL, **cart_zatm=NULL; double **cart_xocn=NULL, **cart_yocn=NULL, **cart_zocn=NULL; double **cart_xlnd=NULL, **cart_ylnd=NULL, **cart_zlnd=NULL; double **cart_xwav=NULL, **cart_ywav=NULL, **cart_zwav=NULL; double **area_ocn = NULL, **area_lnd = NULL, **area_atm = NULL, **area_wav = NULL; double **atm_xarea=NULL; double **omask = NULL; double sea_level = 0.; int clip_method = LEGACY_CLIP; int atm_great_circle_algorithm=0; int lnd_great_circle_algorithm=0; int ocn_great_circle_algorithm=0; int lnd_same_as_atm = 0; int ocn_same_as_atm = 0; int wav_same_as_ocn = 0; int option_index = 0; double axo_area_sum = 0, axl_area_sum = 0; double axo_area_sum_nest = 0, axl_area_sum_nest = 0; double wxo_area_sum = 0; int ocn_south_ext = 0; int tile_nest, is_nest, ie_nest, js_nest, je_nest; int tile_parent, is_parent, ie_parent, js_parent, je_parent; int print_memory=0; static struct option long_options[] = { {"atmos_mosaic", required_argument, NULL, 'a'}, {"land_mosaic", required_argument, NULL, 'l'}, {"ocean_mosaic", required_argument, NULL, 'o'}, {"wave_mosaic", required_argument, NULL, 'w'}, {"ocean_topog", required_argument, NULL, 't'}, {"sea_level", required_argument, NULL, 's'}, {"interp_order", required_argument, NULL, 'i'}, {"mosaic_name", required_argument, NULL, 'm'}, {"area_ratio_thresh", required_argument, NULL, 'r'}, {"check", no_argument, NULL, 'n'}, {"verbose", no_argument, NULL, 'v'}, {"print_memory", no_argument, NULL, 'p'}, {NULL, 0, NULL, 0} }; mpp_init(&argc, &argv); mpp_domain_init(); /* * process command line */ while ((c = getopt_long(argc, argv, "i:", long_options, &option_index) ) != -1) switch (c) { case 'a': amosaic = optarg; break; case 'l': lmosaic = optarg; break; case 'o': omosaic = optarg; break; case 'w': wmosaic = optarg; break; case 't': otopog = optarg; break; case 'i': interp_order = atoi(optarg); break; case 's': sea_level = atof(optarg); break; case 'm': strcpy(mosaic_name,optarg); break; case 'r': area_ratio_thresh = atof(optarg); case 'n': check = 1; break; case 'v': verbose = 1; break; case 'p': print_memory = 1; break; case '?': errflg++; } if (errflg || !amosaic || !omosaic || !otopog) { char **u = usage; while (*u) { fprintf(stderr, "%s\n", *u); u++; } exit(2); } /* interp_order should be 1 or 2 */ if(interp_order != 1 && interp_order !=2 )mpp_error("make_coupler_mosaic: interp_order should be 1 or 2"); strcpy(history,argv[0]); for(i=1;i0) { if( atm_great_circle_algorithm != great_circle_algorithm) mpp_error("make_topog: atribute 'great_circle_algorithm' of field 'tile' have different value for different tile"); } atm_great_circle_algorithm = great_circle_algorithm; } mpp_close(g_fid); if(nxa[n]%x_refine != 0 ) mpp_error("make_coupler_mosaic: atmos supergrid x-size can not be divided by x_refine"); if(nya[n]%y_refine != 0 ) mpp_error("make_coupler_mosaic: atmos supergrid y-size can not be divided by y_refine"); nxa[n] /= x_refine; nya[n] /= y_refine; xatm[n] = (double *)malloc((nxa[n]+1)*(nya[n]+1)*sizeof(double)); yatm[n] = (double *)malloc((nxa[n]+1)*(nya[n]+1)*sizeof(double)); area_atm[n] = (double *)malloc((nxa[n] )*(nya[n] )*sizeof(double)); if(check) { atm_xarea[n]= (double *)malloc((nxa[n] )*(nya[n] )*sizeof(double)); for(i=0; i0) { tile1 = (int *)malloc(ncontacts*sizeof(int)); tile2 = (int *)malloc(ncontacts*sizeof(int)); istart1 = (int *)malloc(ncontacts*sizeof(int)); iend1 = (int *)malloc(ncontacts*sizeof(int)); jstart1 = (int *)malloc(ncontacts*sizeof(int)); jend1 = (int *)malloc(ncontacts*sizeof(int)); istart2 = (int *)malloc(ncontacts*sizeof(int)); iend2 = (int *)malloc(ncontacts*sizeof(int)); jstart2 = (int *)malloc(ncontacts*sizeof(int)); jend2 = (int *)malloc(ncontacts*sizeof(int)); read_mosaic_contact(amosaic, tile1, tile2, istart1, iend1, jstart1, jend1, istart2, iend2, jstart2, jend2); nnest = get_nest_contact( nxa, nya, ncontacts, tile1, tile2, istart1, iend1, jstart1, jend1, istart2, iend2, jstart2, jend2, &tile_nest, &tile_parent, &is_nest, &ie_nest, &js_nest, &je_nest, &is_parent, &ie_parent, &js_parent, &je_parent); free(tile1); free(tile2); free(istart1); free(iend1); free(jstart1); free(jend1); free(istart2); free(iend2); free(jstart2); free(jend2); } } if(print_memory)print_mem_usage("after read atmosphere grid"); /* * Read land grid */ if (strcmp(lmosaic, amosaic) ) { /* land mosaic is different from atmosphere mosaic */ int n, m_fid, g_fid, vid, gid, tid; size_t start[4], nread[4]; char dir[STRING], filename[STRING], file[2*STRING]; for(n=0; n<4; n++) { start[n] = 0; nread[n] = 1; } m_fid = mpp_open(lmosaic, MPP_READ); vid = mpp_get_varid(m_fid, "mosaic"); mpp_get_var_value(m_fid, vid, lmosaic_name); ntile_lnd = mpp_get_dimlen(m_fid, "ntiles"); nxl = (int *) malloc (ntile_lnd*sizeof(int)); nyl = (int *) malloc (ntile_lnd*sizeof(int)); xlnd = (double **) malloc( ntile_lnd*sizeof(double *)); ylnd = (double **) malloc( ntile_lnd*sizeof(double *)); area_lnd = (double **) malloc( ntile_lnd*sizeof(double *)); ltile_name = (char **)malloc(ntile_lnd*sizeof(char *)); /* grid should be located in the same directory of mosaic file */ get_file_path(lmosaic, dir); gid = mpp_get_varid(m_fid, "gridfiles"); tid = mpp_get_varid(m_fid, "gridtiles"); for(n=0; n0) { if( lnd_great_circle_algorithm != great_circle_algorithm) mpp_error("make_topog: atribute 'great_circle_algorithm' of field 'tile' have different value for different tile"); } lnd_great_circle_algorithm = great_circle_algorithm; } mpp_close(g_fid); if(nxl[n]%x_refine != 0 ) mpp_error("make_coupler_mosaic: land supergrid x-size can not be divided by x_refine"); if(nyl[n]%y_refine != 0 ) mpp_error("make_coupler_mosaic: land supergrid y-size can not be divided by y_refine"); nxl[n] /= x_refine; nyl[n] /= y_refine; xlnd[n] = (double *)malloc((nxl[n]+1)*(nyl[n]+1)*sizeof(double)); ylnd[n] = (double *)malloc((nxl[n]+1)*(nyl[n]+1)*sizeof(double)); area_lnd[n] = (double *)malloc((nxl[n] )*(nyl[n] )*sizeof(double)); get_global_grid(file, nxl[n], nyl[n], x_refine, y_refine, xlnd[n], ylnd[n]); /*scale grid from degree to radian, because create_xgrid assume the grid is in radians */ for(i=0; i<(nxl[n]+1)*(nyl[n]+1); i++) { xlnd[n][i] *= D2R; ylnd[n][i] *= D2R; } } /* compute lnd_area */ if(clip_method == GREAT_CIRCLE_CLIP) { cart_xlnd = (double **) malloc( ntile_lnd*sizeof(double *)); cart_ylnd = (double **) malloc( ntile_lnd*sizeof(double *)); cart_zlnd = (double **) malloc( ntile_lnd*sizeof(double *)); for(n=0; n=0)mpp_error("make_coupler_mosaic: land mosaic must be different from atmos mosaic " "when there is nest region in atmosphere mosaic"); ntile_lnd = ntile_atm; nxl = nxa; nyl = nya; if(clip_method == GREAT_CIRCLE_CLIP) { cart_xlnd = cart_xatm; cart_ylnd = cart_yatm; cart_zlnd = cart_zatm; } else { xlnd = xatm; ylnd = yatm; } area_lnd = area_atm; lnd_same_as_atm = 1; strcpy(lmosaic_name, amosaic_name); ltile_name = atile_name; lnd_great_circle_algorithm = atm_great_circle_algorithm; } if(print_memory)print_mem_usage("after read land grid"); if (strcmp(omosaic, amosaic) == 0 ) ocn_same_as_atm = 1; /* * Read ocean grid boundaries and mask (where water is) for each tile within the mosaic. */ { int n, ntiles, m_fid, g_fid, t_fid, vid, gid, tid; size_t start[4], nread[4]; char dir[STRING], filename[STRING], file[2*STRING]; for(n=0; n<4; n++) { start[n] = 0; nread[n] = 1; } m_fid = mpp_open(omosaic, MPP_READ); vid = mpp_get_varid(m_fid, "mosaic"); mpp_get_var_value(m_fid, vid, omosaic_name); ntile_ocn = mpp_get_dimlen(m_fid, "ntiles"); nxo = (int *) malloc(ntile_ocn*sizeof(int)); nyo = (int *) malloc(ntile_ocn*sizeof(int)); xocn = (double **) malloc(ntile_ocn*sizeof(double *)); yocn = (double **) malloc(ntile_ocn*sizeof(double *)); area_ocn = (double **) malloc(ntile_ocn*sizeof(double *)); otile_name = (char **) malloc(ntile_ocn*sizeof(char *)); /* grid should be located in the same directory of mosaic file */ get_file_path(omosaic, dir); gid = mpp_get_varid(m_fid, "gridfiles"); tid = mpp_get_varid(m_fid, "gridtiles"); /* For the purpose of reproducing between processor count, the layout is set to (1, npes). */ for(n=0; n0) { if( ocn_great_circle_algorithm != great_circle_algorithm) mpp_error("make_topog: atribute 'great_circle_algorithm' of field 'tile' have different value for different tile"); } ocn_great_circle_algorithm = great_circle_algorithm; } mpp_close(g_fid); if(nxo[n]%x_refine != 0 ) mpp_error("make_coupler_mosaic: ocean supergrid x-size can not be divided by x_refine"); if(nyo[n]%y_refine != 0 ) mpp_error("make_coupler_mosaic: ocean supergrid y-size can not be divided by y_refine"); nxo[n] /= x_refine; nyo[n] /= y_refine; tmpx = (double *)malloc((nxo[n]+1)*(nyo[n]+1)*sizeof(double)); tmpy = (double *)malloc((nxo[n]+1)*(nyo[n]+1)*sizeof(double)); get_global_grid(file, nxo[n], nyo[n], x_refine, y_refine, tmpx, tmpy); /* sometimes the ocean is only covered part of atmosphere, especially not cover the south pole region. In order to get all the exchange grid between atmosXland, we need to extend one point to cover the whole atmosphere. This need the assumption of one-tile ocean. Also we assume the latitude is the along j=0 */ if(ntile_ocn == 1) { int na; int is_uniform; /* check if the latitude is uniform or not at j=1 */ is_uniform = 1; for(i=1; i<=nxo[n]; i++) { if(tmpy[i] != tmpy[i-1]) { is_uniform = 0; if(mpp_pe()==mpp_root_pe()) printf("\nNOTE from make_coupler_mosaic: ocean grid latitude is not uniform along j = 1\n"); } } /* if latitude is uniform along j = 1, may add row in southmost */ if( is_uniform ) { /* calculate the minimum of latitude of atmosphere grid */ min_atm_lat = 9999; /* dummy large value */ for(na=0; na min_lat) min_atm_lat = min_lat; } min_atm_lat = -90.*D2R; if(tmpy[0]*D2R > min_atm_lat + TINY_VALUE) { /* extend one point in south direction*/ ocn_south_ext = 1; if(mpp_pe()==mpp_root_pe())printf("make_coupler_mosaic: one row is add to the south end to cover the globe\n"); } } } nyo_old = nyo[n]; 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)); /* assign the global latitude */ for(j = 0; j < nyo_old+1; j++) for(i = 0; i < nxo[n]+1; i++) { xocn[n][(j+ocn_south_ext)*(nxo[n]+1)+i] = tmpx[j*(nxo[n]+1)+i] * D2R; yocn[n][(j+ocn_south_ext)*(nxo[n]+1)+i] = tmpy[j*(nxo[n]+1)+i] * D2R; } if(ocn_south_ext==1) { for(i=0; isea_level) omask[n][(j+ocn_south_ext)*nx+i] = 1; } free(depth); } } } if(print_memory)print_mem_usage("after read ocean grid"); /* when atm_great_circle_algorithm is 0, lnd_great_circle_algorithm/ocn_great_circle_algorithm must be 0 */ if( atm_great_circle_algorithm == 0 && (lnd_great_circle_algorithm || ocn_great_circle_algorithm)) mpp_error("make_coupler_mosaic: when atm does not use great_circle_algorithm, lnd/ocn can not use great_circle_algorithm"); /* * Read wave grid */ if ( wmosaic ) { int n, m_fid, g_fid, vid, gid, tid; size_t start[4], nread[4]; char dir[STRING], filename[STRING], file[2*STRING]; if (strcmp(wmosaic, omosaic) == 0 ) wav_same_as_ocn = 1; for(n=0; n<4; n++) { start[n] = 0; nread[n] = 1; } m_fid = mpp_open(wmosaic, MPP_READ); vid = mpp_get_varid(m_fid, "mosaic"); mpp_get_var_value(m_fid, vid, wmosaic_name); ntile_wav = mpp_get_dimlen(m_fid, "ntiles"); nxw = (int *) malloc (ntile_wav*sizeof(int)); nyw = (int *) malloc (ntile_wav*sizeof(int)); xwav = (double **) malloc( ntile_wav*sizeof(double *)); ywav = (double **) malloc( ntile_wav*sizeof(double *)); area_wav = (double **) malloc( ntile_wav*sizeof(double *)); wtile_name = (char **)malloc(ntile_wav*sizeof(char *)); /* grid should be located in the same directory of mosaic file */ get_file_path(wmosaic, dir); gid = mpp_get_varid(m_fid, "gridfiles"); tid = mpp_get_varid(m_fid, "gridtiles"); for(n=0; n ya_max) ya_max = ya[0]; if(ya[1] > ya_max) ya_max = ya[1]; if(ya[2] > ya_max) ya_max = ya[2]; if(ya[3] > ya_max) ya_max = ya[3]; if(ya[0] < ya_min) ya_min = ya[0]; if(ya[1] < ya_min) ya_min = ya[1]; if(ya[2] < ya_min) ya_min = ya[2]; if(ya[3] < ya_min) ya_min = ya[3]; } for(nl=0; nl ya_min ) { if(jl < js_lnd[nl] ) js_lnd[nl] = jl; } if( yy < ya_max ) { if(jl > je_lnd[nl] ) je_lnd[nl] = jl; } } js_lnd[nl] = max(0, js_lnd[nl]-1); je_lnd[nl] = min(nyl[nl]-1, je_lnd[nl]+1); if(nl==na || !lnd_same_as_atm ) { is_lnd[nl] = 0; ie_lnd[nl] = nxl[nl]-1; } else { is_lnd[nl] = nxl[nl]-1; ie_lnd[nl] = 0; } } for(no=0; no ya_min ) { if(jo < js_ocn[no] ) js_ocn[no] = jo; } if( yy < ya_max ) { if(jo > je_ocn[no] ) je_ocn[no] = jo; } } js_ocn[no] = max(0, js_ocn[no]-1); je_ocn[no] = min(nyo[no]-1, je_ocn[no]+1); if(no==na || !ocn_same_as_atm ) { is_ocn[no] = 0; ie_ocn[no] = nxo[no]-1; } else { is_ocn[no] = nxo[no]-1; ie_ocn[no] = 0; } } } if(mpp_pe()==mpp_root_pe() && verbose)printf("na = %d, la = %d, is=%d, ie = %d\n", na, la, is, ie); for(la=is;la<=ie;la++) { ia = la%nxa[na]; ja = la/nxa[na]; if(print_grid) { n0 = ja *(nxa[na]+1) + ia; n1 = ja *(nxa[na]+1) + ia+1; n2 = (ja+1)*(nxa[na]+1) + ia+1; n3 = (ja+1)*(nxa[na]+1) + ia; xa[0] = xatm[na][n0]; ya[0] = yatm[na][n0]; xa[1] = xatm[na][n1]; ya[1] = yatm[na][n1]; xa[2] = xatm[na][n2]; ya[2] = yatm[na][n2]; xa[3] = xatm[na][n3]; ya[3] = yatm[na][n3]; printf("atm grid is \n"); printf("%15.11f, %15.11f \n", xa[0]*R2D, ya[0]*R2D); printf("%15.11f, %15.11f \n", xa[1]*R2D, ya[1]*R2D); printf("%15.11f, %15.11f \n", xa[2]*R2D, ya[2]*R2D); printf("%15.11f, %15.11f \n", xa[3]*R2D, ya[3]*R2D); printf("%15.11f, %15.11f \n", xa[0]*R2D, ya[0]*R2D); } if(clip_method == GREAT_CIRCLE_CLIP) { /*clockwise*/ n0 = ja *(nxa[na]+1) + ia; n1 = (ja+1)*(nxa[na]+1) + ia; n2 = (ja+1)*(nxa[na]+1) + ia+1; n3 = ja *(nxa[na]+1) + ia+1; xa[0] = cart_xatm[na][n0]; ya[0] = cart_yatm[na][n0]; za[0] = cart_zatm[na][n0]; xa[1] = cart_xatm[na][n1]; ya[1] = cart_yatm[na][n1]; za[1] = cart_zatm[na][n1]; xa[2] = cart_xatm[na][n2]; ya[2] = cart_yatm[na][n2]; za[2] = cart_zatm[na][n2]; xa[3] = cart_xatm[na][n3]; ya[3] = cart_yatm[na][n3]; za[3] = cart_zatm[na][n3]; } else { n0 = ja *(nxa[na]+1) + ia; n1 = ja *(nxa[na]+1) + ia+1; n2 = (ja+1)*(nxa[na]+1) + ia+1; n3 = (ja+1)*(nxa[na]+1) + ia; xa[0] = xatm[na][n0]; ya[0] = yatm[na][n0]; xa[1] = xatm[na][n1]; ya[1] = yatm[na][n1]; xa[2] = xatm[na][n2]; ya[2] = yatm[na][n2]; xa[3] = xatm[na][n3]; ya[3] = yatm[na][n3]; ya_min = minval_double(4, ya); ya_max = maxval_double(4, ya); na_in = fix_lon(xa, ya, 4, M_PI); xa_min = minval_double(na_in, xa); xa_max = maxval_double(na_in, xa); xa_avg = avgval_double(na_in, xa); } count = 0; for(nl=0; nl= 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; n_out = clip_2dx2d( xa, ya, na_in, xl, yl, nl_in, x_out, y_out ); } if ( n_out > 0 ) { if(clip_method == GREAT_CIRCLE_CLIP) xarea=great_circle_area ( n_out, x_out, y_out, z_out); else 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 ) { if(print_grid) { double xtmp[20],ytmp[20]; printf("n_axl is %d\n", n_out); /* convert to lon-lat */ xyz2latlon(n_out, x_out, y_out, z_out, xtmp, ytmp); 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 MIN_AREA_FRAC) { /* over sea/ice */ /* xo should in the same range as xa after lon_fix, so no need to consider cyclic condition */ if( clip_method == GREAT_CIRCLE_CLIP ) { n_out = clip_2dx2d_great_circle(xa, ya, za, 4, xo, yo, zo, 4, x_out, y_out, z_out); } else { if(xa_min >= xo_max || xa_max <= xo_min || yo_min >= ya_max || yo_max <= ya_min ) continue; n_out = clip_2dx2d( xa, ya, na_in, xo, yo, no_in, x_out, y_out ); } if ( n_out > 0) { if( clip_method == GREAT_CIRCLE_CLIP ) xarea=great_circle_area ( n_out, x_out, y_out, z_out)*ocn_frac; else xarea = poly_area(x_out, y_out, n_out )*ocn_frac; 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)*ocn_frac; atmxocn_clat[na][no][naxo[na][no]] = poly_ctrlat ( x_out, y_out, n_out )*ocn_frac; } ++(naxo[na][no]); if(naxo[na][no] > MAXXGRID) mpp_error("naxo is greater than MAXXGRID, increase MAXXGRID"); } } } if(lnd_frac > MIN_AREA_FRAC) { /* 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; n_out = clip_2dx2d( atmxlnd_x[l], atmxlnd_y[l], num_v[l], xo, yo, no_in, x_out, y_out ); } if( n_out > 0) { if( clip_method == GREAT_CIRCLE_CLIP ) xarea=great_circle_area ( n_out, x_out, y_out, z_out)*lnd_frac; else xarea = poly_area(x_out, y_out, n_out )*lnd_frac; 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) { if(print_grid) { double xtmp[20],ytmp[20]; printf("num exchange grid between ocean and axl is %d\n", n_out); /* convert to lon-lat */ xyz2latlon(n_out, x_out, y_out, z_out, xtmp, ytmp); for(n=0; n 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); free(is_lnd); free(ie_lnd); free(js_lnd); free(je_lnd); free(is_ocn); free(ie_ocn); free(js_ocn); free(je_ocn); if(print_memory) { sprintf(mesg, "end of loop na=%d", na); print_mem_usage(mesg); } } /* end of na loop */ if(print_memory)print_mem_usage("after calcuting exchange grid"); time_end = time(NULL); if(verbose) printf("one pe %d, The loop used %f seconds.\n", mpp_pe(), difftime(time_end, time_start)); /* 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 ) { nbad++; 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); } 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); if( clip_method == GREAT_CIRCLE_CLIP) mpp_def_global_att(fid, "great_circle_algorithm", "TRUE"); 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); } if(nbad>0) { printf("make_coupler_mosaic: number of points with omask != ofrac is %d\n", nbad); mpp_error("make_coupler_mosaic: omask is not equal ocn_frac"); } } /* 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); if( clip_method == GREAT_CIRCLE_CLIP) mpp_def_global_att(fid, "great_circle_algorithm", "TRUE"); 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); if( clip_method == GREAT_CIRCLE_CLIP) mpp_def_global_att(fid, "great_circle_algorithm", "TRUE"); 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) { int *gdata_ia=NULL, *gdata_ja=NULL; int ia, ja; gdata_ia = (int *)malloc(nxgrid*sizeof(int)); gdata_ja = (int *)malloc(nxgrid*sizeof(int)); mpp_gather_field_int(naxl[na][nl], atmxlnd_ia[na][nl], gdata_ia); mpp_gather_field_int(naxl[na][nl], atmxlnd_ja[na][nl], gdata_ja); 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); if( clip_method == GREAT_CIRCLE_CLIP) mpp_def_global_att(fid, "great_circle_algorithm", "TRUE"); 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) { int *gdata_ia=NULL, *gdata_ja=NULL; int ia, ja; gdata_ia = (int *)malloc(nxgrid*sizeof(int)); gdata_ja = (int *)malloc(nxgrid*sizeof(int)); mpp_gather_field_int(naxo[na][no], atmxocn_ia[na][no], gdata_ia); mpp_gather_field_int(naxo[na][no], atmxocn_ja[na][no], gdata_ja); for(n=0; n yl_max) yl_max = yl[0]; if(yl[1] > yl_max) yl_max = yl[1]; if(yl[2] > yl_max) yl_max = yl[2]; if(yl[3] > yl_max) yl_max = yl[3]; if(yl[0] < yl_min) yl_min = yl[0]; if(yl[1] < yl_min) yl_min = yl[1]; if(yl[2] < yl_min) yl_min = yl[2]; if(yl[3] < yl_min) yl_min = yl[3]; } for(no=0; no yl_min && yy < yl_max ) { if(jo > je_ocn[no] ) je_ocn[no] = jo; if(jo < js_ocn[no] ) js_ocn[no] = jo; } } js_ocn[no] = max(0, js_ocn[no]-1); je_ocn[no] = min(nyo[no]-1, je_ocn[no]+1); is_ocn[no] = 0; ie_ocn[no] = nxo[no] - 1; } } for(ll=is;ll<=ie;ll++) { il = ll%nxl[nl]; jl = ll/nxl[nl]; if(clip_method == GREAT_CIRCLE_CLIP) { /* clockwise */ n0 = jl *(nxl[nl]+1) + il; n1 = (jl+1)*(nxl[nl]+1) + il; n2 = (jl+1)*(nxl[nl]+1) + il+1; n3 = jl *(nxl[nl]+1) + il+1; xl[0] = cart_xlnd[nl][n0]; yl[0] = cart_ylnd[nl][n0]; zl[0] = cart_zlnd[nl][n0]; xl[1] = cart_xlnd[nl][n1]; yl[1] = cart_ylnd[nl][n1]; zl[1] = cart_zlnd[nl][n1]; xl[2] = cart_xlnd[nl][n2]; yl[2] = cart_ylnd[nl][n2]; zl[2] = cart_zlnd[nl][n2]; xl[3] = cart_xlnd[nl][n3]; yl[3] = cart_ylnd[nl][n3]; zl[3] = cart_zlnd[nl][n3]; } else { n0 = jl *(nxl[nl]+1) + il; n1 = jl *(nxl[nl]+1) + il+1; n2 = (jl+1)*(nxl[nl]+1) + il+1; n3 = (jl+1)*(nxl[nl]+1) + il; xl[0] = xlnd[nl][n0]; yl[0] = ylnd[nl][n0]; xl[1] = xlnd[nl][n1]; yl[1] = ylnd[nl][n1]; xl[2] = xlnd[nl][n2]; yl[2] = ylnd[nl][n2]; xl[3] = xlnd[nl][n3]; yl[3] = ylnd[nl][n3]; yl_min = minval_double(4, yl); yl_max = maxval_double(4, yl); nl_in = fix_lon(xl, yl, 4, M_PI); xl_min = minval_double(nl_in, xl); xl_max = maxval_double(nl_in, xl); xl_avg = avgval_double(nl_in, xl); } for(no=0; no MIN_AREA_FRAC) { double ocn_frac; if(clip_method == GREAT_CIRCLE_CLIP) { /* clockwise */ n0 = jo *(nxo[no]+1) + io; n1 = (jo+1)*(nxo[no]+1) + io; n2 = (jo+1)*(nxo[no]+1) + io+1; n3 = jo *(nxo[no]+1) + io+1; xo[0] = cart_xocn[no][n0]; yo[0] = cart_yocn[no][n0]; zo[0] = cart_zocn[no][n0]; xo[1] = cart_xocn[no][n1]; yo[1] = cart_yocn[no][n1]; zo[1] = cart_zocn[no][n1]; xo[2] = cart_xocn[no][n2]; yo[2] = cart_yocn[no][n2]; zo[2] = cart_zocn[no][n2]; xo[3] = cart_xocn[no][n3]; yo[3] = cart_yocn[no][n3]; zo[3] = cart_zocn[no][n3]; } else { 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; } ocn_frac = omask[no][jo*nxo[no]+io]; if(clip_method == GREAT_CIRCLE_CLIP) { n_out = clip_2dx2d_great_circle(xl, yl, zl, 4, xo, yo, zo, 4, x_out, y_out, z_out); } else n_out = clip_2dx2d( xl, yl, nl_in, xo, yo, no_in, x_out, y_out ); if ( n_out > 0 ){ if(clip_method == GREAT_CIRCLE_CLIP) xarea=great_circle_area ( n_out, x_out, y_out, z_out)*ocn_frac; else xarea = poly_area(x_out, y_out, n_out )*ocn_frac; 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)*ocn_frac; lndxocn_clat[nl][no][nlxo[nl][no]] = poly_ctrlat ( x_out, y_out, n_out )*ocn_frac; } ++(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); free(js_ocn); free(je_ocn); free(is_ocn); free(ie_ocn); }/* 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); if( clip_method == GREAT_CIRCLE_CLIP) mpp_def_global_att(fid, "great_circle_algorithm", "TRUE"); 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 yw_max) yw_max = yw[0]; if(yw[1] > yw_max) yw_max = yw[1]; if(yw[2] > yw_max) yw_max = yw[2]; if(yw[3] > yw_max) yw_max = yw[3]; if(yw[0] < yw_min) yw_min = yw[0]; if(yw[1] < yw_min) yw_min = yw[1]; if(yw[2] < yw_min) yw_min = yw[2]; if(yw[3] < yw_min) yw_min = yw[3]; } /* printf("yw_min=%f, yw_max=%f\n", yw_min, yw_max); */ for(no=0; no yw_min && yy < yw_max ) { if(jo > je_ocn[no] ) je_ocn[no] = jo; if(jo < js_ocn[no] ) js_ocn[no] = jo; } } js_ocn[no] = max(0, js_ocn[no]-1); je_ocn[no] = min(nyo[no]-1, je_ocn[no]+1); if(no==nw || !wav_same_as_ocn ) { is_ocn[no] = 0; ie_ocn[no] = nxo[no]-1; } else { is_ocn[no] = nxo[no]-1; ie_ocn[no] = 0; } } } for(lw=is;lw<=ie;lw++) { if(mpp_pe()==mpp_root_pe() && verbose)printf("nw = %d, lw = %d, is=%d, ie = %d\n", nw, lw, is, ie); iw = lw%nxw[nw]; jw = lw/nxw[nw]; if(clip_method == GREAT_CIRCLE_CLIP) { /* clockwise */ n0 = jw *(nxw[nw]+1) + iw; n1 = (jw+1)*(nxw[nw]+1) + iw; n2 = (jw+1)*(nxw[nw]+1) + iw+1; n3 = jw *(nxw[nw]+1) + iw+1; xw[0] = cart_xwav[nw][n0]; yw[0] = cart_ywav[nw][n0]; zw[0] = cart_zwav[nw][n0]; xw[1] = cart_xwav[nw][n1]; yw[1] = cart_ywav[nw][n1]; zw[1] = cart_zwav[nw][n1]; xw[2] = cart_xwav[nw][n2]; yw[2] = cart_ywav[nw][n2]; zw[2] = cart_zwav[nw][n2]; xw[3] = cart_xwav[nw][n3]; yw[3] = cart_ywav[nw][n3]; zw[3] = cart_zwav[nw][n3]; } else { n0 = jw *(nxw[nw]+1) + iw; n1 = jw *(nxw[nw]+1) + iw+1; n2 = (jw+1)*(nxw[nw]+1) + iw+1; n3 = (jw+1)*(nxw[nw]+1) + iw; xw[0] = xwav[nw][n0]; yw[0] = ywav[nw][n0]; xw[1] = xwav[nw][n1]; yw[1] = ywav[nw][n1]; xw[2] = xwav[nw][n2]; yw[2] = ywav[nw][n2]; xw[3] = xwav[nw][n3]; yw[3] = ywav[nw][n3]; yw_min = minval_double(4, yw); yw_max = maxval_double(4, yw); nw_in = fix_lon(xw, yw, 4, M_PI); xw_min = minval_double(nw_in, xw); xw_max = maxval_double(nw_in, xw); xw_avg = avgval_double(nw_in, xw); } /* calculate wave/ocean x-cells */ for(no=0; no MIN_AREA_FRAC) { /* over sea/ice */ /* xo should in the same range as xa after lon_fix, so no need to consider cyclic condition */ if(clip_method == GREAT_CIRCLE_CLIP) { n_out = clip_2dx2d_great_circle(xw, yw, zw, 4, xo, yo, zo, 4, x_out, y_out, z_out); } else { if(xw_min >= xo_max || xw_max <= xo_min || yw_min >= yo_max || yw_max <= yo_min ) continue; n_out = clip_2dx2d( xw, yw, nw_in, xo, yo, no_in, x_out, y_out); } if ( n_out > 0) { xarea = poly_area(x_out, y_out, n_out )*ocn_frac; min_area = min(area_ocn[no][jo*nxo[no]+io], area_wav[nw][lw]); if(xarea/min_area > area_ratio_thresh) { wavxocn_area[nw][no][nwxo[nw][no]] = xarea; wavxocn_io[nw][no][nwxo[nw][no]] = io; wavxocn_jo[nw][no][nwxo[nw][no]] = jo; wavxocn_iw[nw][no][nwxo[nw][no]] = iw; wavxocn_jw[nw][no][nwxo[nw][no]] = jw; if(interp_order == 2) { wavxocn_clon[nw][no][nwxo[nw][no]] = poly_ctrlon ( x_out, y_out, n_out, xw_avg)*ocn_frac; wavxocn_clat[nw][no][nwxo[nw][no]] = poly_ctrlat ( x_out, y_out, n_out )*ocn_frac; } ++(nwxo[nw][no]); if(nwxo[nw][no] > MAXXGRID) mpp_error("nwxo is greater than MAXXGRID, increase MAXXGRID"); } } } } } }/* end of la loop */ mpp_delete_domain2d(&Dom); free(is_ocn); free(ie_ocn); free(js_ocn); free(je_ocn); } /* end of nw loop */ time_end = time(NULL); if(verbose) printf("on pe %d, calculating waveXocean used %f seconds.\n", mpp_pe(), difftime(time_end, time_start)); /* calculate the centroid of model grid, as well as land_mask and ocean_mask */ { double **w_area, **w_area2; int nw, lw; w_area = (double **)malloc(ntile_wav*sizeof(double *)); for(nw =0; nw 0) { double *g_area; int *g_iw, *g_jw; int ii; g_iw = (int *)malloc(nxgrid*sizeof(int )); g_jw = (int *)malloc(nxgrid*sizeof(int )); g_area = (double *)malloc(nxgrid*sizeof(double)); mpp_gather_field_int (nwxo[nw][no], wavxocn_io[nw][no], g_iw); mpp_gather_field_int (nwxo[nw][no], wavxocn_jo[nw][no], g_jw); mpp_gather_field_double(nwxo[nw][no], wavxocn_area[nw][no], g_area); for(i=0; i 0) { double *g_area, *g_clon, *g_clat; int *g_iw, *g_jw, *g_io, *g_jo; int ii; g_iw = (int *)malloc(nxgrid*sizeof(int )); g_jw = (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 (nwxo[nw][no], wavxocn_iw[nw][no], g_iw); mpp_gather_field_int (nwxo[nw][no], wavxocn_jw[nw][no], g_jw); mpp_gather_field_int (nwxo[nw][no], wavxocn_io[nw][no], g_io); mpp_gather_field_int (nwxo[nw][no], wavxocn_jo[nw][no], g_jo); mpp_gather_field_double(nwxo[nw][no], wavxocn_area[nw][no], g_area); mpp_gather_field_double(nwxo[nw][no], wavxocn_clon[nw][no], g_clon); mpp_gather_field_double(nwxo[nw][no], wavxocn_clat[nw][no], g_clat); for(i=0; i 0) { w_clon[lw] /= w_area[nw][lw]; w_clat[lw] /= w_area[nw][lw]; } } for(no=0; no 0) { o_clon[no][lo] /= o_area[no][lo]; o_clat[no][lo] /= o_area[no][lo]; } } for(nw=0; nw 1) sprintf(wav_mask_file, "wave_mask_tile%d.nc", nw+1); else strcpy(wav_mask_file, "wave_mask.nc"); fid = mpp_open(wav_mask_file, MPP_WRITE); mpp_def_global_att(fid, "grid_version", grid_version); mpp_def_global_att(fid, "code_version", tagname); if( clip_method == GREAT_CIRCLE_CLIP) mpp_def_global_att(fid, "great_circle_algorithm", "TRUE"); mpp_def_global_att(fid, "history", history); dims[1] = mpp_def_dim(fid, "nx", nxw[nw]); dims[0] = mpp_def_dim(fid, "ny", nyw[nw]); id_mask = mpp_def_var(fid, "mask", MPP_DOUBLE, 2, dims, 2, "standard_name", "land/sea fraction at T-cell centers", "units", "none"); mpp_end_def(fid); mpp_put_var_value(fid, id_mask, mask); mpp_close(fid); free(mask); free(mask2); } } } /* write out the wavxocn */ for(nw = 0; nw < ntile_wav; nw++) { for(no = 0; no < ntile_ocn; no++) { int nxgrid; /* get total number of exchange grid on all the pes */ nxgrid = nwxo[nw][no]; mpp_sum_int(1, &nxgrid); if(nxgrid >0) { 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(wxo_file[nfile_wxo], "wav_%s_%sXocn_%s_%s.nc", wmosaic_name, wtile_name[nw], omosaic_name, otile_name[no]); else sprintf(wxo_file[nfile_wxo], "%s_%sX%s_%s.nc", wmosaic_name, wtile_name[nw], omosaic_name, otile_name[no]); sprintf(contact, "%s:%s::%s:%s", wmosaic_name, wtile_name[nw], omosaic_name, otile_name[no]); fid = mpp_open(wxo_file[nfile_wxo], MPP_WRITE); mpp_def_global_att(fid, "grid_version", grid_version); mpp_def_global_att(fid, "code_version", tagname); if( clip_method == GREAT_CIRCLE_CLIP) mpp_def_global_att(fid, "great_circle_algorithm", "TRUE"); 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 < nwxo[nw][no]; i++) { ++(wavxocn_iw[nw][no][i]); ++(wavxocn_jw[nw][no][i]); ++(wavxocn_io[nw][no][i]); wavxocn_jo[nw][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(nwxo[nw][no], wavxocn_area[nw][no], gdata_dbl); if(check) { for(n=0; n 1.e-5) { printf("at tile =%d, i=%d, j=%d, ratio=%g, area1=%g, area2=%g\n", n, i%nxa[n], i/nxa[n], cur_ratio, area_atm[n][i], atm_xarea[n][i] ); } if(cur_ratio > max_ratio) { max_ratio = cur_ratio; max_ta = n; max_ja = i/nxa[n]; max_ia = i%nxa[n]; } } } printf("The maximum area is at tile=%d, i=%d, j=%d, ratio=%g, area1=%g, area2=%g\n", max_ta, max_ia, max_ja, max_ratio, area_atm[max_ta][max_ja*nxa[max_ta]+max_ia], atm_xarea[max_ta][max_ja*nxa[max_ta]+max_ia]); /* for cubic grid, when number of model points is odd, some grid cell area will be negative, need to think about how to solve this issue in the future */ /* atm_area_sum = 4*M_PI*RADIUS*RADIUS; */ atm_area_sum = 0; atm_area_sum_nest = 0; for(n=0; n=0) { axo_area_frac_nest = axo_area_sum_nest/atm_area_sum_nest*100; axl_area_frac_nest = axl_area_sum_nest/atm_area_sum_nest*100; tiling_area_nest = (atm_area_sum_nest-axo_area_sum_nest-axl_area_sum_nest)/atm_area_sum_nest*100; printf("\nNOTE: axo_area_sum_nest is %f and ocean fraction is %f%%\n", axo_area_sum_nest, axo_area_frac_nest); printf("NOTE: axl_area_sum_nest is %f and land fraction is %f%%\n", axl_area_sum_nest, axl_area_frac_nest); printf("NOTE: tiling error for nest region is %f%%\n", tiling_area_nest ); } } } /*Fianlly create the coupler mosaic file mosaic_name.nc */ { int fid, dim_string, dim_axo, dim_lxo, dim_axl, dim_axw, dim_wxo, dims[4], n; size_t start[4], nwrite[4]; int id_lmosaic_dir, id_lmosaic_file, id_omosaic_dir, id_omosaic_file, id_wmosaic_dir; int id_amosaic_dir, id_amosaic_file, id_otopog_dir, id_otopog_file, id_wmosaic_file; int id_xgrids_dir, id_axo_file, id_lxo_file, id_axl_file, id_wxo_file; int id_amosaic, id_lmosaic, id_omosaic, id_wmosaic; fid = mpp_open(mosaic_file, MPP_WRITE); mpp_def_global_att(fid, "grid_version", grid_version); mpp_def_global_att(fid, "code_version", tagname); if( clip_method == GREAT_CIRCLE_CLIP) mpp_def_global_att(fid, "great_circle_algorithm", "TRUE"); mpp_def_global_att(fid, "history", history); dim_string = mpp_def_dim(fid, "string", STRING); dim_axo = mpp_def_dim(fid, "nfile_aXo", nfile_axo); dim_axl = mpp_def_dim(fid, "nfile_aXl", nfile_axl); dim_lxo = mpp_def_dim(fid, "nfile_lXo", nfile_lxo); if(nfile_wxo >0) dim_wxo = mpp_def_dim(fid, "nfile_wXo", nfile_wxo); id_amosaic_dir = mpp_def_var(fid, "atm_mosaic_dir", MPP_CHAR, 1, &dim_string, 1, "standard_name", "directory_storing_atmosphere_mosaic"); id_amosaic_file = mpp_def_var(fid, "atm_mosaic_file", MPP_CHAR, 1, &dim_string, 1, "standard_name", "atmosphere_mosaic_file_name"); id_amosaic = mpp_def_var(fid, "atm_mosaic", MPP_CHAR, 1, &dim_string, 1, "standard_name", "atmosphere_mosaic_name"); id_lmosaic_dir = mpp_def_var(fid, "lnd_mosaic_dir", MPP_CHAR, 1, &dim_string, 1, "standard_name", "directory_storing_land_mosaic"); id_lmosaic_file = mpp_def_var(fid, "lnd_mosaic_file", MPP_CHAR, 1, &dim_string, 1, "standard_name", "land_mosaic_file_name"); id_lmosaic = mpp_def_var(fid, "lnd_mosaic", MPP_CHAR, 1, &dim_string, 1, "standard_name", "land_mosaic_name"); id_omosaic_dir = mpp_def_var(fid, "ocn_mosaic_dir", MPP_CHAR, 1, &dim_string, 1, "standard_name", "directory_storing_ocean_mosaic"); id_omosaic_file = mpp_def_var(fid, "ocn_mosaic_file", MPP_CHAR, 1, &dim_string, 1, "standard_name", "ocean_mosaic_file_name"); id_omosaic = mpp_def_var(fid, "ocn_mosaic", MPP_CHAR, 1, &dim_string, 1, "standard_name", "ocean_mosaic_name"); if(wmosaic) { id_wmosaic_dir = mpp_def_var(fid, "wav_mosaic_dir", MPP_CHAR, 1, &dim_string, 1, "standard_name", "directory_storing_wave_mosaic"); id_wmosaic_file = mpp_def_var(fid, "wav_mosaic_file", MPP_CHAR, 1, &dim_string, 1, "standard_name", "wave_mosaic_file_name"); id_wmosaic = mpp_def_var(fid, "wav_mosaic", MPP_CHAR, 1, &dim_string, 1, "standard_name", "wave_mosaic_name"); } id_otopog_dir = mpp_def_var(fid, "ocn_topog_dir", MPP_CHAR, 1, &dim_string, 1, "standard_name", "directory_storing_ocean_topog"); id_otopog_file = mpp_def_var(fid, "ocn_topog_file", MPP_CHAR, 1, &dim_string, 1, "standard_name", "ocean_topog_file_name"); /* since exchange grid is created in this tool, we may add command line option to specify where to store the output */ /* id_xgrids_dir = mpp_def_var(fid, "xgrids_dir", MPP_CHAR, 1, &dim_string, 1, "standard_name", "directory_storing_xgrids"); */ dims[0] = dim_axo; dims[1] = dim_string; id_axo_file = mpp_def_var(fid, "aXo_file", MPP_CHAR, 2, dims, 1, "standard_name", "atmXocn_exchange_grid_file"); dims[0] = dim_axl; dims[1] = dim_string; id_axl_file = mpp_def_var(fid, "aXl_file", MPP_CHAR, 2, dims, 1, "standard_name", "atmXlnd_exchange_grid_file"); dims[0] = dim_lxo; dims[1] = dim_string; id_lxo_file = mpp_def_var(fid, "lXo_file", MPP_CHAR, 2, dims, 1, "standard_name", "lndXocn_exchange_grid_file"); if(nfile_wxo >0) { dims[0] = dim_wxo; dims[1] = dim_string; id_wxo_file = mpp_def_var(fid, "wXo_file", MPP_CHAR, 2, dims, 1, "standard_name", "wavXocn_exchange_grid_file"); } mpp_end_def(fid); for(i=0; i<4; i++) { start[i] = 0; nwrite[i] = 1; } nwrite[0] = strlen(amosaic_dir); mpp_put_var_value_block(fid, id_amosaic_dir, start, nwrite, amosaic_dir); nwrite[0] = strlen(amosaic_file); mpp_put_var_value_block(fid, id_amosaic_file, start, nwrite, amosaic_file); nwrite[0] = strlen(amosaic_name); mpp_put_var_value_block(fid, id_amosaic, start, nwrite, amosaic_name); nwrite[0] = strlen(lmosaic_dir); mpp_put_var_value_block(fid, id_lmosaic_dir, start, nwrite, lmosaic_dir); nwrite[0] = strlen(lmosaic_file); mpp_put_var_value_block(fid, id_lmosaic_file, start, nwrite, lmosaic_file); nwrite[0] = strlen(lmosaic_name); mpp_put_var_value_block(fid, id_lmosaic, start, nwrite, lmosaic_name); nwrite[0] = strlen(omosaic_dir); mpp_put_var_value_block(fid, id_omosaic_dir, start, nwrite, omosaic_dir); nwrite[0] = strlen(omosaic_file); mpp_put_var_value_block(fid, id_omosaic_file, start, nwrite, omosaic_file); nwrite[0] = strlen(omosaic_name); mpp_put_var_value_block(fid, id_omosaic, start, nwrite, omosaic_name); if(wmosaic) { nwrite[0] = strlen(wmosaic_dir); mpp_put_var_value_block(fid, id_wmosaic_dir, start, nwrite, wmosaic_dir); nwrite[0] = strlen(wmosaic_file); mpp_put_var_value_block(fid, id_wmosaic_file, start, nwrite, wmosaic_file); nwrite[0] = strlen(wmosaic_name); mpp_put_var_value_block(fid, id_wmosaic, start, nwrite, wmosaic_name); } nwrite[0] = strlen(otopog_dir); mpp_put_var_value_block(fid, id_otopog_dir, start, nwrite, otopog_dir); nwrite[0] = strlen(otopog_file); mpp_put_var_value_block(fid, id_otopog_file, start, nwrite, otopog_file); nwrite[0] = 1; for(n=0; n