|
ocssw
1.0
|
00001 #include <sys/types.h> 00002 #include <unistd.h> 00003 00004 /* ============================================================================ */ 00005 /* module sssref.c - retrieve salinity from WOA climatology */ 00006 /* ============================================================================ */ 00007 00008 #include "l12_proto.h" 00009 #include "hdf_utils.h" 00010 #include "h5io.h" 00011 00012 static float sssbad = BAD_FLT; 00013 static int32_t format = -1; 00014 00015 00016 /* ----------------------------------------------------------------------------------- */ 00017 /* get_woasssclim() - read and interpolate WOA salinity climatology */ 00018 /* ----------------------------------------------------------------------------------- */ 00019 00020 #define WOASSSNX 360 00021 #define WOASSSNY 180 00022 00023 float get_woasssclim(char *sssfile, float lon, float lat, int day) 00024 { 00025 static int firstCall = 1; 00026 static int nx = WOASSSNX; 00027 static int ny = WOASSSNY; 00028 static float dx = 360.0/WOASSSNX; 00029 static float dy = 180.0/WOASSSNY; 00030 static float sssref[WOASSSNY+2][WOASSSNX+2]; 00031 static float sss_sav[ WOASSSNY + 2 ][ WOASSSNX + 2 ]; 00032 00033 float sss = sssbad; 00034 int i,j,ii; 00035 float xx,yy; 00036 float t,u; 00037 00038 if (firstCall) { 00039 00040 float ssstmp[WOASSSNY][WOASSSNX]; 00041 char *tmp_str; 00042 char name [H4_MAX_NC_NAME] = ""; 00043 char sdsname[H4_MAX_NC_NAME] = ""; 00044 int32 sd_id; 00045 int32 sds_id; 00046 int32 rank; 00047 int32 nt; 00048 int32 dims[H4_MAX_VAR_DIMS]; 00049 int32 nattrs; 00050 int32 start[4] = {0,0,0,0}; 00051 int32 status, ix, iy, ct; 00052 int32 lymin, lxmin, lymax, lxmax; 00053 float sum; 00054 00055 int16 mon = (int) day/31 +1; // month of year (no need for perfection) 00056 00057 firstCall = 0; 00058 00059 printf("Loading SSS reference from Climatlogy file: %s\n",sssfile); 00060 printf("\n"); 00061 00062 /* Open the file */ 00063 sd_id = SDstart(sssfile, DFACC_RDONLY); 00064 if (sd_id == FAIL){ 00065 fprintf(stderr,"-E- %s line %d: error reading %s.\n", 00066 __FILE__,__LINE__,sssfile); 00067 exit(1); 00068 } 00069 00070 /* Get the SDS index */ 00071 sprintf(sdsname,"sss%2.2i",mon); 00072 sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname)); 00073 status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs); 00074 status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) &ssstmp[0][0]); 00075 if (status != 0) { 00076 fprintf(stderr,"-E- %s Line %d: Error reading SDS %s from %s.\n", 00077 __FILE__,__LINE__,sdsname,sssfile); 00078 exit(1); 00079 } 00080 status = SDendaccess(sds_id); 00081 status = SDend(sd_id); 00082 00083 /* rotate 180-deg and add wrapping border to simplify interpolation */ 00084 /* new grid is -180.5,180.5 in i=0,361 and -90.5,90.5 in j=0,181 */ 00085 00086 for (j=0; j<ny; j++) { 00087 for (i=0; i<nx; i++) { 00088 ii = (i < nx/2) ? i+nx/2 : i-nx/2; 00089 if (ssstmp[j][i] > 0) 00090 sssref[j+1][ii+1] = ssstmp[j][i]; 00091 else 00092 sssref[j+1][ii+1] = sssbad; 00093 } 00094 sssref[j+1][0] = sssref[j+1][nx]; 00095 sssref[j+1][nx+1] = sssref[j+1][1]; 00096 } 00097 for (i=0; i<nx+2; i++) { 00098 sssref[0] [i] = sssref[1][i]; 00099 sssref[ny+1][i] = sssref[ny][i]; 00100 } 00101 /* 00102 * Extend the grid by 1 point (use only full grid boxes later) 00103 */ 00104 memcpy( sss_sav, sssref, 00105 ( WOASSSNY + 2 ) * ( WOASSSNX + 2 ) * sizeof( float ) ); 00106 for( iy = 0; iy < ny+2; iy++ ) 00107 { 00108 lymin = ( iy == 0 ) ? 0 : iy - 1; 00109 lymax = ( iy == ny+1 ) ? ny+1 : iy + 1; 00110 00111 for( ix = 0; ix < nx+2; ix++ ) 00112 { 00113 if( sssref[iy][ix] < sssbad + 1 ) 00114 { 00115 lxmin = ( ix == 0 ) ? 0 : ix - 1; 00116 lxmax = ( ix == nx+1 ) ? nx+1 : ix + 1; 00117 00118 sum = 0.; 00119 ct = 0; 00120 for( j = lymin; j < lymax + 1; j++ ) 00121 for( i = lxmin; i < lxmax + 1; i++ ) 00122 { 00123 if( sss_sav[j][i] > sssbad + 1 ) 00124 { 00125 sum += sss_sav[j][i]; 00126 ct ++; 00127 } 00128 } 00129 if( ct > 0 ) 00130 sssref[iy][ix] = sum / ct; 00131 } 00132 } 00133 } 00134 } /* end 1-time grid set up */ 00135 00136 /* locate LL position within reference grid */ 00137 i = MAX(MIN((int) ((lon+180.0+dx/2)/dx),WOASSSNX+1),0); 00138 j = MAX(MIN((int) ((lat+ 90.0+dy/2)/dy),WOASSSNY+1),0); 00139 00140 /* compute longitude and latitude of that grid element */ 00141 xx = i*dx - 180.0 - dx/2; 00142 yy = j*dy - 90.0 - dy/2; 00143 00144 /* bilinearly interpolate, only using full grid boxes */ 00145 t = (lon - xx)/dx; 00146 u = (lat - yy)/dy; 00147 00148 if( ( sssref[j][i] > sssbad+1 ) && ( sssref[j][i+1] > sssbad+1 ) && 00149 ( sssref[j+1][i] > sssbad+1 ) && ( sssref[j+1][i+1] > sssbad+1 ) ) { 00150 00151 sss = (1-t)*(1-u) * sssref[j ][i ] 00152 + t*(1-u) * sssref[j ][i+1] 00153 + t*u * sssref[j+1][i+1] 00154 + (1-t)*u * sssref[j+1][i ]; 00155 00156 } else 00157 sss = sssbad; 00158 00159 return(sss); 00160 } 00161 00162 /* 00163 * for the HYCOM file usage... 00164 */ 00165 #define HYCOMNX 1440 00166 #define HYCOMNY 720 00167 00168 float get_hycom_sss( char *sssfile, float lon, float lat, float *sss ) 00169 /******************************************************************* 00170 00171 get_hycom_sss 00172 00173 purpose: read in the daiy HYCOM sss and give the interpolated value 00174 for the input lat, lon 00175 00176 Returns type: int - return status: 0 is good, 1 if any other error 00177 will fail with -1 status if an h5 io error found (not associated with 00178 the file not being h5 or correct format of h5) 00179 00180 Parameters: (in calling order) 00181 Type Name I/O Description 00182 ---- ---- --- ----------- 00183 char * sssfile I sss file name 00184 float lon I longitude: -180 -> 180 00185 float lat I latitude: -90 -> 90 00186 float * sss O sss found - will be sssbad 00187 if over land or otherwise 00188 bad value 00189 00190 Modification history: 00191 Programmer Date Description of change 00192 ---------- ---- --------------------- 00193 W. Robinson, SAIC 12 Oct 2012 adapted from get_woasssclim in this file 00194 W. Robinson, SAIC 16 Oct 2012 this version 2 will interpolate using 00195 0 for any missing data 00196 00197 Note that the HYCOM file only contains a dataset of salinity, so 00198 the size of 1440, 720 probably should be an ID requirement as is the 00199 existance of the dataset name 'salinity' 00200 00201 *******************************************************************/ 00202 { 00203 static int firstCall = 1; 00204 static int nx = HYCOMNX; 00205 static int ny = HYCOMNY; 00206 static float dx = 360.0 / HYCOMNX; 00207 static float dy = 180.0 / HYCOMNY; 00208 static float sssref[ HYCOMNY + 2 ][ HYCOMNX + 2 ]; 00209 static float sss_sav[ HYCOMNY + 2 ][ HYCOMNX + 2 ]; 00210 00211 int i, j, iter_ct, iter_max; 00212 float t, u, xx,yy; 00213 00214 *sss = sssbad; 00215 if( firstCall ) 00216 { 00217 float ssstmp[ny][nx], sum; 00218 h5io_str fid, dsid, grpid; 00219 int nobj, *types, ndim, dim_siz[10], sto_len, ix, iy, ct; 00220 int lymin, lxmin, lymax, lxmax; 00221 char **o_names; 00222 H5T_class_t class; 00223 hid_t native_typ; 00224 00225 firstCall = 0; 00226 00227 printf("Loading SSS reference from HYCOM file: %s\n",sssfile); 00228 printf("\n"); 00229 00230 /* 00231 * check the file to see that it is correct - just return to caller if 00232 * not right format, but exit if something that should be performed can't be 00233 */ 00234 if( h5io_openr( sssfile, 0, &fid ) != 0 ) 00235 return 1; 00236 /* the group may need setting to work, see if '/' will do */ 00237 if( h5io_set_grp( &fid, "/", &grpid ) != 0 ) 00238 { 00239 fprintf( stderr, "-E- %s, %d: Failed to set trivial group\n", 00240 __FILE__, __LINE__ ); 00241 exit(-1); 00242 } 00243 if( h5io_grp_contents( &grpid, &nobj, &o_names, &types ) != 0 ) 00244 { 00245 fprintf( stderr, "-E- %s, %d: Failed to find contents of sss file\n", 00246 __FILE__, __LINE__ ); 00247 exit(-1); 00248 } 00249 /* first name should be salinity */ 00250 if( strcmp( o_names[0], "salinity" ) != 0 ) 00251 return 1; 00252 00253 /* set to the salinity dataset and get info */ 00254 if( h5io_set_ds( &fid, "salinity", &dsid ) != 0 ) 00255 return 1; 00256 if( h5io_info( &dsid, NULL, &class, &native_typ, &ndim, dim_siz, 00257 &sto_len ) != 0 ) 00258 { 00259 fprintf( stderr, "-E- %s, %d: Failed to find info on salinity dataset\n", 00260 __FILE__, __LINE__ ); 00261 exit(-1); 00262 } 00263 00264 if( ( ndim != 2 ) || ( dim_siz[0] != 720 ) || ( dim_siz[1] != 1440 ) ) 00265 return 1; 00266 /* 00267 * at this point, we have identified the salinity file, just read 00268 * the dataset now 00269 */ 00270 if( h5io_rd_ds( &dsid, ( void *) ssstmp ) != 0 ) 00271 { 00272 fprintf( stderr, "-E- %s, %d: Failed to read salinity dataset\n", 00273 __FILE__, __LINE__ ); 00274 exit(-1); 00275 } 00276 /* close and go */ 00277 h5io_close( &dsid ); 00278 h5io_close( &fid ); 00279 00280 /* add wrapping border to simplify interpolation */ 00281 /* new grid is -180.125,180.125 in i=0,1441 and -90.125,90.125 in j=0,721 */ 00282 for( j = 0; j < ny; j++ ) 00283 { 00284 for( i = 0; i < nx; i++ ) 00285 { 00286 if( ( ssstmp[j][i] > 0) && ( ssstmp[j][i] <= 1000. ) ) 00287 sssref[j+1][i+1] = ssstmp[j][i]; 00288 else 00289 sssref[j+1][i+1] = sssbad; 00290 } 00291 sssref[j+1][0] = sssref[j+1][nx]; 00292 sssref[j+1][nx+1] = sssref[j+1][1]; 00293 } 00294 for (i=0; i<nx+2; i++) 00295 { 00296 sssref[0] [i] = sssref[1][i]; 00297 sssref[ny+1][i] = sssref[ny][i]; 00298 } 00299 /* 00300 * extend data to have good interpolation points at edge = coast 00301 */ 00302 iter_ct = 0; 00303 iter_max = 4; 00304 do 00305 { 00306 iter_ct++; 00307 memcpy( sss_sav, sssref, 00308 ( HYCOMNY + 2 ) * ( HYCOMNX + 2 ) * sizeof( float ) ); 00309 for( iy = 0; iy < ny+2; iy++ ) 00310 { 00311 lymin = ( iy == 0 ) ? 0 : iy - 1; 00312 lymax = ( iy == ny+1 ) ? ny+1 : iy + 1; 00313 00314 for( ix = 0; ix < nx+2; ix++ ) 00315 { 00316 if( sssref[iy][ix] < sssbad + 1 ) 00317 { 00318 lxmin = ( ix == 0 ) ? 0 : ix - 1; 00319 lxmax = ( ix == nx+1 ) ? nx+1 : ix + 1; 00320 00321 sum = 0.; 00322 ct = 0; 00323 for( j = lymin; j < lymax + 1; j++ ) 00324 for( i = lxmin; i < lxmax + 1; i++ ) 00325 { 00326 if( sss_sav[j][i] > sssbad + 1 ) 00327 { 00328 sum += sss_sav[j][i]; 00329 ct ++; 00330 } 00331 } 00332 if( ct > 0 ) 00333 sssref[iy][ix] = sum / ct; 00334 } 00335 } 00336 } 00337 } while ( iter_ct < iter_max ); 00338 /* and end set up */ 00339 } 00340 00341 /* locate LL position within reference grid */ 00342 i = MAX( MIN( (int) ( ( lon + 180.0 + dx / 2 ) / dx ), nx + 1 ), 0 ); 00343 j = MAX( MIN( (int) ( ( lat + 90.0 + dy / 2 ) / dy ), ny + 1 ), 0 ); 00344 00345 /* compute longitude and latitude of that grid element */ 00346 xx = i * dx - 180.0 - dx / 2; 00347 yy = j * dy - 90.0 - dy / 2; 00348 00349 /* bilinearly interpolate, replacing missing (land) values with 00350 average of valid values in box */ 00351 t = ( lon - xx ) / dx; 00352 u = ( lat - yy ) / dy; 00353 00354 if( ( sssref[j][i] > sssbad+1 ) && ( sssref[j][i+1] > sssbad+1 ) && 00355 ( sssref[j+1][i] > sssbad+1 ) && ( sssref[j+1][i+1] > sssbad+1 ) ) 00356 *sss = (1-t)*(1-u) * sssref[j ][i ] 00357 + t*(1-u) * sssref[j ][i+1] 00358 + t*u * sssref[j+1][i+1] 00359 + (1-t)*u * sssref[j+1][i ]; 00360 else 00361 *sss = sssbad; 00362 00363 return 0; 00364 } 00365 00366 /* ----------------------------------------------------------------------------------- */ 00367 /* get_sssref() - retrieves reference sea surface salinity . */ 00368 /* ----------------------------------------------------------------------------------- */ 00369 #define WOACLIM 1 00370 float get_sssref(char *sssfile, float lon, float lat, int day) 00371 { 00372 float sss = sssbad; 00373 int32_t sd_id; 00374 00375 if (format < 0) { 00376 00377 /* Does the file exist? */ 00378 if (access(sssfile, F_OK) || access(sssfile, R_OK)) { 00379 printf("-E- %s: SSS input file '%s' does not exist or cannot open.\n", 00380 __FILE__, sssfile); 00381 exit(1); 00382 } 00383 00384 /* What is it? */ 00385 sd_id = SDstart(sssfile, DFACC_RDONLY); 00386 if (sd_id != FAIL) { 00387 /* Format is HDF-like */ 00388 char title[255] = ""; 00389 if (SDreadattr(sd_id,SDfindattr(sd_id,"Title"),(VOIDP)title) == 0) { 00390 if (strstr(title,"WOA Sea Surface Salinity Monthly Climatology") != NULL) { 00391 format = WOACLIM; 00392 } 00393 } else { 00394 printf("-E- %s: unable to read SSS input file %s.\n",__FILE__, sssfile); 00395 exit(1); 00396 } 00397 SDend(sd_id); 00398 } 00399 else 00400 { 00401 if( get_hycom_sss( sssfile, lon, lat, &sss ) == 0 ) 00402 format = 2; 00403 } 00404 } 00405 00406 switch (format) { 00407 case WOACLIM: 00408 sss = get_woasssclim(sssfile,lon,lat,day); 00409 break; 00410 case 2: 00411 get_hycom_sss( sssfile, lon, lat, &sss ); 00412 break; 00413 default: 00414 printf("-E- %s: unknown SSS input file format for %s.\n",__FILE__, sssfile); 00415 exit(1); 00416 break; 00417 } 00418 00419 // assume fresh water if no data 00420 00421 if (sss < 0.0) { 00422 sss = 0.0; 00423 } 00424 00425 return(sss); 00426 } 00427 00428
1.7.6.1