share/landread.c

   1 #ifndef CRAY
   2 # ifdef NOUNDERSCORE
   3 #      define GET_TERRAIN get_terrain
   4 #      define GET_LANDUSE get_landuse
   5 # else
   6 #   ifdef F2CSTYLE
   7 #      define GET_TERRAIN get_terrain__
   8 #      define GET_LANDUSE get_landuse__
   9 #   else
  10 #      define GET_TERRAIN get_terrain_
  11 #      define GET_LANDUSE get_landuse_
  12 #   endif
  13 # endif
  14 #endif
  15
  16 #ifdef LANDREAD_STUB
  17 #ifndef MS_SUA
  18 # include <stdio.h>
  19 #endif
  20
  21 int GET_TERRAIN (        float *adx,
  22                          float *xlat,
  23                          float *xlon,
  24                          float       *terrain,
  25                          int   *mix,
  26                          int   *mjx,
  27                          int   *iyyn,
  28                          int   *jxxn,
  29                          int   *ipath , int * ipathlen)  /* integer coded ASCII string from Funtran and len */
  30
  31 {
  32 #ifndef MS_SUA
  33  fprintf(stderr, "***************************************************************\n" ) ;
  34  fprintf(stderr, "Access to RSMAS Topo Ingest Code is by Special Arrangement Only\n" ) ;
  35  fprintf(stderr, "in WRF.  Please contact wrfhelp@ucar.edu.                      \n" ) ;
  36  fprintf(stderr, "***************************************************************\n" ) ;
  37 #endif
  38  return(0) ;
  39 }
  40
  41 int GET_LANDUSE (        float *adx,
  42                          float *xlat,
  43                          float *xlon,
  44                          float       *landuse,
  45                          int   *mix,
  46                          int   *mjx,
  47                          int   *iyyn,
  48                          int   *jxxn,
  49                          int   *ipath , int * ipathlen ) /* integer coded ASCII string from Funtran and len */
  50
  51 {
  52 #ifndef MS_SUA
  53  fprintf(stderr, "***************************************************************\n" ) ;
  54  fprintf(stderr, "Access to RSMAS Topo Ingest Code is by Special Arrangement Only\n" ) ;
  55  fprintf(stderr, "in WRF.  Please contact wrfhelp@ucar.edu.                      \n" ) ;
  56  fprintf(stderr, "***************************************************************\n" ) ;
  57 #endif
  58  return(0) ;
  59 }
  60 #else
  61
  62 #ifdef FSEEKO_OK
  63 #  define _FILE_OFFSET_BITS 64
  64 #endif
  65 #ifndef MS_SUA
  66 # include <stdio.h>
  67 #endif
  68 #if (RPC_TYPES == 1)
  69 # include <rpc/types.h>
  70 # include <rpc/xdr.h>
  71 #elif (RPC_TYPES == 2)
  72 # include <tirpc/rpc/types.h>
  73 # include <tirpc/rpc/xdr.h>
  74 #endif
  75 #include <math.h>
  76 #ifndef MACOS
  77 # include <malloc.h>
  78 #else
  79 # include <malloc/malloc.h>
  80 #endif
  81 #include <unistd.h>
  82 #include <string.h>
  83 #define MAXTOPOFILES  100
  84 #define MAXLEN        4096
  85
  86
  87 typedef struct
  88 {
  89   /* Filenames. */
  90   char  fn[MAXTOPOFILES][MAXLEN];
  91
  92   /* Grid spacings in km. */
  93   float dx[MAXTOPOFILES];
  94
  95   /* Number of entries. */
  96   int num;
  97 } TsFileInfo;
  98
  99 static TsFileInfo tsfTopo;
 100 static TsFileInfo tsfOcean;
 101 static TsFileInfo tsfLU;
 102
 103 static int tsFileInfo_initialized = 0;
 104 /* static float last_adx = 0.0 ; */
 105
 106 static char tsfTopo_fn[MAXLEN];
 107 static char tsfLU_fn[MAXLEN];
 108 static char tsfOcean_fn[MAXLEN];
 109
 110 static float vmiss;
 111
 112 static int    numHeaderBytes;
 113 static int    globalNx;
 114 static int    globalNy;
 115 static int    tileNx;
 116 static int    tileNy;
 117 static int    extraNx;
 118 static int    extraNy;
 119 static int    numTilesX;
 120 static int    numTilesY;
 121 static double dlat;
 122 static double dlon;
 123 static double lat0;
 124 static double lon0;
 125 static int    ntiles;
 126 static int    wrapx;
 127 static int    wrapy;
 128
 129 /* File information. */
 130 static XDR  *xdrs;
 131 static FILE *fp;
 132
 133 int nint(double x)
 134 {
 135   if ( x > 0.0 ) { return( (int)(x + 0.5) ) ; }
 136   return((int)(x - 0.5));
 137 }
 138
 139 double aint(double x)
 140 {
 141   int ix = (int)(x);
 142   return((double)(ix));
 143 }
 144
 145 double anint(double x)
 146 {
 147   if (x > 0.0) return((double)((int)(x + 0.5)));
 148   return((double)((int)(x - 0.5)));
 149 }
 150
 151 static double normalizeAngle(double ang)
 152 {
 153   for (;;)
 154     {
 155       if (ang >= 360.0)
 156         {
 157           ang -= 360.0;
 158         }
 159       else if (ang < 0.0)
 160         {
 161           ang += 360.0;
 162         }
 163       else
 164         {
 165           break;
 166         }
 167     }
 168
 169   return(ang);
 170 }
 171
 172 static double lonDistNowrap(double lon1, double lon2)
 173 {
 174   double lon11 = normalizeAngle(lon1);
 175   double lon22 = normalizeAngle(lon2);
 176   if (lon22 < lon11) lon22 += 360.0;
 177   return(fabs(lon22 - lon11));
 178 }
 179
 180 int tsLatLonToGridpoint(double  lat,
 181                         double  lon,
 182                         double       *ix,
 183                         double       *iy)
 184 {
 185   *ix = lonDistNowrap(lon0, lon) / dlon;
 186   *iy = (lat - lat0) / dlat;
 187   return(0);
 188 }
 189
 190 static int areEqual(double v1, double v2)
 191 {
 192   if (fabs(v1-v2) < 0.001) return(1);
 193   return(0);
 194 }
 195
 196 static int setWrapAroundFlags(void)
 197 {
 198   /* Compute the end gridpoint location in x. */
 199   double lon1  = lon0 + dlon*(globalNx);
 200   double lon2  = lon0 + dlon*(globalNx-1);
 201   double lat1  = lat0 + dlat*(globalNy);
 202   double lon0n = normalizeAngle(lon0);
 203   double lon1n = normalizeAngle(lon1);
 204   double lon2n = normalizeAngle(lon2);
 205
 206   wrapx = 0;
 207   if (areEqual(lon0n, lon1n))
 208     {
 209       /* Here the first and last indices in x are one grid interval
 210          apart. */
 211       wrapx = 1;
 212     }
 213   else if (areEqual(lon0n, lon2n))
 214     {
 215       /* Here the first and last indices in x are coincident. */
 216       wrapx = 2;
 217     }
 218
 219   wrapy = 0;
 220   if (areEqual(lat0, -90.0))
 221     {
 222       /* Here the first and last indices in x are one grid interval
 223          apart. */
 224       wrapy += 1;
 225     }
 226   if (areEqual(lat1, 90.0))
 227     {
 228       /* Here the first and last indices in x are coincident. */
 229       wrapy += 2;
 230     }
 231
 232   return(0);
 233 }
 234
 235 static int isMissing(float v)
 236 {
 237   if (fabs(vmiss - v) < 0.1) return(1);
 238   return(0);
 239 }
 240
 241 float tsGetValueInt(int aix, int aiy)
 242 {
 243   float f = vmiss;
 244
 245   int iy = aiy;
 246   int ix = aix;
 247
 248   /* Perform bounds checking. */
 249   if (iy < 0)
 250     {
 251       return(f);
 252     }
 253   else if (iy > globalNy - 1)
 254     {
 255       return(f);
 256     }
 257
 258   if (aix < 0)
 259     {
 260       if (wrapx == 1)
 261         {
 262           int n  = -(aix - (globalNx - 1)) / globalNx;
 263           ix += n*globalNx;
 264         }
 265       else if (wrapx == 2)
 266         {
 267           int nx = globalNx - 1;
 268           int n  = -(aix - (nx - 1)) / nx;
 269           ix += n*nx;
 270         }
 271       else
 272         {
 273           return(f);
 274         }
 275     }
 276
 277   if (ix > globalNx-1)
 278     {
 279       if (wrapx == 1)
 280         {
 281           int n  = aix / globalNx;
 282           ix -= n*globalNx;
 283         }
 284       else if (wrapx == 2)
 285         {
 286           int nx = globalNx - 1;
 287           int n  = aix / nx;
 288           ix -= n*nx;
 289         }
 290       else
 291         {
 292           return(f);
 293         }
 294     }
 295
 296   int tx  = ix / tileNx;
 297   int ty  = iy / tileNy;
 298   int tn  = tx + ty*numTilesX;
 299   int txg = ix - tx*tileNx;
 300   int tyg = iy - ty*tileNy;
 301   int gn  = txg + tyg*tileNx;
 302
 303   long long ll_gn = gn;
 304   long long ll_numHeaderBytes  = numHeaderBytes;
 305   long long ll_tileNx = tileNx;
 306   long long ll_tileNy = tileNy;
 307
 308 #ifdef FSEEKO64_OK
 309   /* This is used on machines that support fseeko64. Tested for in ./configure script */
 310   long long loc = ll_numHeaderBytes + ll_tileNx*ll_tileNy*sizeof(float)*tn +
 311     ll_gn*sizeof(float);
 312
 313   /* Seek to the proper location in the file and get the data value. */
 314   fseeko64(fp, loc, SEEK_SET);
 315 #else
 316 #  ifdef FSEEKO_OK
 317   /* This is used on machines that support _FILE_OFFSET_BITS=64 which makes
 318      off_t be 64 bits, and for which fseeko can handle 64 bit offsets.  This
 319      is tested in the ./configure script */
 320   off_t loc = ll_numHeaderBytes + ll_tileNx*ll_tileNy*sizeof(float)*tn +
 321     ll_gn*sizeof(float);
 322
 323   fseeko(fp, loc, SEEK_SET);
 324 #  else
 325   /* Note, this will not work correctly for very high resolution terrain input
 326      because the offset is only 32 bits.   */
 327   off_t loc = ll_numHeaderBytes + ll_tileNx*ll_tileNy*sizeof(float)*tn +
 328     ll_gn*sizeof(float);
 329
 330   fseek(fp, loc, SEEK_SET);
 331 #  endif
 332 #endif
 333   xdr_float(xdrs, (float *) &f);
 334
 335   return(f);
 336 }
 337
 338 float tsGetValue(double ix, double iy)
 339 {
 340   int i0 = (int)(floor(ix));
 341   int j0 = (int)(floor(iy));
 342   int i1 = (int)(ceil(ix));
 343   int j1 = (int)(ceil(iy));
 344
 345   /* Interpolate linearly to (oiloc, ojloc). */
 346   float v0 = tsGetValueInt(i0,j0);
 347   float v1 = tsGetValueInt(i0,j1);
 348   float v2 = tsGetValueInt(i1,j0);
 349   float v3 = tsGetValueInt(i1,j1);
 350
 351   if (isMissing(v0)) return(vmiss);
 352   if (isMissing(v1)) return(vmiss);
 353   if (isMissing(v2)) return(vmiss);
 354   if (isMissing(v3)) return(vmiss);
 355
 356   double w0 = ix - i0;
 357   double w1 = iy - j0;
 358
 359   float v4 = v2*w0 + v0*(1.0-w0);
 360   float v5 = v3*w0 + v1*(1.0-w0);
 361   float v6 = w1*v5 + (1.0-w1)*v4;
 362   float val = v6;
 363
 364   return(val);
 365 }
 366
 367 float tsGetValueLatLon(double lat, double lon)
 368 {
 369   double ix, iy;
 370   tsLatLonToGridpoint(lat,lon,&ix,&iy);
 371   return(tsGetValue(ix,iy));
 372 }
 373
 374 int tsCloseTileSet(void)
 375 {
 376   if (xdrs)
 377     {
 378       xdr_destroy(xdrs);
 379       free(xdrs);
 380       xdrs = NULL;
 381     }
 382
 383   if (fp)
 384     {
 385       fclose(fp);
 386       fp = NULL;
 387     }
 388
 389   return(0);
 390 }
 391
 392 int tsInitTileSet(char *fn)
 393 {
 394   vmiss = -100000000.00;
 395
 396   xdrs = NULL;
 397   fp   = NULL;
 398
 399 # if 0
 400   fprintf(stderr,"Open %s\n", fn) ;
 401 # endif
 402
 403   /* fp = (FILE *) fopen64(fn, "r"); */
 404   if (( fp = fopen(fn, "r")) == NULL ) {
 405 #ifndef MS_SUA
 406     fprintf(stderr,"tsInitTileSet: can not open %s\n",fn) ;
 407 #endif
 408     return(1) ;
 409   }
 410   xdrs = (XDR *) malloc(sizeof(XDR));
 411   xdrstdio_create(xdrs, fp, XDR_DECODE);
 412
 413   numHeaderBytes = 5000;
 414
 415   xdr_int(xdrs,    (int *)    &globalNx);
 416   xdr_int(xdrs,    (int *)    &globalNy);
 417   xdr_int(xdrs,    (int *)    &tileNx);
 418   xdr_int(xdrs,    (int *)    &tileNy);
 419   xdr_int(xdrs,    (int *)    &extraNx);
 420   xdr_int(xdrs,    (int *)    &extraNy);
 421   xdr_int(xdrs,    (int *)    &numTilesX);
 422   xdr_int(xdrs,    (int *)    &numTilesY);
 423   xdr_double(xdrs, (double *) &dlat);
 424   xdr_double(xdrs, (double *) &dlon);
 425   xdr_double(xdrs, (double *) &lat0);
 426   xdr_double(xdrs, (double *) &lon0);
 427   xdr_int(xdrs,    (int *)    &ntiles);
 428
 429   setWrapAroundFlags();
 430
 431   return(0);
 432 }
 433
 434 int tsPrintTileSetInfo()
 435 {
 436   return(0);
 437 }
 438
 439 int tsInitFileInfo (char path[])
 440 {
 441   int i, n;
 442   char type[MAXLEN];
 443   char  res[MAXLEN];
 444   char   fn[MAXLEN];
 445   char buff[MAXLEN];
 446   float dx;
 447
 448   tsfTopo.num  = 0;
 449   tsfOcean.num = 0;
 450   tsfLU.num    = 0;
 451
 452   tsfLU.dx[0] = 0.;
 453   tsfTopo.dx[0] = 0.;
 454   tsfOcean.dx[0] = 0. ;
 455
 456   if (access("RSMAS_Topo_Land.TBL", F_OK) == 0 ) {
 457   /* Read in the list of topography/land use filenames. */
 458     fp = fopen("RSMAS_Topo_Land.TBL", "r");
 459     if ( fp == NULL ) {
 460 #ifndef MS_SUA
 461       fprintf(stderr, "tsInitFileInfo : can not open RSMAS_Topo_Land.TBL\n");
 462 #endif
 463       return(-1);
 464     }
 465
 466     /* skipps header */
 467     fgets(buff, MAXLEN, fp);
 468
 469     while (fscanf(fp, "%s %s %s", type, res, fn) != EOF) {
 470         sscanf(res, "%f", &dx);
 471         if (strcmp(type, "landuse") == 0)
 472           {
 473             if ( tsfLU.num >= MAXTOPOFILES ) {continue;}
 474             n = tsfLU.num ;
 475             for ( i = 0 ; i < tsfLU.num ; i++ ) {
 476               if ( tsfLU.dx[i] > dx ) {
 477                 n = i ;
 478                 break;
 479               }
 480             }
 481             for ( i = tsfLU.num ; i > n ; i-- ) {
 482                tsfLU.dx[i]=tsfLU.dx[i-1];
 483                strcpy(tsfLU.fn[i], tsfLU.fn[i-1]);
 484             }
 485             tsfLU.dx[n] = dx;
 486             strcpy(tsfLU.fn[n], fn);
 487             tsfLU.num++;
 488           }
 489         else if (strcmp(type, "topography") == 0)
 490           {
 491             if ( tsfTopo.num >= MAXTOPOFILES ) {continue;}
 492             n = tsfTopo.num;
 493             for ( i = 0 ; i < tsfTopo.num ; i++ ) {
 494               if ( tsfTopo.dx[i] > dx ) {
 495                 n = i ;
 496                 break;
 497               }
 498             }
 499             for ( i = tsfTopo.num ; i > n ; i-- ) {
 500                tsfTopo.dx[i]=tsfTopo.dx[i-1];
 501                strcpy(tsfTopo.fn[i], tsfTopo.fn[i-1]);
 502             }
 503             tsfTopo.dx[n] = dx;
 504             strcpy(tsfTopo.fn[n], fn);
 505             tsfTopo.num++;
 506           }
 507         else if (strcmp(type, "bathymetry") == 0)
 508           {
 509             if ( tsfOcean.num >= MAXTOPOFILES ) {continue;}
 510             n = tsfOcean.num;
 511             for ( i = 0 ; i < tsfOcean.num ; i++ ) {
 512               if ( tsfOcean.dx[i] > dx ) {
 513                 n = i ;
 514                 break;
 515               }
 516             }
 517             for ( i = tsfOcean.num ; i > n ; i-- ) {
 518                tsfOcean.dx[i]=tsfOcean.dx[i-1];
 519                strcpy(tsfOcean.fn[i], tsfOcean.fn[i-1]);
 520             }
 521             tsfOcean.dx[n] = dx;
 522             strcpy(tsfOcean.fn[n], fn);
 523             tsfOcean.num++;
 524           }
 525     }
 526
 527     fclose(fp);
 528
 529   } else {
 530
 531     for ( i = 1; i < 10 ; i++ ) {
 532       sprintf(tsfTopo.fn[tsfTopo.num], "%s/topo.%02dkm.ts", path, i);
 533       if (access(tsfTopo.fn[tsfTopo.num], F_OK) == 0) { tsfTopo.dx[tsfTopo.num] = i; tsfTopo.num++ ; }
 534     }
 535     for ( i = 10; i<=40 ; i += 10 ) {
 536       sprintf(tsfTopo.fn[tsfTopo.num], "%s/topo.%02dkm.ts", path, i);
 537       if (access(tsfTopo.fn[tsfTopo.num], F_OK) == 0) { tsfTopo.dx[tsfTopo.num] = i; tsfTopo.num++ ; }
 538     }
 539
 540     for ( i = 1; i < 10 ; i++ ) {
 541       sprintf(tsfLU.fn[tsfLU.num], "%s/glcc.usgs20.%02dkm.ts", path, i);
 542       if (access(tsfLU.fn[tsfLU.num], F_OK) == 0) { tsfLU.dx[tsfLU.num] = i; tsfLU.num++ ; }
 543     }
 544     for ( i = 10; i<=40 ; i += 10 ) {
 545       sprintf(tsfLU.fn[tsfLU.num], "%s/glcc.usgs20.%02dkm.ts", path, i);
 546       if (access(tsfLU.fn[tsfLU.num], F_OK) == 0) { tsfLU.dx[tsfLU.num] = i; tsfLU.num++ ; }
 547     }
 548
 549     for ( i = 1; i < 10 ; i++ ) {
 550       sprintf(tsfOcean.fn[tsfOcean.num], "%s/tbase.%02dkm.ts", path, i);
 551       if (access(tsfOcean.fn[tsfOcean.num], F_OK) == 0) { tsfOcean.dx[tsfLU.num] = i; tsfOcean.num++ ; }
 552     }
 553     for ( i = 10; i<=40 ; i += 10 ) {
 554       sprintf(tsfOcean.fn[tsfOcean.num], "%s/tbase.%02dkm.ts", path, i);
 555       if (access(tsfOcean.fn[tsfOcean.num], F_OK) == 0) { tsfOcean.dx[tsfOcean.num] = i; tsfOcean.num++ ; }
 556     }
 557   }
 558 # if 0
 559   for ( i = 0 ; i < tsfTopo.num ; i++ ) {
 560     fprintf(stderr,"%02d. %s\n",i, tsfTopo.fn[i]) ;
 561   }
 562   for ( i =0 ; i < tsfLU.num ; i++ ) {
 563     fprintf(stderr,"%02d. %s\n",i, tsfLU.fn[i]) ;
 564   }
 565   for ( i =0 ; i < tsfOcean.num ; i++ ) {
 566     fprintf(stderr,"%02d. %s\n",i, tsfOcean.fn[i]) ;
 567   }
 568 # endif
 569
 570   return(1);
 571
 572 }
 573
 574 int GET_LANDUSE (        float *adx,
 575                          float *xlat,
 576                          float *xlon,
 577                          float       *landuse,
 578                          int   *mix,
 579                          int   *mjx,
 580                          int   *iyyn,
 581                          int   *jxxn,
 582                          int   *ipath , int * ipathlen ) /* integer coded ASCII string from Funtran and len */
 583 {
 584   int i, j ;
 585   char path[256];
 586   int ix,iy, offset ;
 587   float lat, lon, tv;
 588   double fix, fiy;
 589
 590   static float last_adx = 0.0 ;
 591
 592   if ( tsFileInfo_initialized == 0 ) {
 593      for (i = 0 ; i < *ipathlen ; i++ ) {
 594        path[i] = ipath[i] ;
 595      }
 596      path[*ipathlen] = '\0' ;
 597      tsFileInfo_initialized = tsInitFileInfo(path);
 598      if ( tsFileInfo_initialized == -1 ) { return(1); }
 599   }
 600
 601 /* if ( fabs(last_adx - *adx)  > 1.0E-6 ) { */
 602 /*   last_adx = *adx; */
 603     if ( tsfLU.num > 0 ) {
 604       strcpy(tsfLU_fn, tsfLU.fn[tsfLU.num-1]);
 605       for ( i = 0; i < tsfLU.num; i++) {
 606 # if 0
 607           fprintf(stderr,"%d fn %s dx %f adx %f\n",i,tsfLU.fn[i],tsfLU.dx[i],*adx) ;
 608 # endif
 609           if (tsfLU.dx[i] > *adx) {
 610               strcpy(tsfLU_fn, tsfLU.fn[i-1]);
 611               break;
 612             }
 613         }
 614     } else {
 615 # ifndef MS_SUA
 616         fprintf(stderr, "Not found LANDUSE datasets!\n");
 617 # endif
 618         return(1);
 619     }
 620 /* } */
 621
 622   /* Get the land use. */
 623   if (tsInitTileSet(tsfLU_fn)) { return(1); }
 624
 625   for ( j = 0; j < *jxxn; j++) {
 626     offset = *mix*j;
 627     for ( i = 0; i < *iyyn; i++) {
 628       lat = xlat[offset + i];
 629       lon = xlon[offset + i];
 630
 631       tsLatLonToGridpoint(lat,lon,&fix,&fiy);
 632       ix = nint(fix);
 633       iy = nint(fiy);
 634       tv = tsGetValueInt(ix, iy);
 635
 636       /* Set out-of-range values to water. */
 637       if (tv < 0.9 || tv > 24.1) tv = 16.0;
 638
 639       landuse[offset + i] = tv;
 640     }
 641   }
 642   tsCloseTileSet();
 643 }
 644
 645 int GET_TERRAIN (        float *adx,
 646                          float *xlat,
 647                          float *xlon,
 648                          float       *terrain,
 649                          int   *mix,
 650                          int   *mjx,
 651                          int   *iyyn,
 652                          int   *jxxn,
 653                          int   *ipath , int * ipathlen)  /* integer coded ASCII string from Funtran and len */
 654 {
 655   int i, j ;
 656   char path[256];
 657   int ix, iy, offset ;
 658   float lon, lat, tv;
 659   double fix, fiy;
 660
 661   static float last_adx = 0.0 ;
 662
 663   if ( tsFileInfo_initialized == 0 ) {
 664      for (i = 0 ; i < *ipathlen ; i++ ) {
 665        path[i] = ipath[i] ;
 666      }
 667      path[*ipathlen] = '\0' ;
 668      tsFileInfo_initialized = tsInitFileInfo(path);
 669      if ( tsFileInfo_initialized == -1 ) { return(1); }
 670   }
 671
 672   /* Use the data with the largest spacing less than the grid
 673      spacing specified in the argument list. */
 674 /* if ( fabs(last_adx - *adx)  > 1.0E-6 ) { */
 675 /*  last_adx = *adx; */
 676     if ( tsfTopo.num > 0 ) {
 677       strcpy(tsfTopo_fn, tsfTopo.fn[tsfTopo.num-1]);
 678       for (i = 0; i < tsfTopo.num; i++) {
 679 #if 0
 680 fprintf(stderr,"%d fn %s dx %f adx %f\n",i,tsfTopo.fn[i],tsfTopo.dx[i],*adx ) ;
 681 #endif
 682           if ( tsfTopo.dx[i] > *adx) {
 683               strcpy(tsfTopo_fn, tsfTopo.fn[i-1]);
 684               break;
 685             }
 686         }
 687     } else {
 688 #ifndef MS_SUA
 689       fprintf(stderr,"Not found GTOPO datasets\n");
 690 #endif
 691       return(1);
 692     }
 693
 694 #ifdef TERRAIN_TBASE
 695     if ( tsfOcean.num > 0 ) {
 696       strcpy(tsfOcean_fn, tsfOcean.fn[tsfOcean.num-1]);
 697       for ( i = 0; i < tsfOcean.num; i++) {
 698           if (tsfOcean.dx[i] > *adx) {
 699             strcpy(tsfOcean_fn, tsfOcean.fn[i-1]);
 700             break;
 701           }
 702       }
 703     } else {
 704 # ifndef MS_SUA
 705         fprintf(stderr, "Not found TBASE datasets!\n");
 706 # endif
 707         return(1);
 708     }
 709 #endif
 710 /* } */
 711
 712   /* First get the terrain from GTOPO30. */
 713     if (tsInitTileSet(tsfTopo_fn)) { return(1); }
 714
 715     for ( j = 0; j < *jxxn; j++)
 716       { offset = *mix*j;
 717         for ( i = 0; i < *iyyn; i++)
 718           {
 719             lat = xlat[offset + i];
 720             lon = xlon[offset + i];
 721
 722             tsLatLonToGridpoint(lat,lon,&fix,&fiy);
 723             tv = tsGetValue(fix, fiy);
 724             terrain[offset + i] = tv;
 725           }
 726       }
 727     tsCloseTileSet();
 728
 729 #ifdef TERRAIN_TBASE
 730   /* Next get the terrain from TBASE. */
 731     if (tsInitTileSet(tsfOcean_fn)) { return(1); }
 732
 733     for ( j = 0; j < *jxxn; j++)
 734       { offset = *mix*j;
 735         for ( i = 0; i < *iyyn; i++)
 736           {
 737             lat = xlat[offset + i];
 738             lon = xlon[offset + i];
 739
 740             tsLatLonToGridpoint(lat,lon,&fix,&fiy);
 741             tv = tsGetValue(fix, fiy);
 742             if (isMissing(terrain[offset+i]))
 743               {
 744                 if (tv < 0.0) tv = 0.0;
 745                 terrain[offset + i] = tv;
 746               }
 747           }
 748       }
 749     tsCloseTileSet();
 750 #endif
 751   return(0);
 752 }
 753
 754 #ifdef BATHYMETRY
 755 int get_bathymetry_(float *tadx,
 756                     float *xlat,
 757                     float *xlon,
 758                     float       *depth,
 759                     int   *mix,
 760                     int   *mjx,
 761                     int   *iyyn,
 762                     int   *jxxn,
 763                     float *mindepth,
 764                     float *zlimww3,
 765                     int   *ipath , int * ipathlen)  /* integer coded ASCII string from Funtran and len */
 766 {
 767   int i, j;
 768   char fn[MAXLEN];
 769   char path[1024];
 770   float maxdx, tv;
 771   double fix, fiy, lat, lon;
 772   int ix,iy,nx,ny;
 773   int offset;
 774
 775   nx=*mix;
 776   ny=*mjx;
 777
 778   /* Set grid resolution to .1 km to get highest resolution data possible. */
 779   float adx = 0.1;
 780
 781   if ( tsFileInfo_initialized == 0 ) {
 782      for (i = 0 ; i < *ipathlen ; i++ ) {
 783        path[i] = ipath[i] ;
 784      }
 785      path[*ipathlen] = '\0' ;
 786      tsFileInfo_initialized = tsInitFileInfo(path);
 787      if ( tsFileInfo_initialized == -1 ) { return(1); }
 788   }
 789
 790   /* Get the water depth from TBASE. */
 791   if ( tsfOcean.num > 0 ) {
 792     /* Use the data with highest resolution possible. */
 793     maxdx = 0.0;
 794     strcpy(fn, tsfOcean.fn[tsfOcean.num-1]);
 795     for ( i = 0; i < tsfOcean.num; i++)
 796       {
 797         if (tsfOcean.dx[i] < maxdx) continue;
 798         if ( tsfOcean.dx[i] < adx)
 799           {
 800             maxdx = tsfOcean.dx[i];
 801             strcpy(fn, tsfOcean.fn[i]);
 802           }
 803       }
 804
 805     if (tsInitTileSet(fn))
 806       {
 807         return(1);
 808       }
 809
 810     for ( i = 0; i < nx*ny; i++)
 811       {
 812         depth[i] = vmiss;
 813       }
 814
 815     for ( j = 0; j < *jxxn; j++)
 816       { offset = nx * j;
 817         for ( i = 0; i < *iyyn; i++)
 818           {
 819             lat = xlat[offset + i];
 820             lon = xlon[offset + i];
 821
 822             tsLatLonToGridpoint(lat,lon,&fix,&fiy);
 823             tv = tsGetValue(fix, fiy);
 824             if (isMissing(depth[offset+i]))
 825               {
 826                 depth[offset + i] = -tv;
 827               }
 828           }
 829       }
 830     tsCloseTileSet();
 831   } else {
 832 # ifndef MS_SUA
 833         fprintf(stderr, "Not found TBASE datasets!\n");
 834 # endif
 835         return(1);
 836   }
 837
 838   /* Next get the land use. */
 839   if ( tsfLU.num > 0 ) {
 840     /* Use the data with the largest spacing less than the grid
 841        spacing specified in the argument list. */
 842     maxdx = 0.0;
 843     strcpy(fn, tsfLU.fn[tsfLU.num-1]);
 844     for ( i = 0; i < tsfLU.num; i++)
 845       {
 846         if (tsfLU.dx[i] < maxdx) continue;
 847         if (tsfLU.dx[i] < adx)
 848           {
 849             maxdx = tsfLU.dx[i];
 850             strcpy(fn, tsfLU.fn[i]);
 851           }
 852       }
 853
 854     if (tsInitTileSet(fn))
 855       {
 856         return(1);
 857       }
 858
 859     for ( j = 0; j < *jxxn; j++)
 860       { offset = nx*j;
 861         for ( i = 0; i < *iyyn; i++)
 862           {
 863             lat = xlat[offset + i];
 864             lon = xlon[offset + i];
 865
 866             tsLatLonToGridpoint(lat,lon,&fix,&fiy);
 867             ix = nint(fix);
 868             iy = nint(fiy);
 869             tv = tsGetValueInt(ix, iy);
 870
 871             /* Set out-of-range values to water. */
 872             if (tv < 0.9 || tv > 24.1) tv = 16.0;
 873
 874             if (fabs(tv - 16.0) < 0.1)
 875               {
 876                 /* Water. */
 877                 if (1)
 878                   {
 879                     if (depth[offset + i] < *mindepth) depth[offset + i] = *mindepth;
 880                   }
 881                 else
 882                   {
 883                     if (depth[offset + i] < -(*zlimww3) )
 884                       {
 885                         /* Water depth below zlimww3, so turn this point
 886                            into land. */
 887                         depth[offset + i] = -0.1;
 888                       }
 889                     else if (depth[offset + i] < *mindepth)
 890                       {
 891                         depth[offset + i] = *mindepth;
 892                       }
 893                   }
 894               }
 895             else
 896               {
 897                 /* Land. Set depth to 0.0. */
 898                 depth[offset + i] = 0.0;
 899               }
 900           }
 901       }
 902     tsCloseTileSet();
 903   } else {
 904 # ifndef MS_SUA
 905         fprintf(stderr, "Not found LANDUSE datasets!\n");
 906 # endif
 907         return(1);
 908   }
 909   return(0);
 910 }
 911 #endif
 912
 913 #endif