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