/* * 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). Code here is non-reentrant for speed. * - Mike Winton (4/01) */ /* * to compile using libMPI: * * cc -O2 -o make_xgrids_parallel make_xgrids.c -Duse_libMPI -I/usr/local/include -L/usr/local/lib -lnetcdf -lmpi -lm * * to compile without using libMPI: * cc -O2 -o make_xgrids make_xgrids.c -I/usr/local/include -L/usr/local/lib -lnetcdf -lm */ #include #include #include #include #include #include #ifdef use_libMPI #include #endif /* This constant is used to control the area of exchange grids. */ #define EPS1 (1e-10) #define EPS3 (1e-4) #define MAXLOCAL 1e7 int pe, npes, root_pe; #ifdef use_libMPI MPI_Request request=MPI_REQUEST_NULL; #endif char *usage[] = { "", " __________________ MAKE_XGRIDS __________________", "", "Make_xgrids generates two exchange grids for the FMS coupler, one for fluxes", "between the atmosphere and surface (sea ice and land), the other for runoff", "between the land and sea ice. Make_xgrids expects a NetCDF format input", "specification for the ice/ocean grid. Three fields are required to be in", "this file:", "", " (1) wet - a 2D array of double precision numbers set to 1.0 where the ice", " and ocean models are active and 0.0 else where. WET has im indices in", " the i (pseudo east-west) direction and jm indices in the j (pseudo", " north-south) direction. These correspond to the global arrays of", " temperature, salinity and ice thickness in the coupled climate model.", "", " (2) geolon_vert_t and (3) geolat_vert_t - 2D double precision arrays", " (dimensioned im+1 by jm+1) that contain the longitudes and latitudes", " (respectively) of the corners of the sea ice and ocean model grid cells", " in degrees.", "", "Since the atmosphere and land models are required to be latitude/longitude,", "1D arrays of the latitude and longitude boundaries suffice to specify their", "grids. Make_xgrids expects this information to be provided in text files with", "a single latitude or longitude (in degrees) per line. In both the 2D and 1D", "boundary specification arrays, the coordinate values must increase", "monotonically with index.", "", "A sample call to make_xgrids that makes exchange grids for a tripolar grid", "sea ice/ocean, a T42 atmosphere, and a one-degree land model might look like:", "", "If compiled with -Duse_libMPI, the command will be:", "", " make_xgrids -o ocean_grid_spec.nc -a atmos_grid_spec.nc -l land_grid_spec.nc", "or", " make_xgrids -o tripolar_ocean.nc -a t42_x.dat,t42_y.dat -l 1deg_x.dat,1deg_y.dat", "", "", "If compiled without -Duse_libMPI, the command will be:", "", " mpirun -np # make_xgrids -o ocean_grid_spec.nc -a atmos_grid_spec.nc -l land_grid_spec.nc", "or", " mpirun -np # make_xgrids -o tripolar_ocean.nc -a t42_x.dat,t42_y.dat -l 1deg_x.dat,1deg_y.dat", "", "Make_xgrids copies all fields of the ocean grid specification file to its output", "file, \"grid_spec.nc\", and appends fields that specify the atmosphere and land", "model grids, and the surface and runoff exchange grids. The ocean model, sea ice", "model, land model and coupler expect to find the grid_spec.nc file in the INPUT", "subdirectory at run time. The coupler checks that the run time atmosphere model", "has the same resolution as that in grid_spec.nc.", "", " __________________ MAKE_XGRIDS __________________", "", NULL }; #define NC_CALL(F) \ { int status; \ if ((status=(F))!=NC_NOERR) { \ fprintf(stderr, "netCDF error: %s\n", nc_strerror(status)); \ exit(1); \ } \ } int n_max = 0; /* maximum number of vertices in x-cell */ #define MV 20 typedef struct { double x, y; } VTX; void error_handler(char *str) { fprintf( stderr, "%s\n", str); #ifdef use_libMPI MPI_Abort(MPI_COMM_WORLD, -1); #else exit(1); #endif } /* * Sutherland-Hodgeman algorithm sequentially clips parts outside 4 boundaries */ int clip(VTX v_in[], int n_in, VTX lolf, VTX uprt, VTX v_out[]) { VTX v_tmp[MV], v_last, v_isect; int i_in, i_out, n_out, inside_last, inside; /* clip polygon with LEFT boundary - clip V_IN to V_TMP */ v_last = v_in[n_in-1]; inside_last = (v_last.x >= lolf.x); for (i_in=0,i_out=0;i_in= lolf.x))!=inside_last) { v_isect.x = lolf.x; v_isect.y = v_last.y + (lolf.x - v_last.x) * (v_in[i_in].y - v_last.y) / (v_in[i_in].x - v_last.x); v_tmp[i_out++] = v_isect; } /* if "to" point is right of LEFT boundary, output it */ if (inside) v_tmp[i_out++] = v_in[i_in]; v_last = v_in[i_in]; inside_last = inside; } if (!(n_out=i_out)) return(0); /* clip polygon with RIGHT boundary - clip V_TMP to V_OUT */ v_last = v_tmp[n_out-1]; inside_last = (v_last.x <= uprt.x); for (i_in=0,i_out=0;i_in= lolf.y); for (i_in=0,i_out=0;i_in= lolf.y))!=inside_last) { v_isect.y = lolf.y; v_isect.x = v_last.x + (lolf.y - v_last.y) * (v_out[i_in].x - v_last.x) / (v_out[i_in].y - v_last.y); v_tmp[i_out++] = v_isect; } /* if "to" point is above BOTTOM boundary, output it */ if (inside) v_tmp[i_out++] = v_out[i_in]; v_last = v_out[i_in]; inside_last = inside; } if (!(n_out=i_out)) return(0); /* clip polygon with TOP boundary - clip V_TMP to V_OUT */ v_last = v_tmp[n_out-1]; inside_last = (v_last.y <= uprt.y); for (i_in=0,i_out=0;i_inn_max ? n_out : n_max; return(i_out); } /* clip */ #define TRAPEZOID_RULE 0 /* * poly_area - obtains area of input polygon by line integrating -sin(lat)d(lon) * * Vertex coordinates must be in degrees. * * Vertices must be listed counter-clockwise around polygon. * * Returns area in fraction of total sphere surface area. */ double poly_area(VTX v[], int n) { double area = 0.0; double D2R = M_PI/180; int i; 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 -= dx*sin(lat1); else #if TRAPEZOID_RULE area -= dx*(sin(lat1)+sin(lat2))/2; #else area += dx*(cos(lat1)-cos(lat2))/(lat1-lat2); #endif } return (area/(4*M_PI)); } /* poly_area */ /* This routine is used to calculate the integral which equal the product of latitude of the centroid and area of any grid box */ double poly_int3(VTX v[], int n) { double int3 = 0.0; double D2R = M_PI/180; int i; 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 */ int3 -= dx*(cos(lat1) + lat1*sin(lat1)); else #if TRAPEZOID_RULE int3 -= dx*(cos(lat1) + lat1*sin(lat1)+cos(lat2) + lat2*sin(lat2))/2.0; #else int3 -= dx*(2*sin(lat1)-lat1*cos(lat1) - 2*sin(lat2)+lat2*cos(lat2) )/(lat1-lat2); #endif } return (int3/(4*M_PI)); } /* poly_int3 */ /* this routine can be used to calculate the latitude of the centroid of any grid box */ double poly_ctrlat(VTX v[], int n) { double area; area = poly_area(v,n); if(area < EPS1) return (0.0); else return (poly_int3(v,n)/area); } /* poly_ctrlat */ /* This routine is used to calculate the integral which equal the product of lontitude of the centroid and area of any grid box */ double poly_int2(VTX v[], int n, double clon) { double int2 = 0.0; double D2R = M_PI/180; int i; clon = clon * D2R; for (i=0;i M_PI) dphi = dphi - 2.0*M_PI; if(dphi < -M_PI) dphi = dphi + 2.0*M_PI; dphi1 = phi1 - clon; if( dphi1 > M_PI) dphi1 -= 2.0*M_PI; if( dphi1 <-M_PI) dphi1 += 2.0*M_PI; dphi2 = phi2 -clon; if( dphi2 > M_PI) dphi2 -= 2.0*M_PI; if( dphi2 <-M_PI) dphi2 += 2.0*M_PI; if(abs(dphi2 -dphi1) < M_PI) { int2 -= dphi * (dphi1*f1+dphi2*f2)/2.0; } else { if(dphi1 > 0.0) fac = M_PI; else fac = -M_PI; fint = f1 + (f2-f1)*(fac-dphi1)/abs(dphi); int2 -= 0.5*dphi1*(dphi1-fac)*f1 - 0.5*dphi2*(dphi2+fac)*f2 + 0.5*fac*(dphi1+dphi2)*fint; } } return (int2/(4*M_PI)); } /* poly_int2 */ /* this routine can be used to calculate the lontitude of the centroid of any grid box */ double poly_ctrlon(VTX v[], int n, double clon) { double area; area = poly_area(v,n); if(area < EPS1) return (0.0); else return (poly_int2(v,n,clon)/area); } /* poly_ctrlon */ /* end change */ double box_area(VTX ll, VTX ur) { double D2R = M_PI/180; double dx = (ur.x-ll.x)*D2R; if(dx > M_PI) dx = dx - 2.0*M_PI; if(dx < -M_PI) dx = dx + 2.0*M_PI; return ( dx*(sin(ur.y*D2R)-sin(ll.y*D2R))/(4*M_PI) ); } /* box_area */ /* This routine is used to calculate the integral which equal the product of latitude of the centroid and area of a uniform grid box */ double box_int3(VTX ll, VTX ur) { double D2R = M_PI/180; double dphi = (ur.x-ll.x)*D2R; if(dphi > M_PI) dphi = dphi - 2.0*M_PI; if(dphi < -M_PI) dphi = dphi + 2.0*M_PI; return ( dphi*(cos(ur.y*D2R) + ur.y*D2R*sin(ur.y*D2R)-(cos(ll.y*D2R) + ll.y*D2R*sin(ll.y*D2R)))/(4*M_PI) ); } /* box_int3 */ /* this routine can be used to calculate the latitude of the centroid of a uniform grid box */ double box_ctrlat(VTX ll, VTX ur) { double area; area = box_area(ll,ur); if(area < EPS1) return (0.0); else return ( box_int3(ll,ur)/area); } /* box_ctrlat */ /* This routine is used to calculate the integral which equal the product of lontitude of the centroid and area of a uniform grid box */ double box_int2(VTX ll, VTX ur, double clon) { double phi1, phi2, dphi, lat1, lat2, dphi1, dphi2; double f1, f2, fac, fint; double D2R = M_PI/180; double int2 = 0.0; int i; clon = clon * D2R; for( i =0; i<2; i++) { if(i == 0) { phi1 = ur.x * D2R; phi2 = ll.x * D2R; lat1 = lat2 = ll.y * D2R; } else { phi1 = ll.x * D2R; phi2 = ur.x * D2R; lat1 = lat2 = ur.y * D2R; } dphi = phi1 - phi2; f1 = 0.5*(cos(lat1)*sin(lat1)+lat1); f2 = 0.5*(cos(lat2)*sin(lat2)+lat2); if(dphi > M_PI) dphi = dphi - 2.0*M_PI; if(dphi < -M_PI) dphi = dphi + 2.0*M_PI; /* make sure the center is in the same grid box. */ dphi1 = phi1 - clon; if( dphi1 > M_PI) dphi1 -= 2.0*M_PI; if( dphi1 <-M_PI) dphi1 += 2.0*M_PI; dphi2 = phi2 -clon; if( dphi2 > M_PI) dphi2 -= 2.0*M_PI; if( dphi2 <-M_PI) dphi2 += 2.0*M_PI; if(abs(dphi2 -dphi1) < M_PI) { int2 -= dphi * (dphi1*f1+dphi2*f2)/2.0; } else { if(dphi1 > 0.0) fac = M_PI; else fac = -M_PI; fint = f1 + (f2-f1)*(fac-dphi1)/abs(dphi); int2 -= 0.5*dphi1*(dphi1-fac)*f1 - 0.5*dphi2*(dphi2+fac)*f2 + 0.5*fac*(dphi1+dphi2)*fint; } } return (int2/(4*M_PI)); } /* box_int2 */ /* this routine can be used to calculate the lontitude of the centroid of a uniform grid box */ double box_ctrlon(VTX ll, VTX ur, double clon) { double area; area = box_area(ll,ur); if(area < EPS1) return (0.0); else return(box_int2(ll,ur,clon)/area); } /* box_ctrlon */ int vtx_delete(VTX v[], int n, int n_del) { for (;n_del=n_ins;i--) v[i+1] = v[i]; v[n_ins].x = x; v[n_ins].y = y; return (n+1); } /* vtx_insert */ void v_print(VTX v[], int n) { int i; for (i=0;i=90.0-EPS3) pole = 1; if (0&&pole) { printf("fixing pole cell\n"); v_print(v,nn); printf("---------"); } /* all pole points must be paired */ for (i=0;i=90.0-EPS3) { int im=(i+nn-1)%nn, ip=(i+1)%nn; if (v[im].y==v[i].y && v[ip].y==v[i].y) { nn = vtx_delete(v, nn, i); i--; } else if (v[im].y!=v[i].y && v[ip].y!=v[i].y) { nn = vtx_insert(v, nn, i, v[i].x, v[i].y); i++; } } /* first of pole pair has longitude of previous vertex */ /* second of pole pair has longitude of subsequent vertex */ for (i=0;i=90.0-EPS3) { int im=(i+nn-1)%nn, ip=(i+1)%nn; if (v[im].y!=v[i].y) v[i].x = v[im].x; if (v[ip].y!=v[i].y) v[i].x = v[ip].x; } if (nn) x_sum = v[0].x; else return(0); for (i=1;i 180) dx = dx - 360; x_sum += (v[i].x = v[i-1].x + dx); } dx = (x_sum/nn)-tlon; if (dx < -180) for (i=0;i 180) for (i=0;i= nfiles ) return 0; if(n == 0) { strcpy(curfile, filename); } else { pathlength = strlen(filename) - 3; if(strncmp( filename+pathlength, ".nc", 3) != 0) error_handler("The suffix of the input file is not .nc"); strncpy(curfile, filename, pathlength); curfile[pathlength] = '\0'; sprintf(str, "%d", n); strcat(curfile, str); strcat(curfile, ".nc"); } NC_CALL(nc_open(curfile, NC_NOWRITE, ncid)); status = nc_inq_varid(*ncid, fieldname, varid); if (status == NC_NOERR) return 1; NC_CALL(nc_close(*ncid)); n++; } return 0; } /* get_ncid_varid */ void *get_var2d_field3d(fname, name, nfiles, size1ptr, size2ptr) char *fname; char *name; int nfiles; int *size1ptr, *size2ptr; { double *var; double *result; int ncid, ncvid, dims[4]; int status, i,j,m,n; size_t im, jm, km; status = get_ncid_varid( fname, name, nfiles, &ncid, &ncvid ); if(status == 0) return 0; NC_CALL(nc_inq_vardimid(ncid, ncvid, dims)) NC_CALL(nc_inq_dimlen(ncid, dims[0], &km)) NC_CALL(nc_inq_dimlen(ncid, dims[1], &jm)) NC_CALL(nc_inq_dimlen(ncid, dims[2], &im)) if( km != 4) { fprintf(stderr,"each T-cell should have 4 vertices %s\n", "get_nc_var2d"); exit(1); } *size1ptr = im; *size2ptr = jm; var = (double *) malloc ( im*jm*4*sizeof(double) ); result = (double *) malloc ( (im+1)*(jm+1)*sizeof(double) ); if (nc_get_var_double(ncid, ncvid, var)==-1) { fprintf(stderr,"Couldn't get variable %s\n", name); exit(1); } NC_CALL(nc_close(ncid)) m=0; for(j=0;j(Y)?(X):(Y)) int nums_from_file(fname, nums) char *fname; double *nums[]; { FILE *fp = fopen(fname, "r"); double d[MN]; int i, n; for (n=0;n npts ) error_handler("define_domain: number of divisions is larger than number of points"); npts_left = npts; start = 0; for(n=0; n MAXLOCAL) error_handler(" number of atmos points on this pe is greater MAXLOCAL, increase MAXLOCAL"); atm_area = (double *) malloc (MAXLOCAL*sizeof(double)); atm_ctrlon = (double *) malloc (MAXLOCAL*sizeof(double)); atm_ctrlat = (double *) malloc (MAXLOCAL*sizeof(double)); atm_area_left = (double *) malloc (MAXLOCAL*sizeof(double)); for (na = 0; na< natm_local; na ++) { VTX ll, ur; ia = (na+isc)%IMA; ja = (na+isc)/IMA; ll.x = atmlonb[ia ]; ll.y = atmlatb[ja ]; ur.x = atmlonb[ia+1]; ur.y = atmlatb[ja+1]; atm_area_left[na] = atm_area[na] = box_area(ll, ur); atm_ctrlon[na] = box_ctrlon(ll, ur, 0.5*(ll.x+ur.x)); atm_ctrlat[na] = box_ctrlat(ll, ur); } /* get global atm_area */ gatm_area = (double *)malloc(IMA*JMA*sizeof(double)); get_global_data_double(atm_area, gatm_area, natm_local); lnd_area = (double *) malloc(IML*JML*sizeof(double)); for (i=0;i0.5) { VTX v[MV]; int n_in; /* list ocean cell vertices in counter-clockwise order */ v[0].y = ocnlatb[ jo *(IMO+1)+io ]; /* SW - y */ v[1].y = ocnlatb[ jo *(IMO+1)+io+1]; /* SE - y */ v[2].y = ocnlatb[(jo+1)*(IMO+1)+io+1]; /* NE - y */ v[3].y = ocnlatb[(jo+1)*(IMO+1)+io ]; /* NW - y */ v[0].x = ocnlonb[ jo *(IMO+1)+io ]; /* SW - x */ v[1].x = ocnlonb[ jo *(IMO+1)+io+1]; /* SE - x */ v[2].x = ocnlonb[(jo+1)*(IMO+1)+io+1]; /* NE - x */ v[3].x = ocnlonb[(jo+1)*(IMO+1)+io ]; /* NW - x */ n_in = lon_fix(v, 4, 180.0); if ((ocn_minus_runoff[jo*IMO+io] = ocn_area_left[jo*IMO+io] = ocn_area[jo*IMO+io] = poly_area(v, n_in))<0) { fprintf(stderr, "\nOcean cell at i=%d j=%d has negative area!\n", io, jo); v_print(v, 4); fprintf(stderr,"\n"); exit(1); } } else ocn_minus_runoff[jo*IMO+io] = ocn_area_left[jo*IMO+io] = ocn_area[jo*IMO+io] = 0.0; /* * Define output file variables for atmos/ocean and atmos/land X-grids */ NC_CALL(nc_create(axo_file, NC_NOCLOBBER, &ncidAO)) NC_CALL(nc_def_dim(ncidAO, "i_atmXocn", NC_UNLIMITED, &dim_AO)) NC_CALL(nc_def_var(ncidAO, "AREA_ATMxOCN", NC_DOUBLE,1, &dim_AO, &area_AO_id)) NC_CALL(nc_def_var(ncidAO, "DI_ATMxOCN", NC_DOUBLE,1, &dim_AO, &di_AO_id)) NC_CALL(nc_def_var(ncidAO, "DJ_ATMxOCN", NC_DOUBLE,1, &dim_AO, &dj_AO_id)) NC_CALL(nc_def_var(ncidAO, "I_ATM_ATMxOCN", NC_INT, 1, &dim_AO, &ia_AO_id)) NC_CALL(nc_def_var(ncidAO, "J_ATM_ATMxOCN", NC_INT, 1, &dim_AO, &ja_AO_id)) NC_CALL(nc_def_var(ncidAO, "I_OCN_ATMxOCN", NC_INT, 1, &dim_AO, &io_AO_id)) NC_CALL(nc_def_var(ncidAO, "J_OCN_ATMxOCN", NC_INT, 1, &dim_AO, &jo_AO_id)) NC_CALL(nc_enddef(ncidAO)) NC_CALL(nc_create(axl_file, NC_NOCLOBBER, &ncidAL)) NC_CALL(nc_def_dim(ncidAL, "I_ATMxLND", NC_UNLIMITED, &dim_AL)) NC_CALL(nc_def_var(ncidAL, "AREA_ATMxLND", NC_DOUBLE,1, &dim_AL, &area_AL_id)) NC_CALL(nc_def_var(ncidAL, "DI_ATMxLND", NC_DOUBLE,1, &dim_AL, &di_AL_id)) NC_CALL(nc_def_var(ncidAL, "DJ_ATMxLND", NC_DOUBLE,1, &dim_AL, &dj_AL_id)) NC_CALL(nc_def_var(ncidAL, "I_ATM_ATMxLND", NC_INT, 1, &dim_AL, &ia_AL_id)) NC_CALL(nc_def_var(ncidAL, "J_ATM_ATMxLND", NC_INT, 1, &dim_AL, &ja_AL_id)) NC_CALL(nc_def_var(ncidAL, "I_LND_ATMxLND", NC_INT, 1, &dim_AL, &il_AL_id)) NC_CALL(nc_def_var(ncidAL, "J_LND_ATMxLND", NC_INT, 1, &dim_AL, &jl_AL_id)) NC_CALL(nc_enddef(ncidAL)) /* * Define output file variables for land/ocean X-grid */ NC_CALL(nc_create(lxo_file, NC_NOCLOBBER, &ncidLO)) NC_CALL(nc_def_dim(ncidLO, "i_lndXocn", NC_UNLIMITED, &dim_LO)) NC_CALL(nc_def_var(ncidLO, "AREA_LNDxOCN", NC_DOUBLE,1, &dim_LO, &area_LO_id)) NC_CALL(nc_def_var(ncidLO, "DI_LNDxOCN", NC_DOUBLE,1, &dim_LO, &di_LO_id)) NC_CALL(nc_def_var(ncidLO, "DJ_LNDxOCN", NC_DOUBLE,1, &dim_LO, &dj_LO_id)) NC_CALL(nc_def_var(ncidLO, "I_LND_LNDxOCN", NC_INT, 1, &dim_LO, &il_LO_id)) NC_CALL(nc_def_var(ncidLO, "J_LND_LNDxOCN", NC_INT, 1, &dim_LO, &jl_LO_id)) NC_CALL(nc_def_var(ncidLO, "I_OCN_LNDxOCN", NC_INT, 1, &dim_LO, &io_LO_id)) NC_CALL(nc_def_var(ncidLO, "J_OCN_LNDxOCN", NC_INT, 1, &dim_LO, &jo_LO_id)) NC_CALL(nc_enddef(ncidLO)) /* * Find overlaps of atmosphere cells (the window) and surface cells * Find the information about land. */ { int atmXlnd_n = 0; int atmXocn_n = 0; int atmXlnd_start, atmXlnd_n_prev; double *atmXlnd_area, *atmXlnd_int2, *atmXlnd_int3; double *atmXlnd_di, *atmXlnd_dj; double *atmXocn_area, *atmXocn_di, *atmXocn_dj; double *gal_area=0, *gal_di=0, *gal_dj=0; int *gal_il=0, *gal_jl=0, *gal_ia=0, *gal_ja=0; double *gao_area=0, *gao_di=0, *gao_dj=0; int *gao_io=0, *gao_jo=0, *gao_ia=0, *gao_ja=0; int *atmXocn_io, *atmXocn_jo; int *atmXocn_ia, *atmXocn_ja; int *atmXlnd_il, *atmXlnd_jl; int *atmXlnd_ia, *atmXlnd_ja; VTX *atmXlnd_ll, *atmXlnd_ur; atmXlnd_area = (double *) malloc (MAXLOCAL*sizeof(double)); atmXlnd_di = (double *) malloc (MAXLOCAL*sizeof(double)); atmXlnd_dj = (double *) malloc (MAXLOCAL*sizeof(double)); atmXlnd_int2 = (double *) malloc (MAXLOCAL*sizeof(double)); atmXlnd_int3 = (double *) malloc (MAXLOCAL*sizeof(double)); atmXlnd_ll = (VTX *) malloc (MAXLOCAL*sizeof( VTX)); atmXlnd_ur = (VTX *) malloc (MAXLOCAL*sizeof( VTX)); atmXlnd_il = (int *) malloc (MAXLOCAL*sizeof( int)); atmXlnd_jl = (int *) malloc (MAXLOCAL*sizeof( int)); atmXlnd_ia = (int *) malloc (MAXLOCAL*sizeof( int)); atmXlnd_ja = (int *) malloc (MAXLOCAL*sizeof( int)); atmXocn_area = (double *) malloc (MAXLOCAL*sizeof(double)); atmXocn_di = (double *) malloc (MAXLOCAL*sizeof(double)); atmXocn_dj = (double *) malloc (MAXLOCAL*sizeof(double)); atmXocn_io = (int *) malloc (MAXLOCAL*sizeof( int)); atmXocn_jo = (int *) malloc (MAXLOCAL*sizeof( int)); atmXocn_ia = (int *) malloc (MAXLOCAL*sizeof( int)); atmXocn_ja = (int *) malloc (MAXLOCAL*sizeof( int)); for (na = 0; na < natm_local; na ++) { VTX atm_ll, atm_ur; atmXlnd_start = atmXlnd_n; ia = (na+isc)%IMA; ja = (na+isc)/IMA; /* atmosphere cell clipping window */ atm_ll.x = atmlonb[ia ]; atm_ll.y = atmlatb[ja ]; atm_ur.x = atmlonb[ia+1]; atm_ur.y = atmlatb[ja+1]; for (jl=0;jl atm_ur.x) { lnd_ll.x -= 360; lnd_ur.x -= 360; } if(atm_ll.x > lnd_ur.x) { lnd_ll.x += 360; lnd_ur.x += 360; } atmXlnd_ll[atmXlnd_n].x = MAX(atm_ll.x,lnd_ll.x); atmXlnd_ur[atmXlnd_n].x = MIN(atm_ur.x,lnd_ur.x); if (atmXlnd_ll[atmXlnd_n].x>atmXlnd_ur[atmXlnd_n].x) continue; atmXlnd_ll[atmXlnd_n].y = MAX(atm_ll.y,lnd_ll.y); atmXlnd_ur[atmXlnd_n].y = MIN(atm_ur.y,lnd_ur.y); if (atmXlnd_ll[atmXlnd_n].y>atmXlnd_ur[atmXlnd_n].y) continue; atmXlnd_area[atmXlnd_n] = box_area(atmXlnd_ll[atmXlnd_n], atmXlnd_ur[atmXlnd_n]); atmXlnd_int2[atmXlnd_n] = box_int2(atmXlnd_ll[atmXlnd_n], atmXlnd_ur[atmXlnd_n], 0.5*(atm_ll.x+atm_ur.x)); atmXlnd_int3[atmXlnd_n] = box_int3(atmXlnd_ll[atmXlnd_n], atmXlnd_ur[atmXlnd_n]); atmXlnd_il[atmXlnd_n] = il+1; atmXlnd_jl[atmXlnd_n] = jl+1; atmXlnd_ia[atmXlnd_n] = ia+1; atmXlnd_ja[atmXlnd_n] = ja+1; atmXlnd_n++; if(atmXlnd_n > MAXLOCAL) error_handler("atmXlnd_n is greater than MAXLOCAL, increase MAXLOCAL"); }; /* calculate atmos/ocean x-cells */ for (jo=0;jo0.5) { int n_in, n_out; VTX v_in[MV], v_out[MV]; double Xarea, Xctrlon, Xctrlat; /* list ocean cell vertices in counter-clockwise order */ v_in[0].y = ocnlatb[ jo *(IMO+1)+io ]; /* SW - y */ v_in[1].y = ocnlatb[ jo *(IMO+1)+io+1]; /* SE - y */ v_in[2].y = ocnlatb[(jo+1)*(IMO+1)+io+1]; /* NE - y */ v_in[3].y = ocnlatb[(jo+1)*(IMO+1)+io ]; /* NW - y */ if ( (v_in[0].y<=atm_ll.y) && (v_in[1].y<=atm_ll.y) && (v_in[2].y<=atm_ll.y) && (v_in[3].y<=atm_ll.y) ) continue; if ( (v_in[0].y>=atm_ur.y) && (v_in[1].y>=atm_ur.y) && (v_in[2].y>=atm_ur.y) && (v_in[3].y>=atm_ur.y) ) continue; v_in[0].x = ocnlonb[ jo *(IMO+1)+io ]; /* SW - x */ v_in[1].x = ocnlonb[ jo *(IMO+1)+io+1]; /* SE - x */ v_in[2].x = ocnlonb[(jo+1)*(IMO+1)+io+1]; /* NE - x */ v_in[3].x = ocnlonb[(jo+1)*(IMO+1)+io ]; /* NW - x */ n_in = lon_fix(v_in, 4, (atm_ll.x+atm_ur.x)/2); if ( ((n_out = clip ( v_in, n_in, atm_ll, atm_ur, v_out )) > 0) && ((Xarea = poly_area ( v_out, n_out )) > EPS1) ) { if (ocn_area[jo*IMO+io] < EPS1) { fprintf(stderr, "\nOcean cell at i=%d j=%d has invalid area!\n", io, jo); fprintf(stderr,"\n"); exit(1); } if (atm_area[na] < EPS1) { fprintf(stderr, "\nAtmosphter cell at i=%d j=%d has invalid area!\n", ia, ja); fprintf(stderr,"\n"); exit(1); } Xctrlon = poly_ctrlon ( v_out, n_out, 0.5*(atm_ll.x + atm_ur.x)); Xctrlat = poly_ctrlat ( v_out, n_out ); atmXocn_di[atmXocn_n] = Xctrlon - atm_ctrlon[na]; atmXocn_dj[atmXocn_n] = Xctrlat - atm_ctrlat[na]; atmXocn_area[atmXocn_n] = Xarea; atmXocn_io[atmXocn_n] = io+1; atmXocn_jo[atmXocn_n] = jo+1; atmXocn_ia[atmXocn_n] = ia+1; atmXocn_ja[atmXocn_n] = ja+1; atmXocn_n++; if(atmXocn_n > MAXLOCAL) error_handler("atmXocn_n is greater than MAXLOCAL, increase MAXLOCAL"); ocn_area_left[jo*IMO+io] -= Xarea; atm_area_left[na] -= Xarea; tot_area_AXO += Xarea; } /* remove ocean parts from atmos/land x-cells */ for (i=atmXlnd_start;i 0) && ( (Xarea = poly_area ( v_out, n_out )) > EPS1) ) { atmXlnd_area[i] -= Xarea; atmXlnd_int2[i] -= poly_int2 ( v_out, n_out, 0.5*(atm_ll.x+ atm_ur.x)); atmXlnd_int3[i] -= poly_int3 ( v_out, n_out ); } } /* ocean loop */ atmXlnd_n_prev = atmXlnd_n; atmXlnd_n = atmXlnd_start; for (i=atmXlnd_start;iEPS1) { /*calculate the centroid of amt/lnd xgrid. */ atmXlnd_di[atmXlnd_n] = atmXlnd_int2[i] /atmXlnd_area[i] - atm_ctrlon[na]; atmXlnd_dj[atmXlnd_n] = atmXlnd_int3[i] /atmXlnd_area[i] - atm_ctrlat[na]; atm_area_left[na] -= atmXlnd_area[i]; lnd_area[(atmXlnd_jl[i]-1)*IML+atmXlnd_il[i]-1] += atmXlnd_area[i]; tot_area_AXL += atmXlnd_area[i]; atmXlnd_area[atmXlnd_n] = atmXlnd_area[i]; atmXlnd_il[atmXlnd_n] = atmXlnd_il[i]; atmXlnd_jl[atmXlnd_n] = atmXlnd_jl[i]; atmXlnd_ia[atmXlnd_n] = atmXlnd_ia[i]; atmXlnd_ja[atmXlnd_n] = atmXlnd_ja[i]; atmXlnd_n++; } if( npes == 1 ) { /* print this information in parallel run is complicated */ if (ia+1==IMA) printf( "to %6.2f: %d AL x-cells (area=%g); %d AO x-cells (area=%g)\n", atmlatb[ja+1], atmXlnd_n, tot_area_AXL, atmXocn_n, tot_area_AXO); } } /* atmosphere loop */ iAXL = get_global_count( atmXlnd_n ); iAXO = get_global_count( atmXocn_n ); sum_over_pe(&tot_area_AXL, 1); sum_over_pe(&tot_area_AXO, 1); if(pe == root_pe) printf( "\n%d AL x-cells (area=%g); %d AO x-cells (area=%g)\n\n", (int)iAXL, tot_area_AXL, (int)iAXO, tot_area_AXO ); /*write out atm/land xgrid, first need to send the data to root pe. write out only from root pe*/ if ( pe == root_pe ) { gal_area = (double *) malloc( iAXL * sizeof(double) ); gal_di = (double* ) malloc( iAXL * sizeof(double) ); gal_dj = (double* ) malloc( iAXL * sizeof(double) ); gal_ia = (int* ) malloc( iAXL * sizeof(int) ); gal_ja = (int* ) malloc( iAXL * sizeof(int) ); gal_il = (int* ) malloc( iAXL * sizeof(int) ); gal_jl = (int* ) malloc( iAXL * sizeof(int) ); } get_global_data_double( atmXlnd_area, gal_area, atmXlnd_n ); get_global_data_double( atmXlnd_di, gal_di, atmXlnd_n ); get_global_data_double( atmXlnd_dj, gal_dj, atmXlnd_n ); get_global_data_int( atmXlnd_il, gal_il, atmXlnd_n ); get_global_data_int( atmXlnd_jl, gal_jl, atmXlnd_n ); get_global_data_int( atmXlnd_ia, gal_ia, atmXlnd_n ); get_global_data_int( atmXlnd_ja, gal_ja, atmXlnd_n ); /* write out from root pe */ if( pe == root_pe ) { size_t start, nwrite; start = 0; nwrite = iAXL; nc_put_vara_double( ncidAL, area_AL_id, &start, &nwrite, gal_area ); nc_put_vara_double( ncidAL, di_AL_id, &start, &nwrite, gal_di ); nc_put_vara_double( ncidAL, dj_AL_id, &start, &nwrite, gal_dj ); nc_put_vara_int( ncidAL, ia_AL_id, &start, &nwrite, gal_ia ); nc_put_vara_int( ncidAL, ja_AL_id, &start, &nwrite, gal_ja ); nc_put_vara_int( ncidAL, il_AL_id, &start, &nwrite, gal_il ); nc_put_vara_int( ncidAL, jl_AL_id, &start, &nwrite, gal_jl ); /* release memory */ free(gal_area); free(gal_di); free(gal_dj); free(gal_ia); free(gal_ja); free(gal_il); free(gal_jl); } /*write out atm/ocn xgrid, first need to send the data to root pe. write out only from root pe*/ if ( pe == root_pe ) { gao_area = (double *) malloc( iAXO * sizeof(double) ); gao_di = (double* ) malloc( iAXO * sizeof(double) ); gao_dj = (double* ) malloc( iAXO * sizeof(double) ); gao_io = (int* ) malloc( iAXO * sizeof(int) ); gao_jo = (int* ) malloc( iAXO * sizeof(int) ); gao_ia = (int* ) malloc( iAXO * sizeof(int) ); gao_ja = (int* ) malloc( iAXO * sizeof(int) ); } get_global_data_double( atmXocn_area, gao_area, atmXocn_n ); get_global_data_double( atmXocn_di, gao_di, atmXocn_n ); get_global_data_double( atmXocn_dj, gao_dj, atmXocn_n ); get_global_data_int( atmXocn_ia, gao_ia, atmXocn_n ); get_global_data_int( atmXocn_ja, gao_ja, atmXocn_n ); get_global_data_int( atmXocn_io, gao_io, atmXocn_n ); get_global_data_int( atmXocn_jo, gao_jo, atmXocn_n ); /* write out from root pe */ if( pe == root_pe ) { size_t start, nwrite; start = 0; nwrite = iAXO; nc_put_vara_double( ncidAO, area_AO_id, &start, &nwrite, gao_area ); nc_put_vara_double( ncidAO, di_AO_id, &start, &nwrite, gao_di ); nc_put_vara_double( ncidAO, dj_AO_id, &start, &nwrite, gao_dj ); nc_put_vara_int( ncidAO, ia_AO_id, &start, &nwrite, gao_ia ); nc_put_vara_int( ncidAO, ja_AO_id, &start, &nwrite, gao_ja ); nc_put_vara_int( ncidAO, io_AO_id, &start, &nwrite, gao_io ); nc_put_vara_int( ncidAO, jo_AO_id, &start, &nwrite, gao_jo ); if(land_same_as_atmos == 1) { iLXO = iAXO; nc_put_vara_double( ncidLO, area_LO_id, &start, &nwrite, gao_area ); nc_put_vara_double( ncidLO, di_LO_id, &start, &nwrite, gao_di ); nc_put_vara_double( ncidLO, dj_LO_id, &start, &nwrite, gao_dj ); nc_put_vara_int( ncidLO, il_LO_id, &start, &nwrite, gao_ia ); nc_put_vara_int( ncidLO, jl_LO_id, &start, &nwrite, gao_ja ); nc_put_vara_int( ncidLO, io_LO_id, &start, &nwrite, gao_io ); nc_put_vara_int( ncidLO, jo_LO_id, &start, &nwrite, gao_jo ); } /* release memory */ free(gao_area); free(gao_di); free(gao_dj); free(gao_ia); free(gao_ja); free(gao_io); free(gao_jo); } /* release memory */ free(atmXocn_area); free(atmXocn_di); free(atmXocn_dj); free(atmXocn_ia); free(atmXocn_ja); free(atmXocn_io); free(atmXocn_jo); free(atmXlnd_area); free(atmXlnd_di); free(atmXlnd_dj); free(atmXlnd_ia); free(atmXlnd_ja); free(atmXlnd_il); free(atmXlnd_jl); /* get global lnd_area by summing over all the pe. */ sum_over_pe(lnd_area, IML*JML); } /* * Find overlaps of land cells (the window) and ocean cells for runoff x-grid * when land and atmosphere are on different grid. */ if( land_same_as_atmos == 0 ) { int lndXocn_n=0; int nlnd_local, nl; double *lndXocn_area, *lndXocn_di, *lndXocn_dj; double *glo_area, *glo_di, *glo_dj; int *lndXocn_il, *lndXocn_jl, *lndXocn_io, *lndXocn_jo; int *glo_il, *glo_jl, *glo_io, *glo_jo; define_domain( IML*JML, npes, startlist, endlist); isc = startlist[pe]; iec = endlist[pe]; nlnd_local = iec - isc + 1; if(nlnd_local > MAXLOCAL) error_handler("number of local land points is greater than MAXLOCAL, increase MAXLOCAL"); lndXocn_area = (double *) malloc (MAXLOCAL*sizeof(double)); lndXocn_di = (double *) malloc (MAXLOCAL*sizeof(double)); lndXocn_dj = (double *) malloc (MAXLOCAL*sizeof(double)); lndXocn_io = (int *) malloc (MAXLOCAL*sizeof( int)); lndXocn_jo = (int *) malloc (MAXLOCAL*sizeof( int)); lndXocn_il = (int *) malloc (MAXLOCAL*sizeof( int)); lndXocn_jl = (int *) malloc (MAXLOCAL*sizeof( int)); for (nl=0;nl0.5) { int n_in, n_out; VTX v_in[MV], v_out[MV]; double Xarea, Xctrlon, Xctrlat; /* list ocean cell vertices in counter-clockwise order */ v_in[0].y = ocnlatb[ jo *(IMO+1)+io ]; /* SW - y */ v_in[1].y = ocnlatb[ jo *(IMO+1)+io+1]; /* SE - y */ v_in[2].y = ocnlatb[(jo+1)*(IMO+1)+io+1]; /* NE - y */ v_in[3].y = ocnlatb[(jo+1)*(IMO+1)+io ]; /* NW - y */ if ( (v_in[0].y<=lnd_ll.y) && (v_in[1].y<=lnd_ll.y) && (v_in[2].y<=lnd_ll.y) && (v_in[3].y<=lnd_ll.y) ) continue; if ( (v_in[0].y>=lnd_ur.y) && (v_in[1].y>=lnd_ur.y) && (v_in[2].y>=lnd_ur.y) && (v_in[3].y>=lnd_ur.y) ) continue; v_in[0].x = ocnlonb[ jo *(IMO+1)+io ]; /* SW - x */ v_in[1].x = ocnlonb[ jo *(IMO+1)+io+1]; /* SE - x */ v_in[2].x = ocnlonb[(jo+1)*(IMO+1)+io+1]; /* NE - x */ v_in[3].x = ocnlonb[(jo+1)*(IMO+1)+io ]; /* NW - x */ n_in = lon_fix(v_in, 4, (lnd_ll.x+lnd_ur.x)/2); if ( ((n_out = clip ( v_in, n_in, lnd_ll, lnd_ur, v_out )) > 0) && ((Xarea = poly_area ( v_out, n_out )) > EPS1) ) { /* calculate the centroid of lnd/ocn runoff */ Xctrlon = poly_ctrlon ( v_out, n_out, 0.5*(lnd_ll.x+lnd_ur.x )); Xctrlat = poly_ctrlat ( v_out, n_out ); lndXocn_area[lndXocn_n] = Xarea; lndXocn_di[lndXocn_n] = Xctrlon - lnd_cell_ctrlon[jl*IML+il]; lndXocn_dj[lndXocn_n] = Xctrlat - lnd_cell_ctrlat[jl*IML+il]; lndXocn_il[lndXocn_n] = il + 1; lndXocn_jl[lndXocn_n] = jl + 1; lndXocn_io[lndXocn_n] = io + 1; lndXocn_jo[lndXocn_n] = jo + 1; lndXocn_n++; ocn_minus_runoff[jo*IMO+io] -= Xarea; } } /* ocean loop */ } /* land loop */ iLXO = get_global_count( lndXocn_n ); if ( pe == root_pe ) { glo_area = (double *) malloc( iLXO * sizeof(double) ); glo_di = (double* ) malloc( iLXO * sizeof(double) ); glo_dj = (double* ) malloc( iLXO * sizeof(double) ); glo_il = (int* ) malloc( iLXO * sizeof(int) ); glo_jl = (int* ) malloc( iLXO * sizeof(int) ); glo_io = (int* ) malloc( iLXO * sizeof(int) ); glo_jo = (int* ) malloc( iLXO * sizeof(int) ); } get_global_data_double( lndXocn_area, glo_area, lndXocn_n ); get_global_data_double( lndXocn_di, glo_di, lndXocn_n ); get_global_data_double( lndXocn_dj, glo_dj, lndXocn_n ); get_global_data_int( lndXocn_il, glo_il, lndXocn_n ); get_global_data_int( lndXocn_jl, glo_jl, lndXocn_n ); get_global_data_int( lndXocn_io, glo_io, lndXocn_n ); get_global_data_int( lndXocn_jo, glo_jo, lndXocn_n ); /* write out from root pe */ if( pe == root_pe ) { size_t start, nwrite; start = 0; nwrite = iLXO; nc_put_vara_double( ncidLO, area_LO_id, &start, &nwrite, glo_area ); nc_put_vara_double( ncidLO, di_LO_id, &start, &nwrite, glo_di ); nc_put_vara_double( ncidLO, dj_LO_id, &start, &nwrite, glo_dj ); nc_put_vara_int( ncidLO, il_LO_id, &start, &nwrite, glo_il ); nc_put_vara_int( ncidLO, jl_LO_id, &start, &nwrite, glo_jl ); nc_put_vara_int( ncidLO, io_LO_id, &start, &nwrite, glo_io ); nc_put_vara_int( ncidLO, jo_LO_id, &start, &nwrite, glo_jo ); /* release memory */ free(glo_area); free(glo_di); free(glo_dj); free(glo_io); free(glo_jo); free(glo_il); free(glo_jl); } free(lndXocn_area); free(lndXocn_di); free(lndXocn_dj); free(lndXocn_il); free(lndXocn_jl); free(lndXocn_io); free(lndXocn_jo); printf( "Land/Ocean x-grid complete (%d x-cells)\n", (int) iLXO); } { double a_err = 0.0, o_err = 0.0; double a_max_err = 0.0, o_max_err = 0.0; int o_cnt = 0; double *ga_area=NULL, *ga_area_left=NULL; ga_area = (double *)malloc(IMA*JMA*sizeof(double) ); ga_area_left = (double *)malloc(IMA*JMA*sizeof(double) ); get_global_data_double( atm_area_left, ga_area_left, natm_local ); get_global_data_double( atm_area, ga_area, natm_local ); if( pe == root_pe ) { for (i=0;ia_max_err*a_max_err) { a_max_err = -err; ja_max_err = i/IMA; ia_max_err = i-ja_max_err*IMA; } a_err += err2; } printf( "\nRMS Atmosphere cell tiling error = %lg\n", sqrt(a_err/(IMA*JMA))); printf( "MAX Atmosphere cell tiling error = %lg ", a_max_err); printf( "at (i=%d,j=%d), (x=%6.2f,y=%6.2f)\n", ia_max_err+1, ja_max_err+1, (atmlonb[ia_max_err]+atmlonb[ia_max_err+1])/2, (atmlatb[ja_max_err]+atmlatb[ja_max_err+1])/2); } free(ga_area); free(ga_area_left); if (ofile) { if( npes > 1 ) { /* need to sum over all the pe and minus (n-1)*ocn_area */ sum_over_pe( ocn_area_left, IMO*JMO ); for(i=0; i0.5) { double err = ocn_area_left[i]/ocn_area[i]; double err2 = err*err; if (err2>o_max_err*o_max_err) { o_max_err = -err; jo_max_err = i/IMO; io_max_err = i-jo_max_err*IMO; } o_err += err2; o_cnt++; } if(o_cnt>0) printf( "\nRMS Ocean cell tiling error (atmos/ocean x-cells) = %lg\n", sqrt(o_err/o_cnt)); printf( "MAX Ocean cell tiling error (atmos/ocean x-cells) = %lg ", o_max_err); printf( "at (i=%d,j=%d), (x=%6.2f,y=%6.2f)\n", io_max_err+1, jo_max_err+1, (ocnlonb[ jo_max_err *(IMO+1)+io_max_err ] +ocnlonb[ jo_max_err *(IMO+1)+io_max_err+1] +ocnlonb[(jo_max_err+1)*(IMO+1)+io_max_err+1] +ocnlonb[(jo_max_err+1)*(IMO+1)+io_max_err ])/4, (ocnlatb[ jo_max_err *(IMO+1)+io_max_err ] +ocnlatb[ jo_max_err *(IMO+1)+io_max_err+1] +ocnlatb[(jo_max_err+1)*(IMO+1)+io_max_err+1] +ocnlatb[(jo_max_err+1)*(IMO+1)+io_max_err ])/4); } if(land_same_as_atmos == 0) { if( npes > 1 ) { /* need to sum over all the pe and minus (n-1)*ocn_area */ sum_over_pe( ocn_minus_runoff, IMO*JMO ); for(i=0; i0.5) { double err = ocn_minus_runoff[i]/ocn_area[i]; double err2 = err*err; if (err2>o_max_err*o_max_err) { o_max_err = -err; jo_max_err = i/IMO; io_max_err = i-jo_max_err*IMO; } o_err += err2; o_cnt++; } } } if( pe == root_pe ) { printf( "MAX Ocean cell tiling error (land/ocean runoff) = %lg ", o_max_err); printf( "at (i=%d,j=%d), (x=%6.2f,y=%6.2f)\n", io_max_err+1, jo_max_err+1, (ocnlonb[ jo_max_err *(IMO+1)+io_max_err ] +ocnlonb[ jo_max_err *(IMO+1)+io_max_err+1] +ocnlonb[(jo_max_err+1)*(IMO+1)+io_max_err+1] +ocnlonb[(jo_max_err+1)*(IMO+1)+io_max_err ])/4, (ocnlatb[ jo_max_err *(IMO+1)+io_max_err ] +ocnlatb[ jo_max_err *(IMO+1)+io_max_err+1] +ocnlatb[(jo_max_err+1)*(IMO+1)+io_max_err+1] +ocnlatb[(jo_max_err+1)*(IMO+1)+io_max_err ])/4); } } } if( pe == root_pe ) { printf("\nGlobe tiling error = %lg\n", tot_area_AXL+tot_area_AXO-1); printf("\nMax # x-cell vertices = %d\n", n_max); } if( pe == root_pe ) { char cmd[BUFSIZ]; int ncid; int area_AO_id2, ia_AO_id2, ja_AO_id2, io_AO_id2, jo_AO_id2; int di_AO_id2, dj_AO_id2, di_AL_id2, dj_AL_id2, di_LO_id2, dj_LO_id2; int area_AL_id2, ia_AL_id2, ja_AL_id2, il_AL_id2, jl_AL_id2; int area_LO_id2, il_LO_id2, jl_LO_id2, io_LO_id2, jo_LO_id2; int dim_xba, dim_yba, dim_xbl, dim_ybl, xba_id, yba_id, xbl_id, ybl_id; int dim_xta, dim_yta, xta_id, yta_id, atm_area_id, a_dims[2]; int dim_xtl, dim_ytl, xtl_id, ytl_id, lnd_area_id, l_dims[2]; int lnd_cell_area_id, dim_null, null_id; int dim_xto, dim_yto, xto_id, yto_id, ocn_area_id, o_dims[2]; size_t i; if (ofile && nofile == 1) { sprintf(cmd, "cp %s %s", ofile, gfile); system(cmd); NC_CALL(nc_open(gfile, NC_WRITE, &ncid)) NC_CALL(nc_redef(ncid)) } else { NC_CALL(nc_create(gfile, NC_NOCLOBBER, &ncid)) } dim_null = -1; if (ofile) { if(iAXO ==0) { NC_CALL(nc_def_dim(ncid, "null_dimension", iAXO, &dim_null)) NC_CALL(nc_def_var(ncid, "null_dimension", NC_DOUBLE,1,&dim_null, &null_id)) dim_AO = dim_null; } else { NC_CALL(nc_def_dim(ncid, "i_atmXocn", iAXO, &dim_AO)) } NC_CALL(nc_def_var(ncid, "AREA_ATMxOCN",NC_DOUBLE,1,&dim_AO,&area_AO_id2)) NC_CALL(nc_def_var(ncid, "DI_ATMxOCN",NC_DOUBLE,1,&dim_AO,&di_AO_id2)) NC_CALL(nc_def_var(ncid, "DJ_ATMxOCN",NC_DOUBLE,1,&dim_AO,&dj_AO_id2)) NC_CALL(nc_def_var(ncid, "I_ATM_ATMxOCN", NC_INT, 1, &dim_AO, &ia_AO_id2)) NC_CALL(nc_def_var(ncid, "J_ATM_ATMxOCN", NC_INT, 1, &dim_AO, &ja_AO_id2)) NC_CALL(nc_def_var(ncid, "I_OCN_ATMxOCN", NC_INT, 1, &dim_AO, &io_AO_id2)) NC_CALL(nc_def_var(ncid, "J_OCN_ATMxOCN", NC_INT, 1, &dim_AO, &jo_AO_id2)) } if(iAXL ==0) { if(dim_null == -1) { NC_CALL(nc_def_dim(ncid, "null_dimension", iAXL, &dim_null)) NC_CALL(nc_def_var(ncid, "null_dimension", NC_DOUBLE,1,&dim_null, &null_id)) dim_AL = dim_null; } else { dim_AL = dim_null; } } else { NC_CALL(nc_def_dim(ncid, "i_atmXlnd", iAXL, &dim_AL)) } NC_CALL(nc_def_var(ncid, "AREA_ATMxLND", NC_DOUBLE,1,&dim_AL, &area_AL_id2)) NC_CALL(nc_def_var(ncid, "DI_ATMxLND", NC_DOUBLE,1,&dim_AL, &di_AL_id2)) NC_CALL(nc_def_var(ncid, "DJ_ATMxLND", NC_DOUBLE,1,&dim_AL, &dj_AL_id2)) NC_CALL(nc_def_var(ncid, "I_ATM_ATMxLND", NC_INT, 1, &dim_AL, &ia_AL_id2)) NC_CALL(nc_def_var(ncid, "J_ATM_ATMxLND", NC_INT, 1, &dim_AL, &ja_AL_id2)) NC_CALL(nc_def_var(ncid, "I_LND_ATMxLND", NC_INT, 1, &dim_AL, &il_AL_id2)) NC_CALL(nc_def_var(ncid, "J_LND_ATMxLND", NC_INT, 1, &dim_AL, &jl_AL_id2)) if(iLXO ==0) { if(dim_null == -1) { NC_CALL(nc_def_dim(ncid, "null_dimension", iLXO, &dim_null)) NC_CALL(nc_def_var(ncid, "null_dimension", NC_DOUBLE,1,&dim_null, &null_id)) dim_LO = dim_null; } else { dim_LO = dim_null; } } else { NC_CALL(nc_def_dim(ncid, "i_lndXocn", iLXO, &dim_LO)) } NC_CALL(nc_def_var(ncid, "AREA_LNDxOCN", NC_DOUBLE,1,&dim_LO, &area_LO_id2)) NC_CALL(nc_def_var(ncid, "DI_LNDxOCN", NC_DOUBLE,1,&dim_LO, &di_LO_id2)) NC_CALL(nc_def_var(ncid, "DJ_LNDxOCN", NC_DOUBLE,1,&dim_LO, &dj_LO_id2)) NC_CALL(nc_def_var(ncid, "I_LND_LNDxOCN", NC_INT, 1, &dim_LO, &il_LO_id2)) NC_CALL(nc_def_var(ncid, "J_LND_LNDxOCN", NC_INT, 1, &dim_LO, &jl_LO_id2)) NC_CALL(nc_def_var(ncid, "I_OCN_LNDxOCN", NC_INT, 1, &dim_LO, &io_LO_id2)) NC_CALL(nc_def_var(ncid, "J_OCN_LNDxOCN", NC_INT, 1, &dim_LO, &jo_LO_id2)) NC_CALL(nc_def_dim(ncid, "xba", IMA+1, &dim_xba)) NC_CALL(nc_def_var(ncid, "xba", NC_DOUBLE, 1, &dim_xba, &xba_id)) NC_CALL(nc_def_dim(ncid, "yba", JMA+1, &dim_yba)) NC_CALL(nc_def_var(ncid, "yba", NC_DOUBLE, 1, &dim_yba, &yba_id)) NC_CALL(nc_def_dim(ncid, "xta", IMA, &dim_xta)) NC_CALL(nc_def_var(ncid, "xta", NC_DOUBLE, 1, &dim_xta, &xta_id)) NC_CALL(nc_def_dim(ncid, "yta", JMA, &dim_yta)) NC_CALL(nc_def_var(ncid, "yta", NC_DOUBLE, 1, &dim_yta, &yta_id)) a_dims[1] = dim_xta; a_dims[0] = dim_yta; NC_CALL(nc_def_var(ncid, "AREA_ATM", NC_DOUBLE, 2, a_dims, &atm_area_id)) NC_CALL(nc_def_dim(ncid, "xbl", IML+1, &dim_xbl)) NC_CALL(nc_def_var(ncid, "xbl", NC_DOUBLE, 1, &dim_xbl, &xbl_id)) NC_CALL(nc_def_dim(ncid, "ybl", JML+1, &dim_ybl)) NC_CALL(nc_def_var(ncid, "ybl", NC_DOUBLE, 1, &dim_ybl, &ybl_id)) NC_CALL(nc_def_dim(ncid, "xtl", IML, &dim_xtl)) NC_CALL(nc_def_var(ncid, "xtl", NC_DOUBLE, 1, &dim_xtl, &xtl_id)) NC_CALL(nc_def_dim(ncid, "ytl", JML, &dim_ytl)) NC_CALL(nc_def_var(ncid, "ytl", NC_DOUBLE, 1, &dim_ytl, &ytl_id)) l_dims[1] = dim_xtl; l_dims[0] = dim_ytl; NC_CALL(nc_def_var(ncid, "AREA_LND", NC_DOUBLE, 2, l_dims, &lnd_area_id)) NC_CALL(nc_def_var(ncid, "AREA_LND_CELL", NC_DOUBLE, 2, l_dims, &lnd_cell_area_id)) if (ofile) { NC_CALL(nc_def_dim(ncid, "xto", IMO, &dim_xto)) NC_CALL(nc_def_var(ncid, "xto", NC_DOUBLE, 1, &dim_xto, &xto_id)) NC_CALL(nc_def_dim(ncid, "yto", JMO, &dim_yto)) NC_CALL(nc_def_var(ncid, "yto", NC_DOUBLE, 1, &dim_yto, &yto_id)) o_dims[1] = dim_xto; o_dims[0] = dim_yto; NC_CALL(nc_def_var(ncid, "AREA_OCN", NC_DOUBLE, 2, o_dims, &ocn_area_id)) } NC_CALL(nc_enddef(ncid)) if (ofile) { nc_copy_1d_double(ncidAO, area_AO_id, ncid, area_AO_id2, (int) iAXO); nc_copy_1d_double(ncidAO, di_AO_id, ncid, di_AO_id2, (int) iAXO); nc_copy_1d_double(ncidAO, dj_AO_id, ncid, dj_AO_id2, (int) iAXO); nc_copy_1d_int(ncidAO, ia_AO_id, ncid, ia_AO_id2, (int) iAXO); nc_copy_1d_int(ncidAO, ja_AO_id, ncid, ja_AO_id2, (int) iAXO); nc_copy_1d_int(ncidAO, io_AO_id, ncid, io_AO_id2, (int) iAXO); nc_copy_1d_int(ncidAO, jo_AO_id, ncid, jo_AO_id2, (int) iAXO); NC_CALL(nc_close(ncidAO)) } nc_copy_1d_double(ncidAL, area_AL_id, ncid, area_AL_id2, (int) iAXL); nc_copy_1d_double(ncidAL, di_AL_id, ncid, di_AL_id2, (int) iAXL); nc_copy_1d_double(ncidAL, dj_AL_id, ncid, dj_AL_id2, (int) iAXL); nc_copy_1d_int(ncidAL, ia_AL_id, ncid, ia_AL_id2, (int) iAXL); nc_copy_1d_int(ncidAL, ja_AL_id, ncid, ja_AL_id2, (int) iAXL); nc_copy_1d_int(ncidAL, il_AL_id, ncid, il_AL_id2, (int) iAXL); nc_copy_1d_int(ncidAL, jl_AL_id, ncid, jl_AL_id2, (int) iAXL); NC_CALL(nc_close(ncidAL)) if (ofile) { nc_copy_1d_double(ncidLO, area_LO_id, ncid, area_LO_id2, (int) iLXO); nc_copy_1d_double(ncidLO, di_LO_id, ncid, di_LO_id2, (int) iLXO); nc_copy_1d_double(ncidLO, dj_LO_id, ncid, dj_LO_id2, (int) iLXO); nc_copy_1d_int(ncidLO, il_LO_id, ncid, il_LO_id2, (int) iLXO); nc_copy_1d_int(ncidLO, jl_LO_id, ncid, jl_LO_id2, (int) iLXO); nc_copy_1d_int(ncidLO, io_LO_id, ncid, io_LO_id2, (int) iLXO); nc_copy_1d_int(ncidLO, jo_LO_id, ncid, jo_LO_id2, (int) iLXO); NC_CALL(nc_close(ncidLO)) } for (i=0;i<=IMA;i++) NC_CALL(nc_put_var1_double(ncid, xba_id, &i, &(atmlonb[i]))); for (i=0;i<=JMA;i++) NC_CALL(nc_put_var1_double(ncid, yba_id, &i, &(atmlatb[i]))); for (i=0;i