/* Routine to output a series of river long profiles past stations /* AND FIX THEM /* /puppet/esri/amls/proffix2.aml /* Based on /puppet/esri/amls/proffix2.aml v. 1.0, which was built as /* /flow/topog/areas/mexico/geo/proffix2.aml /* This version last used on /puppet/topog/alison and Noah's ASTER data. /* 6/5/98: add misc fields /* This version uses an input file of id,x,y for stations. It snaps them to the grid /* version .92 supports distance in meters /* When verion .93 comes to a "lake" filled in on the DEM, it looks ahead /* to the end of the lake (or maxsteps steps). It then proceds to follow /* the path downhill into the lake. /* When it reaches a sink, it looks to the nearest point on the lakesurface /* path. It also looks to the end of the lakesurface path. It can follow /* the lowest cell that lies between those bearings. /* This version has to look at the filled dem /* Calls: realdist.aml , bearing.aml , mdist.aml /* .98 bearend is up to 6 cells away /* .99 bearend is tuned to 4 cells, overtop override and error file added /* 1.00 gets the z increments right. There's a little wandering on flat lakes. /* 1.01 traces upstream first, cleans up some code /* 1.02 allows extra arguments for other coverages to query. /* 1.03 takes command line arguments, and recurses, distances up from mouth=0 /* changed realdist calls to mdist, is this okay? /* 1.04 Has hardwired minimum flowacc value to trace upstream to. /* If a variable %.i% exists, append to dem, filled,flowdir,flowacc /* and prepend to profiles. /* All the dialogs can be replaced with globals. /* 1.05 This applies to ID numbers for recursive calls. /* A new argument is idcov. If this defualts to NONE, the /* program will work in the old way, prepending 9s and other /* digits to the id number to form the id of a tributary profile. /* This is good for a ew tribs and for many mouths. /* However, you can now set this to the name of a calculated /* stream network. The id number of the stream segment will be used. /* This version last used on /puppet/topog/alison and Noah's ASTER data. /* 1.06 Made worldepsilon a function of cellsize. /* e.g. &r proffix2 elv filled flowdir flowacc mouth.xy 99999 _prf.xy NONE NONE CALCSTREAMS &args dem filled flowd flowak stations count outfile misc1 misc2 idcov &type proffix2.aml, version 1.05 minslope=0 &if [var .i] &then &type The value %.i% will appended to input grids and prepended to output. /* the following could be arguments &s dz = .000 /* minimum slope /* &s .xcellworld 85.9306 /* 85.9353 on diagonal (if undefined, it uses DEM units) /* &s .ycellworld 91.8154 /* 91.8154449 on diagonal /* The next three variables control finding paths through sinks. &s deltab 1 /* deltabearing: for rounding errors &s leeway 46 /* degrees ( > 45 call be dangerous) &s baseweight 100 /* ,lower is more center weight, try 100 - 200 &if [show program] ne 'GRID' &then ; GRID &if [length %outfile%] gt 0 &then &GOTO gotparams /*Input DEM, flowdirections grid, stations, and output. &LABEL getinput &do &until [query 'Continue?' .TRUE.] &if [var .dem] &then &s dem %.dem% &else &s dem = [response 'Input DEM name (dem)' dem] &if [locase %dem%] eq 'quit' &then; &RETURN &if [var .filled] &then &s filled %.filled% &else &s filled = [response 'Input FILLED DEM name (filled)' filled] &if [var .flowd] &then &s flowd %.flowd% &else &s flowd = [response 'Input flowdirection grid (FLOWDIR)' FLOWDIR] &if [var .flowak] &then &s flowak %.flowak% &else &s flowak = [response 'Input flowaccumulation grid (FLOWACC)' FLOWACC] &if [var .stations] &then &s stations %.stations% &else ; &do &TY input the name of a file or id,x,y,zmin &s stations = [response 'Input station file (pour.xy)' pour.xy] &end &if [var .count] &then &s count %.count% &else &s count = [response '# of pixels in downstream (99999)' 99999] &if [var .outfile] &then &s outfile %.outfile% &else &s outfile = [response 'Outfile name appendix? (_prf.xy)' _prf.xy] &if [var .misc1] &then &s misc1 %.misc1% &else &s misc1 = [response 'First extra data field? (NONE)' NONE] &if [var .misc2] &then &s misc2 %.misc2% &else &s misc2 = [response 'First extra data field? (NONE)' NONE] &if [var .idcov] &then &s idcov %.idcov% &else &s idcov = [response 'Cover from which to extract recursion IDs? (NONE)' NONE] &if [var .i] &then ; &do &s dem %dem%%.i% &s filled %filled%%.i% &s flowd %flowd%%.i% &s flowak %flowak%%.i% &s stations %stations%%.i% &s outfile _%.i%%outfile% &end &if [var .dem] &then &goto outofinteractiveloop &end &label outofinteractiveloop &s .dot ',' /* This will be changed with recursive calls &s .counter 1 &LABEL gotparams &type Entered values are %dem%, %filled%, %flowd%, %flowak%, %stations%, %count%, and %outfile% and %misc1% and %misc2% %idcov% /* set up conversion from arc/info 2^n flowdirection to Cartesian angle &s dir2bearing1 0 ; &s dir2bearing2 -45 &s dir2bearing4 -90 ; &s dir2bearing8 -135 &s dir2bearing16 180 ; &s dir2bearing32 135 &s dir2bearing64 90 ; &s dir2bearing128 45 &s minflow 0 /* Do only cells with flows > this number &s subthresh 200000 /* do every basin with more cells (6 hits on taiwan) &s subthresh 100000 /* do every basin with more cells setwindow %dem% &s sp ' ' &s m1 ''; &s m2 '' /* in case misc1 and misc2 are NONE &describe %dem% /* input information for dem &s dx = %GRD$DX% &s dy = %GRD$DY% &if [var .xcellworld] AND [var .ycellworld] &then ; &do &ty Using dx and dy of %.xcellworld% and %.ycellworld% &s dx %.xcellworld% &s dy %.ycellworld% &end &s worldepsilon [calc %dx% * 0.01] /* finally fixed 4/2/2004 &s testdist [calc %dx% * 1.01] /* if closer, two cells are adjacent. &s xmax = %GRD$XMAX% &s xmin = %GRD$XMIN% &s ymax = %GRD$YMAX% &s ymin = %GRD$YMIN% &s xminp = %xmin% + %dx% / 2.0 &s yminp = %ymin% + %dy% / 2.0 &type Beginning main program /* begin main loop &sys cat /dev/null >! fixer.gen /* if id,x,y,z was enterd instead of a file name &ty processing station %stations% &if %stations% cn ',' &then ; &do &ty Doing a single station &s idxyworld [unquote %stations%] &s status1 99 /* flag this as a command-line single-station run &GOTO dostation &end /* Read in *.txt Data Files &s lun [open %stations% status1 -read] &if %status1% ne 0 &then &return ERROR opening %stations% : %status1% &LABEL stationloop &s idxyworld [read %lun% status2 ] &if %status2% ne 0 &then; &do &ty end of file read %status2% &goto done &end &LABEL dostation &TYP parsing %idxyworld% &s station [extract 1 %idxyworld%] &s x0 [extract 2 %idxyworld%] &s y0 [extract 3 %idxyworld%] &s zmin [extract 4 %idxyworld%] &if ^ [var zmin] &then ; &s zmin 0.0 &if [type %zmin%] ge 0 &then ; &s zmin 0.0 &type Processing station : %station% at %x0% , %y0% &s xsta = [calc [calc [round [calc [calc %x0% - %xminp% ] / %dx% ] ] * %dx% ] + %xminp% ] &s ysta = [calc [calc [round [calc [calc %y0% - %yminp% ] / %dy% ] ] * %dy% ] + %yminp% ] &if [exists %station%%outfile% -file] &then; &type [delete %station%%outfile%] &sys touch %station%%outfile% /*prepare file and coords for inner loops &if [exists %station%%outfile%.err -file] &then; &type [delete %station%%outfile%.err] &sys touch %station%%outfile%.err /* for overtop warnings &if [type %.dot%] eq 1 &then /* if , &s dist 0.0 /* yes, now mouth = 0 &else &s dist %.distance% /* trib junc location on parent river &s x = %xsta% &s y = %ysta% &s dd = [sqrt [calc %dx% * %dx% + %dy% * %dy%]] /* THE UPSTREAM LOOP------------------------------------------------- &LABEL upstreamloop &s ak = [show cellvalue %flowak% %x% %y% ] /* 1.04 &if %ak% < %minflow% &then &goto continue &s c1 = [show cellvalue %flowak% [calc %x% + %dx%] %y%] &s c2 = [show cellvalue %flowak% [calc %x% + %dx%] [calc %y% - %dy%]] &s c3 = [show cellvalue %flowak% %x% [calc %y% - %dy%]] &s c4 = [show cellvalue %flowak% [calc %x% - %dx%] [calc %y% - %dy%]] &s c5 = [show cellvalue %flowak% [calc %x% - %dx%] %y%] &s c6 = [show cellvalue %flowak% [calc %x% - %dx%] [calc %y% + %dy%]] &s c7 = [show cellvalue %flowak% %x% [calc %y% + %dy%]] &s c8 = [show cellvalue %flowak% [calc %x% + %dx%] [calc %y% + %dy%]] /* reset cells that don't flow into current cell &if [type %c1%] gt 0 &then;&s c1 0;&else &if [show cellvalue %flowd% [calc %x% + %dx%] %y%] ne 16 &then;&s c1 = 0 &if [type %c2%] gt 0 &then;&s c2 0;&else &if [show cellvalue %flowd% [calc %x% + %dx%] [calc %y% - %dy%]] ne 32 &then;&s c2 = 0 &if [type %c3%] gt 0 &then;&s c3 0;&else &if [show cellvalue %flowd% %x% [calc %y% - %dy%]] ne 64 &then;&s c3 = 0 &if [type %c4%] gt 0 &then;&s c4 0;&else &if [show cellvalue %flowd% [calc %x% - %dx%] [calc %y% - %dy%]] ne 128 &then;&s c4 = 0 &if [type %c5%] gt 0 &then;&s c5 0;&else &if [show cellvalue %flowd% [calc %x% - %dx%] %y%] ne 1 &then;&s c5 = 0 &if [type %c6%] gt 0 &then;&s c6 0;&else &if [show cellvalue %flowd% [calc %x% - %dx%] [calc %y% + %dy%]] ne 2 &then;&s c6 = 0 &if [type %c7%] gt 0 &then;&s c7 0;&else &if [show cellvalue %flowd% %x% [calc %y% + %dy%]] ne 4 &then;&s c7 = 0 &if [type %c8%] gt 0 &then;&s c8 0;&else &if [show cellvalue %flowd% [calc %x% + %dx%] [calc %y% + %dy%]] ne 8 &then;&s c8 = 0 /* order flowaks and select largest &s ak2 = [extract 8 [sort %c1% %c2% %c3% %c4% %c5% %c6% %c7% %c8% -NUMERIC]] &s ak3 = [extract 7 [sort %c1% %c2% %c3% %c4% %c5% %c6% %c7% %c8% -NUMERIC]] &if %ak3% > %subthresh% &then; &do &select %ak3% /* determine which cell this is, flowdirectionwise &when %c1% &s us = 1 &when %c2% &s us = 2 &when %c3% &s us = 4 &when %c4% &s us = 8 &when %c5% &s us = 16 &when %c6% &s us = 32 &when %c7% &s us = 64 &when %c8% &s us = 128 &otherwise &ret ERROR finding start of trib /* error trap &end &if %us% in {1,2,128} &then ; &s xtrib [calc %x% + %dx%] &else &if %us% in {8,16,32} &then ; &s xtrib [calc %x% - %dx%] &else &s xtrib %x% &if %us% in {32,64,128} &then ; &s ytrib [calc %y% + %dy%] &else &if %us% in {2,4,8} &then ; &s ytrib [calc %y% - %dy%] &else &s ytrib %y% &s ztrib [show cellvalue %dem% %xtrib% %ytrib%] &if %us% in { 1, 16 } &then &s .distance [calc %dist% + %dx%] &if %us% in {4, 64} &then &s .distance [calc %dist% + %dy%] &if %us% in {2, 8, 32, 128} &then &s .distance [calc %dist% + %dd%] &if [type %.dot%] eq 1 &then ; &do /* , &s .dot 1 &if %idcov% eq NONE &then /* autogenerate ids for tribs &s nextid 9%.counter%9%station% &else ; &do mape [calc %xtrib% - %dx%] [calc %ytrib% - %dx%]~ [calc %xtrib% + %dx%] [calc %ytrib% + %dx%] /* search radius depends on mapex &ty Mapextent is [show mapex] clearsel %idcov% arc resel %idcov% arc one %xtrib% %ytrib% &s nextid %.dot%_[show select %idcov% line 1 item %idcov%#] &TYPE nextid at %xtrib% %ytrib% is _%nextid%_ &end &end &else; &do &s .dot [calc %.dot% + 1] &if %idcov% eq NONE &then /* autogenerate ids for tribs &s nextid 9%.counter%[substr %station% 2] &else ; &do mape [calc %xtrib% - %dx%] [calc %ytrib% - %dx%]~ [calc %xtrib% + %dx%] [calc %ytrib% + %dx%] /* search radius depends on mapex clearsel %idcov% arc resel %idcov% arc one %xtrib% %ytrib% &s nextid %.dot%_[show select %idcov% line 1 item %idcov%#] &TYPE nextid at %xtrib% %ytrib% is -%nextid%- &end &end &ty RECURSING from %x% %y% with dist of %.distance% on %station%/%nextid% &if ^ [exist prolog -file] &then ; &sys touch prolog &sys echo recursing from %x% %y% with dist of %.distance% on %station%/%nextid% >> prolog &s .counter [calc %.counter% + 1] &r proffix2 %dem% %filled% %flowd% %flowak% [quote %nextid%,%xtrib%,%ytrib%,%ztrib%] %count% %outfile% %misc1% %misc2% %idcov% &ty returning to top level from level %.dot% &if %.dot% eq 1 &then &s .dot , &else &s .dot [calc %.dot% - 1] &end /* ak3... make recursive call &select %ak2% /* determine which cell this is, flowdirectionwise &when 0 &s us 0 /* cannot have goto within select block &when %c1% &s us = 1 &when %c2% &s us = 2 &when %c3% &s us = 4 &when %c4% &s us = 8 &when %c5% &s us = 16 &when %c6% &s us = 32 &when %c7% &s us = 64 &when %c8% &s us = 128 &otherwise &ret /* error trap &end &if %us% = 0 &then; &GOTO continue /* &TY goin up %x% %y% %z% %ak2% &if %us% in { 1, 2, 128 } &then &s x = [calc %x% + %dx% ] &if %us% in {8, 16, 32} &then &s x = [calc %x% - %dx% ] &if %us% in {32, 64, 128} &then &s y = [calc %y% + %dy% ] &if %us% in {2, 4, 8} &then &s y = [calc %y% - %dy% ] &if %us% in { 1, 16 } &then &s dist [calc %dist% + %dx%] &if %us% in {4, 64} &then &s dist [calc %dist% + %dy%] &if %us% in {2, 8, 32, 128} &then &s dist [calc %dist% + %dd%] &if %x% > %xmax% or %x% < %xmin% &then &GOTO continue /* needed? &if %y% > %ymax% or %y% < %ymin% &then &GOTO continue /* needed? &if %ak2% = 0 &then &GOTO continue /* needed? &sys echo -n %.dot% /* violation flag /* &if %zf% < %zft% &then &sys echo -n x /* &s zft = %zf% &GOTO upstreamloop &label continue &s z = [show cellvalue %dem% %x% %y% ] &type Found top stream point %x% %y% %z% at distance %dist% /* THE DOWNSTREAM LOOP------------------------------------------------- &s xold -99 ; &s yold -99 &s xoldold -99 ; &s yoldold -99 /* to prevent loops &s pt = 0 &s zold = 99999 &s olddir 1 &s overtopok .false. &set inlake .FALSE. &s z = [show cellvalue %dem% %x% %y% ] &s fixed %z% &s ak = [show cellvalue %flowak% %x% %y% ] /* trouble if we start in a lake &s zold = %fixed% &type /& ---> downstream &if [type %.dot%] eq 1 &then /* mainstem &s needextra .false. &else &s needextra .true. /* let it get an extra cell to connect subbasins &LABEL cellloop &if %inlake% &then ; &do /* ignore flowdir entirely &do i = 1 &to %k% &ty lakepoint %i% [value lx%i%] [value ly%i%] &end &LABEL nextlakepoint &run mdist %x% %y% [value lx%k%] [value ly%k%] &TY distance to end of lake is %.dist% &if %.dist% lt .0001 &then ; &do /* we got to the last cell of the lake &s inlake .FALSE. &goto notinlake &end &s closest 9999 &s closest_to_end 9999 &do i = 1 &to %k% &run mdist %x% %y% [value lx%i%] [value ly%i%] &if [calc %.dist% - %closest%] < 0.0001 &then ; &do /* if closer, more or less &s closest %.dist% &s closek %i% &run mdist [value lx%i%] [value ly%i%] [value lx%k%] [value ly%k%] &s closest_to_end %.dist% &end &end &ty closest is %closek% %closest% %k% %x% %y% [value lx%k%] [value ly%k%] &if [calc %k% - %closek%] > 4 &then /* dont miss a turn on a long lake &run bearing %x% %y% [value lx[calc %closek% + 4]] [value ly[calc %closek% + 4]] &else &run bearing %x% %y% [value lx%k%] [value ly%k%] &s bearend %.bearing% &if %closest% < .0001 &then /*on closest point, look downstream &s bearclose [value dir2bearing[show cellvalue %flowd% %x% %y%]] &else ; &do &run bearing %x% %y% [value lx%closek%] [value ly%closek%] &s bearclose %.bearing% &end &s brange = %bearend% - %bearclose% &if %brange% < 0 &then ; &s brange = [calc %brange% + 360.0] &if %brange% < 180 &then ; &do /* %bearend% to %bearclose% is small &s right = %bearclose% - %deltab% &s left = %bearend% + %deltab% &end &else ; &do &s left = %bearclose% + %deltab% &s right = %bearend% - %deltab% &end &if %bearend% < %right% &then &s bearend [calc %bearend% + 360.] &if %left% < %bearend% &then &s left [calc %left% + 360.] &ty bearclose is %bearclose%, range is %brange%,ler are %left% %bearend% %right% &s br %right% - %leeway% &s bl %left% + %leeway% &s ibr = [trunc [calc %br% / 45]] &if %br% > 0 &then ; &s ibr [calc %ibr% + 1 ] &s ibl = [trunc [calc %bl% / 45]] &if %bl% < 0 &then ; &s ibl [calc %ibl% - 1 ] /* &ty left = %bl% %ibl% , right = %br% %ibr% &label checkneighbors /* jump back here for overtop override &s lowestslope 999999 &do j = %ibr% &to %ibl% &s jdeg [calc %j% * 45] &if %jdeg% < %bearend% &then &s weight = [calc %baseweight% + %jdeg% - %bearend%] &else &if %jdeg% > %bearend% &then &s weight = [calc %baseweight% - %jdeg% + %bearend%] &else &s weight = %baseweight% &if %weight% < 20 &then &s weight 20 /* &if %jdeg% < %right% &then /* &s weight = [calc %baseweight% + %jdeg% - %right%] /* &else &if %jdeg% > %left% &then /* &s weight = [calc %baseweight% - %jdeg% + %left%] /* &else /* &s weight = %baseweight% /*&type %j% weight %weight% %jdeg% %right% %left% /* very verbose &if %j% in {7,0,1,-15,-8,-7,-1,8,9,17} &then &s xt = [calc %x% + %dx% ] &else &if %j% in {3,4,5,-5,-4,-3,11,12,13} &then &s xt = [calc %x% - %dx% ] &else &if %j% in {2,6,10,14,-6,-2} &then &s xt %x% &else &RETURN ERROR, Odd value in ibr or ibl &if %j% in {1,2,3,-7,-6,-5,9,10,11} &then &s yt = [calc %y% + %dy% ] &else &if %j% in {5,6,7,-3,-2,-1,13,14,15} &then &s yt = [calc %y% - %dy% ] &else &if %j% in {0,4,8,12,-8,-4} &then &s yt %y% &else &RETURN ERROR, Odddd value in ibr or ibl &if %j% in {-8,-4,0,4,8,12} &then /* to the downward cell &s dxyworld %dx% &else &if %j% in {-6,-2,2,6,10,14} &then /* to the downward cell &s dxyworld %dy% &else &s dxyworld %dd% &s zt = [show cellvalue %dem% %xt% %yt% ] /* We want to stay in the low part of the lake without wandering in circles. /* Let's try comparing to %z% instead of the %fixed% /* but if z==t, it will always beat out a climb, i.e. if you descend into /* a gorge, you will climb out the other side, but if you are skating /* across a plain the is flat in the unfilled DEM, you are likely to /* be traveling in circles. So, on a plain, we a) Recalc the weight, /* b) add a little slope /* but this can let you avoid a level path to go uphill, then force you /* downhill to reach that cell next time. /* fixed <= z <= filled /* &type xt yt zt = %xt% %yt% %zt% /* very verbose &if %zt% < %fixed% &then /* dropping (negative slope) &s testslope [calc ( %zt% - %z% ) / %dxyworld%] &else &if %zt% < %z% &then /* dropping in a local context &s testslope [calc ( %zt% - %z% ) / %dxyworld%] /* negative &else &if %zt% eq %z% &then &do /* flat spots can be dangerous &s testslope [calc 0.5 / %dxyworld%] /* give it a little slope and &s weight [calc %weight% + ( %weight% - 78 ) * 2.0] /* give more emphasis to the bearing &if %weight% < 0 &then &s weight [calc -1.0 / %weight%] /* never change the sign of the weight &ty revised weight is %weight% &end &else /* when climbing, go to lowest cell, ignoring diag factor &s testslope [calc %zt% - %z% ] &if %testslope% <= 0 &then &s testslope [calc %testslope% * %weight% ] &else &s testslope [calc %testslope% / %weight% ] /* downhill is neg slope /*&ty testslope = %testslope% %zt% %z% %dxyworld% /* very verbose &if %testslope% < %lowestslope% &then ; &do &r mdist %xold% %yold% %xt% %yt% /* &type .distance %.dist% %xold% %yold% %xt% %yt% /* very verbose &if %.dist% > %worldepsilon% &then ; &do /* do not double back &r mdist %xoldold% %yoldold% %xt% %yt% &if %.dist% > %worldepsilon% &then ; &do /* do not loop back &if %zt% le [show cellvalue %filled% %x% %y%] or %overtopok% &then ; &do &s lowestslope %testslope% &s lowestx %xt% &s lowesty %yt% &s lowestz %zt% /* &s lowestj %j% &end &else &ty overtop protection invoked at %x% %y% &end &else &ty loopback protection invoked at %x% %y% &end &else &ty doubleback protection invoked at %x% %y% &end &end /* j &s overtopok .false. &if %lowestslope% eq 999999 &then ; &do &ty overtop override invoked at %x% %y% &sys echo overtop override invoked at %x% %y% >> %station%%outfile%.err &s overtopok .true. &goto checkneighbors /* &RETURN ERROR No lake cells were picked &end &s xoldold %xold% ; &s yoldold %yold% &s xold %x% ; &s yold %y% &s x %lowestx% &s y %lowesty% &s z %lowestz% &s dist [calc %dist% - %dxyworld%] &TY Picked lake cell jxyzZold %j% %x% %y% %z% %Zold% %ibr% %ibl% &s ak = [show cellvalue %flowak% %x% %y% ] &if %z% ge %zold% &then &s fixed [calc %zold% - ( %dz% * %dxyworld% ) ] &else &s fixed = %z% &if %misc1% ne 'NONE' &then ; &s m1 = [show cellvalue %misc1% %x% %y%] &if %misc2% ne 'NONE' &then ; &s m2 = [show cellvalue %misc2% %x% %y%] &sys echo %dist% %fixed% %z% %ak% %x% %y% %m1% %m2% >> %station%%outfile% &sys echo 1 %x% %y% %fixed% >> fixer.gen &s zold = %fixed% &s pt = [calc %pt% + 1] &if %z% eq [show cellvalue %filled% %x% %y%] &then ;&do &s inlake .FALSE. &goto notinlake &end &if %inlake% &then &goto nextlakepoint &end /* inlake &LABEL notinlake &s ds = [show cellvalue %flowd% %x% %y% ] &if [quote %ds%] eq NODATA &then ; &goto breakloop &s xt %x%; &s yt %y% &if %ds% in { 1, 2, 128 } &then; &s xt = [calc %x% + %dx% ] &if %ds% in {8,16,32} &then; &s xt = [calc %x% - %dx% ] &if %ds% in {32, 64, 128} &then; &s yt = [calc %y% + %dy% ] &if %ds% in {2, 4, 8} &then; &s yt = [calc %y% - %dy% ] &if %xt% > %xmax% or %xt% < %xmin% &then; &goto breakloop &if %yt% > %ymax% or %yt% < %ymin% &then; &goto breakloop &s zt = [show cellvalue %dem% %xt% %yt% ] &if [type %zt%] ge 0 &then; &goto breakloop &if ( %zt% < %zmin% ) &then; &do &if %needextra% &then ; &s needextra .false. /* get extra cell &else ; &goto breakloop &end /* &ty %x% %y% %xt% %yt% %ds% %zt% %fixed% &if %zt% < [show cellvalue %filled% %xt% %yt%] &then ; &do /* entering lake &TYPE PREVIEW THE DEFAULT LAKE PATH ----------------------------- &s k 0 &s lakelev [show cellvalue %filled% %xt% %yt%] &if %lakelev% in {-99} &then ; &do /* &SYSTEM drip &PAUSE Lake is %lakelev% &end &LABEL lakeloop &s fillt = [show cellvalue %filled% %xt% %yt% ] /* &ty fillt= %fillt% , lakelev = %lakelev% %xt% %yt% &if %fillt% < %lakelev% &then ; &goto breaklake /* out of this sink &s k [calc %k% + 1] &s lx%k% %xt% &s ly%k% %yt% /* &ty K: %k% [value lx%k%] [value ly%k%] %z% %zt% %fixed% &s ds = [show cellvalue %flowd% [value lx%k%] [value ly%k%]] &if [quote %ds%] eq NODATA &then ; &goto breaklake &s xt [value lx%k%]; &s yt [value ly%k%] &if %ds% in { 1, 2, 128 } &then; &s xt = [calc [value lx%k%] + %dx% ] &if %ds% in {8, 16,32} &then; &s xt = [calc [value lx%k%] - %dx% ] &if %ds% in {32, 64, 128} &then; &s yt = [calc [value ly%k%] + %dy% ] &if %ds% in {2, 4, 8} &then; &s yt = [calc [value ly%k%] - %dy% ] &if %xt% > %xmax% or %xt% < %xmin% &then; &goto breaklake &if %yt% > %ymax% or %yt% < %ymin% &then; &goto breaklake &s zt = [show cellvalue %dem% %xt% %yt% ] &goto lakeloop &LABEL breaklake /* &s k [calc %k% - 1] &if %k% eq 1 &then &GOTO regular /* We aim for pourpoint, so one-point sinks ok &s inlake .TRUE. &goto cellloop &end /* lakeloop &LABEL regular &s z %zt% &s x %xt% &s y %yt% &s ak [show cellvalue %flowak% %x% %y% ] &if %z% ge %zold% &then &s fixed [calc %zold% - ( %dz% * %dxyworld% ) ] &else &s fixed = %z% &s zold %fixed% &if %ds% in {1,16} &then /* to the downward cell &s dxyworld %dx% &else &if %ds% in {4,64} &then /* to the downward cell &s dxyworld %dy% &else &s dxyworld %dd% &s dist [calc %dist% - %dxyworld%] &if %misc1% ne 'NONE' &then ; &s m1 = [show cellvalue %misc1% %x% %y%] &if %misc2% ne 'NONE' &then ; &s m2 = [show cellvalue %misc2% %x% %y%] &sys echo %dist% %fixed% %z% %ak% %x% %y% %m1% %m2% >> %station%%outfile% &s pt = [calc %pt% + 1] &sys echo -n . &if %pt% < %count% &then ; &GOTO cellloop &label breakloop &if %status1% > 0 &then ; &goto done &type Profile for %station% completed, on to the next... &GOTO stationloop /* end of station loop &LABEL done &type That is all folks!