INSTYTUT FIZYKI JĄDROWEJ INSTITUTE OF NUCLEAR PHYSJCS ИНСТИТУТ ЯДЕРНОЙ сризики

INSTYTUT FIZYKI JĄDROWEJ INSTITUTE OF NUCLEAR PHYSJCS ИНСТИТУТ ЯДЕРНОЙ срИЗИКИ

R A P O R T No 920/PM

PROGRAM DO OBLICZEŃ RUCHU ZANIECZYSZCZEŃ W WODACH GRUNTOWYCH

OPARTY NA ANALITYCZNYM ROZWIĄZANIU RÓWNANIA DYSPERSJI

P I O T R M A Ł O S Z E W S K I

KRAKÓW 1976

PROGRAM DO OBLICZEŃ RUCHU ZANIECZYSZCZEŃ W WODACH GRUNTOWYCH OPARU HA ANALITYCZNYM ROZWIĄZANIU

RÓWNANIA DYSPERSJI

A COMPUTER PROGRAM POR CALCULATION OP THE POLLUTANT MOVEMENT IN THE GROUNDWATER BASED ON ANALYTICAL

SOLUTION OF THE DISPERSION EQUATION

ПРОГРАММА ДЛЯ РАСЧЕТОВ ДВИЖЕНИЯ ЗАГРЯЗНЕНИЙ В ГРУНТОВЫХ ВОДАХ ОСНОВАННАЯ НА АНАЛИТИЧЕСКОМ

РЕШЕНИИ УРАВНЕНИЯ ДИСПЕРСИИ

Piotr Małoszewski

Instytut Fizyki Jądrowej, Kraków

M a j 1976

NAKŁADEM INSTYTUTU FIZYKI JĄDROWEJ W KRAKOWIE, UL. RADZIKOWSKIEGO 152

Kopię kserograficzną, druk i oprawę wykonano w IFJ Kraków Wydanie I Zam. 107/76 Nakład 70 egz.

Streszczenie.

W pracy przedstawiono opis programu na maszynę cyfrową, umożliwiającego obliczenie ruchu zanieczysz- czeń w wodach Rruntowych. Program został napisany w opar- ciu o znane rozwiązania równania dyspersji dla przypadków występowania dyspersji podłużnej i poprzecznej.

Obliczenia prowadzone wzdłuż linii prądu przy czym założono, tu można na nich wydzielić odcinki o stałej prędkości przepływu. Koncentracja obliczona na końcu da- nego odcinka jest funkcją wejścia dla następnego odcinka o innej prędkości przepływu. Rozpatrywany jest ruch jonów bez opóźnień tzn. nie podlegających adsorpcji czy wymia- nie jonowej. Obliczone zmiany koncentracji dla poszczegól- nych linii prądu są wykorzystane do wykreślenia map zasię- gu, określonego zanieczyszczenia na przedpolu obszaru inie- łrającego.

She paper contains the description of a computer program for calculation of the pollutant movement in groundwaters. The program is based on the known two-dimen- sional solutions of the dispersion equation. Calculations are performed along the streamlines. Each streamline is divided into sections for which constant velocities may be assumed*

The concentration obtained at the end of a given section is an input function for the next section with a different velocity* The undelayed ion movement (i.e without any adsorption or ion exchange) is considered.

Results of calculations serve for drawing^ maps of iso- lines of chosen pollutant concentration at a given time.

Резюме.

В работе представлено описание программы для цифро- вой вычислительной машины, позваляющей вычислять движе- ние загрязнений в грунтовых водах. Программа была соста- влена на основе известных решений уравнения дисперсии в случаях выступания продольной и поперечной дисперсии.

Расчёты проводились по линии тока с условием, что

.«окно на них выделить отрезки с постоянной скоростью

токопрохокдения. Концентрация, вычисленная на конце дан-

ного отрезка является функцией входа для следующего от-

резка с другой скоростью токопроховдения. Рассматривает-

ся движение ионов без отставаний, т. е. не подвергающизся

адсорбции или ионообмену. Расчеты изменений концентраций

для отдельных линий тока используются для вычерчевия карт

дальности определенного загрязнения перед районом иньек-

ции.

Spis treści.

etr.

1. Wstęp 7 2. Sposób prowadzenia obliczeń 8 2.1 Przypadek dyspersji podłużnej 8 2.2 Występowanie dyspersji poprzecznej . . . . . 10 3. Dyskusja założeń 12 3.1 Warunki iniekcji 12 3.2 Uwzględnianie dyspersji poprzecznej 1²» 3.3 Prędkość przepływu ib 4. Program do obliczeń . . . 16 4.1 Opis programu 16 4.2 Wprowadzanie i wyprowadzanie danych 17 4.3 Zastosowanie programu 17 5. Spis literatury <9 6. Wydruk programu * . . . » . ^{2 O}

1. Wstęp.

Przedstawiony w pracy program do obliczeń ruchu za- nieczyszczeń w wodach gruntowych został napisany dla pew- nego uproszczonego modeilu, który powinien, zdaniem auto- ra, w zadowalający sposób odzwierciedlać warunki często spotykane w rzeczywistości* Dla tego modelu przyjęto na- stępujący założenia:

- pole hydrodynamiczne nie ulega zmianom i jest stałe w czasie trwania iniekcji zanieczyszczenia

- następuje całkowite wymieszanie wód infiltrujących z powierzchni z wodami gruntowymi w wyniku dyspersji pionowej i różnic w ciężarze właściwym

- prędkość ruchu zanieczyszczeń jest równa prędkości ruchu wody, gdyż rozpatrywany jest ruch jonów nie pod- legających adsorpcji rozpadowi czy wymianie (np. chlor- ki, siarczany)

- ruch zanieczyszczeń można analizować jako ruch csąstek po torze wyznaczonym w polu hydrodynamicznym przez linie prądu ze składową prostopadłą do nich w płasz- czyźnie poziomej, wynikającą z dyspersji poprzecznej - iniekcja odbywa się w sposób ciągły w czasie az do za-

kończenia w zadanym momencie czasowym z obszaru o pew- nej określonej długości i szerokości (patrz rozdz.

[3.1]).

W wyniku prowadzonych przy użyciu programu obliczeń uzyskiwane są rezultaty zmian koncentracji zanieczyszczeń w czasie w zadanych punktach na podanych liniach prądu lub rozkłady zanieczyszczeń wzdłuż linii prądu w zada- nych odstępach czasowych. W oparciu o te rezultaty można wykreślić mapy zasięgu zadanej koncentracji w podanych odstępach czasowych na przedpolu obszaru iniekującego.

2. Sposób prowadzenia obliczeń»

2.1 Pxzypadek dyspersji podłużnej.

Rówu.enie transportu materiału (zanieczyszczenia) rozpuszczonego v wodzie dla przypadku jednowymiarowego przed3tswia się następująco:

dt

'1вг С = С (x,t) koncentracja zanieczyszczenia w pewnym punkcie x w chwili czasu t,

D - współczynnik dyspersji podłużnej, v - prędkość przepływu wody,

x,b - współrzędne położenia i czasu.

jest wiele rozwiązań tego równania przy różnych, brzegowych i początkowych* ¥ niniejszej pracy zastoscwane rozwiązanie o postaci całkowej:

•g. Л ^ ~ - exc ( -

0

° M

które jest słuszne przy następujących warunkach:

C(x,t)| -s Co dla t > 0 (5.1)

|x=0

C(x,t)J = 0 dla x > 0 (3.2) (tsO

C(x,t)j = O (3.5)

Rozwiązanie to podali między innymi Banks i Ogata [i]

oraz Hubert i wsp. [3]. Dla przypadku zmiennej w czasie iniekcji t.j. przy warunku (?.1) o postaci:

C(z,t)l = Co(t) 2=0

rozwiązanie podał Hubert [2]:

C(x,t) = / C o ( t X - = = = = = - exp{- t x

O ib D C t ) ^ } d

4 D (t-т)

Rozwiązanie (2) i (4) są słuszne dla pewnej stałej pręd- kości przepływu v. Z uw££i jednak na zmienność tej wiel- kości wzdłuż linii prądu rozwiązania te. musiały zostać zmodyfikowane.

Przyjęto, że wzdłuż zadanej linii prądu można wy- różnić odcinki o pewnych długościach x^, J^, ... x^, którym odpowiadają stałe prędkości przepływu u>j, u

, ...

u,. Należy przeto analizować każdy odcinek oddzielnie.

Pierwszy odcinek 2c, ma swój początek w punkcie iniekcji (patrz rozdz. [3.1]) dla całego profilu. 7!m1 any koncen*- tracji w poszczególnych punktach na tym odcinku obliczy- my posługując się wzorem (2) dla u = u^. Funkcję C(x,t) dla x = x» nazwiemy funkcją iniekcji dla odcinka drugie- go. Posługując się teraz wyrażeniem (4-) i przyjmując, że

Co(t) = C ^ i t ) 1 u = u

, oraz x, < x < Xg obliczymy zmiany czasowe koncentracji w punktach na odcinku Xg«

Procedurę powtarzamy dla następnych odcinków. Wyrażając

to w postaci wzorów dla.i-tego odcinka otrzymamy:

t x r [x-ut(t-T)j C(x,t) s / C f a ^ . t V - exp{ -*

0

V F W (t-т)

dla zⁱ_¹< x < x^l (5.1) gdzie

(5.2)

D (t-т)

Przy czym z uwagi na to, że S jest funkcją u wygodnie jest wyrazić je w postaci następującej (np. Zuber

D = - . v = const o v (6)

2*2 Występowanie dyspersji poprzecznej.

W tym przypadlcu równanie transportu przyjmuje postać:

D 3-S - v 3S = .32 (7) 7 a r ак at

gdzie С в C(x,y^tt)

В - współczynnik dyspersji poprzecznej»

у - współrzędna położenia.

10

Pozostałe jak w równaniu (1).

Rozwiązanie tego równania dla przypadku punktowej iniekcji opisanej deltą Diraca podaje Zuber [*] opiera- jąc się na stochastycznych wzorach, wyprowadzonych wcześ- niej przez Scheideggera [б].

Uogólniając je dla przypadku gdy iniekcja odbywa się ze strefy o pewnej szerokości y^Q i długości r^Q = 0,

linią iniekcji, uzyskujemy:

(8.1)

przy czym warunek początkowy można zapisać jako:

C(x,t) dla О ^ у ^ y⁰

C(x,y,t)j = { (8.2) I**0 0 dla у < 0 , у > yQ

gdzie: C(x,t) jest przedstawione w postaci wyrażenia (2) lub (4) w zależności od warunków iniekcji.

Obliczenia dla określonej wartości у (uderzonej po pro- stopadłej do danej linii prądu) prowadzono analogicznie jak to opisano w rozdziale [2.1] przemnażając wartości C(x,t) przez odpowiednie wartości współczynników A(y,v,t) gdyż:

= C(x,t) . A(y,u,t) (8.3)

(9)

11.

Zależność A od u jest wynikiem wprowadzenia w miej=- sce D wzoru (6).

3. D?skus.ia założeń.

3.1 Warunki iniekcji.

Przedstawione w poprzednich rozdziałach wyrażenia przedstawiają rozwiązania równani a dyspersji dla przypad- ków gdy iniekcja jest punktowa (x⁰ = O, y^Q = O) [(3.1) i (3.4)] oraz gdy iniekcja odbywa się na linii s^Q = 0»

у * 0) [(8.2)]. Zazwyczaj jednak iniekcja odbywa się z ob- szaru pewnej długości хр Ф 0 i aby mogły być wykorzystane podane rozwiązania trzeba przyjąć pewien określony model procesu infiltracji zanieczyszczenia z powierzchni do wód gruntowych. V? niniejszej pracy uczyniono to w dwojaki sposób. W pierwszym przypadku zakładamy, że iniekcja od- bywa się v/ strefie o największej przepuszczalności piono- wej na odcinku o tak małej długości w stosunku do odleg- łości, dla której rozpatrujemy zasięg zanieczyszczenia

( x^o« x ) , że można ją uznać za równą zero (z⁰ « 0 ) . Zakładamy więc, że strefę o największej przepuszczalności możD"» uznać za linię iniekcyjną (8.2). W tej sytuacji roz«- wią2._iie tego problemu podaje wyrażenie (6). Drugi model jest bardziej złożony. Zakładamy, że wody skażone przecho- dzą z powierzchni do wód gruntowych na całej długości od- cinka xo. W tej sytuacji z chwilą zakończenia infiltracji z powierzchni pod dnem obszaru na odcinku x^Q będzie pewien przestrzenny układ koncentracji g(§) (-x⁰<£O). Obecność tych wód spowoduje, że mimo zakończenia iniekcji w chwili t = Ц w punkcie x = 0 (zwanym punktem iniekcji) w dal- szym ciągu będą występowały zmiany koncentracji zanieczy- szczenia. W związku z tym należy przeliczyć funkcję iniek- cji dla czasów t > t^. Aby to wykonaćskorzystano z roz-

12

wiązania równania dyspersji dla warunku początkowego opi- sanego:

C(x,t)| = g(£) dla - r

< C < 0 , (10)

Co daje dla t > t>, rozwiązanie:

0 g U ) CC - (x-t

6(x,t) = / j • ' expi

(10.1) punkcie x s 0 rozwiązanie to daje:

c(x,t)i - ?

1

(10.2)

gdzie: v

jest średnią prędkością przepływu wód grunto- wych aa odcinku x

Całkowita funkcja iniekcji będzie wiec sumą dwóch wyra- żeń:

с (t) = C(x,t)| + C(x,t)j (11)

° |x=0 jx=O

Pierwsze z nich odpowiada warunkowi początkowemu opisane-

mu wyrażeniem (5.1) lub (5.4) przy założeniu zakończenia

iniekcji v chwili t = t.,, natomiast drugie odpowiada wyra-

żeniu (10.2).

3.2 Uwzględnienie dyspersji poprzecznej.

Analizę przemieszczania się zanieczyszczenia na boki prostopadle do linii prądu wywołanego obecnością dyspersji poprzecznej przedstawiono w rozdziale (2*2).

Wyprowadzone wzory są słuszne jedynie dla przypadku gdy obszar iniekcji rozciąga się na pewną szerokość y

. Hoz- ważmy jednak przypadki w których wpływ dyspersji poprzecz- nej można zaniedbać. Na rys. 1

oprócz obszaru iniekcyj- nego zaznaczono również poglądowo umownie przyjmowane

skrajnie linie prądu tj. osie у = O oraz у = у • Dodatko- wo na tym rysunku zaznaczono zmiany funkcji A(y,u,t) (9) wzdłuż osi y.

Można nie wprowadzać poprawek (9) na dyspersję po- przeczną przy prowadzeniu obliczań zmian koncentracji za- nieczyszczeń na liniach prądu (profilach) leżących w pew- nej odległości od skrajnych linii. Poprawka A(y,u,t) jest bowiem wtedy równa jedności* Wynika to stąd, że rozkłady poprzeczne zmian koncentracji zanieczyszczeń od poszcze- gólnych profili będą się wzajemnie na siebie nakładały.

Powstałe w ten sposób sumaryczne pole koncentracji zanie- czyszczeń zachowuje się tak jakby zjawisko dyspersji po- przecznej nie występowało. Natomiast dla skrajnych linii prądu oraz całego obszaru leżącego poza nia uwzględnianie dyspersji poprzecznej jest konieczne.

3.3 Prędkość przepływu*

Prędkość przepływu wód gruntowych można określić do-

świadczalnie bezpośrednio przez pomiary wskaźnikowe lub

pośrednio na podstawie prawa Darcy. Mamy wtedy:

С/Со

Kierunek przepływu wód grunt.

Rys. nr1 Poglądowe przedstawienie Unii iniekcji (Xo=O,y=yo) z zaznaczeniem skrajnych linii prqdu (y=O,y=yo).

Dodatkowo pokazano zmiany funkcji A(y,v«t) dla danego punktu X wzdłuż osi y.

gdzie: - - stosunek współczynnika filtracji do poro- watości efektywnej,

j5 _ spadek hydrauliczny.

W tym przypadka dokładność prognozy wynikająca z dokład- ności oszacowania v zależy w dużym stopniu od znajomości wartości k/n.

4. Program do obliczeń.

4.1 Opis programu.

W programie do obliczeń ruchy zanieczyszczeń w wo- dach gruntowych można wyróżnić Trlilm zasadniczych części.

Na początku prowadzone są obliczenia funkcji iniekcji zgodnie z tym co zostało powiedziane w rozdziale [j.t]

w oparciu o zależności (10.2) i (10.5). W programie przy- jęto, że g(C) = C

na całej długości z

. Hie wyklucza to możliwości użycia innego rozkładu g(£). Halety tylko wtedy w SUBROUMNIE K>F0(S) w karcie SO 09 pomnożyć występujące tam wyrażenie przez funkcję g(§). Parametr § odpowiada w programie literze S.

9 przypadku gdy zadajemy z

= 0 następuje skokowe za- kończenie iniekcji ciągłej w punksle x = 0 w podanym mo- mencie czasowym t*•

Do realizacji tych obliczeń używa się podprogramu całkującego metodą kwadratury GUJSSA (10-cio punktową) wywoływanego jako QG10 z funkcją podcałkową nazywaną ТОЮ.

¥ kolejnym etapie obliczane są poprawki na dyspersję po- przeczną A(y,u,t) zgodnie z wyrażeniem (9) przy użyciu SUBBOTTXH QG19 oraz ruru. Wprowadzenie tych poprawek na- stępuje zgodnie z tym co powiedziano w rozdziale (3.2).

Program realizuje przypadek gdy y

-» » , gdy nie trzeba

wprowadzać poprawek (9) poprzez odpowiednie wczytanie

danych.

Przyjęcie yQ, któremu odpowiada w programie zmienna YB równego zero powoduje, że A(y,u,t) = 1 i obliczenia podawane są bez uwzględniania dyspersji poprzecznej, na- tomiast podanie dowolnej ujemnej wartości za yQ oznacza, że szerokość stawu można przyjąć za nieskończoną, Dalej w programie wyróżnić można grupę instrukcji obliczających zagadnienie transportu zanieczyszczenia dla przypadku jednowymiarowego zgodnie z wyrażeniami przedstawionymi w rozdziale [2.1], posługując się subrutynani QG1O, PIFI i FEFE oraz całkując metodą prostokątów (Ralston [5]).

Następnie następuje wprowadzanie poprawek na dyspersję poprzeczną (9) jeżeli uznany to za konieczne.

4.2 Opis wprowadzania i wyprowadzania danych.

Realizowany jest poprzez blok instrukcji wejściowych (READ) oraz wyjściowych (WBITE).

Instrukcje wejściowe są poprzedzielane instrukcjami typu IF, które umożliwiają prowadzenie obliczeń dla wie- lu kompletów danych. Szczegółowy ich opis wraz z podaniem oznaczeń zmiennych zamieszczono w części opisowej wydruku programu (patrz rozdział [б])« Wyprowadzanie danych w po- staci wykresów następuje przy użyciu podprogramu PRTEDX.

4.5 Zastosowanie programu.

Program oblicza rozkłady czasowe zmian koncentracji zanieczyszczenia w ЮТКК punktach każdego profilu (linii prądu). Przy czym KEK oznacza tu ilość odcisków o stałej prędkości na jakie została podzielona dana linia prądu, natomiast L jest to zadana ilość punktów obliczeniowych na każdym odcinku. Następnie w oparciu o uzyskane w tych punktach rozkłady czasowe program zestawia rozkłady prze- strzenne koncentracji wzdłuż całego profilu dla zadanych

HH momentów czasowych (liczonych od początku iniekcji).

Mając rozkłady czasowe w punktach pomiarowych można weryfikować dobór takich danych jak 2 oraz u (lub -), natomiast w oparciu o rozkłady przestrzenne dla wszyst- kich profili analizowanego obszaru można wyznaczyć stre- fy, gdzie koncentracja zanieczyszczeń przekracza wartoś- ci dopuszczalne» Uzyskuje się w ten sposób mapy zasięgu stref skażonych dla analizowanego obszaru dla wybranych momentów czasowych. Program był praktycznie wykorzystany dla obliczenia prognozy ruchu zanieczyszczeń na przedpo- lu stawu osadnikowego w Gilowie. Szczegółowe rezultaty obliczeń podane są w niepublikowanej pracy Uałoszewskie- go i wsp. [?]. Sks6t tej pracy będzie opublikowany w ter- minie późniejszym* Przyjęte obecnie w programie założenia iniekcji ciągłej do chwili t^ oraz ciągłego rozkładu przestrzennego g(g) = C^Q pod dnem stawu nie są wolne od wielu zastrzeżeń jednakże dotychczasowe obserwacje tere- nowe zgodne są z przyjętymi założeni.ami matematycznymi.

Nie ulega jednak wątpliwości, że najistotniejszym czyn- nikiem obok właściwego rozpoznania prędkości przepływu u jest właściwe rozpoznanie i opis funkcji iniekcji.

18

5. Spis literatury.

[i] A. Ogata, R.B. Banks;

A Solution of the Differential Equation of Longi- tudinal Dispersion in Porous Media.

U.S. Geol, Survey Prof. Paper 4-11A, Washington 1967.

[2] J. Hubert»

Effect of injection speed on the tracer concentra- tion curve. Hukleonika—Tom XV-Nr V70 (str 571-376) [5] J. Hubert, A. Lenda, A. Zuber;

A Solution of the dispersion adsorption equation with linear adsorption isotherm.

Nukleonika - Tom XVI - Nr 5-6/1971 (str 271-278) [4] A. Zuberj

Dyspersja wskaźnika przy przepływach przez ośrodki porowate w aspekcie zastosowań hydrogeologicznych.

Zesz. Nauk. AGH, KPCh z. 7^J Kraków 1971.

[5] A. Ralston;

Wstęp do analizy numerycznej, PWN Warszawa 1971.

[б] А.Е. Scheideggerj

App. Physics, tom 25: 1954, s 997.

[7] P. Małoszewski, S. Witcaak, A. Zuber;

Analiza przemieszczania się wód skażonych na przed- polu stanu osadnikowego "Gilów™ oraz ocenrs wpływu wód osadnika na własności geologiczno-inżynierskie gruntów w podłożu.

Cz. IX. Analiza przemieszczania się frontu wód akażonych na przedpolu osadnika "Gilów".

Instytut Hydrogeologii i Geologii Inżynierskiej AGH Kraków, 1975. Praca niepublikowana.

6. Wydruk ргокгатди.

PROGRAM S I L O W d N P O T . O U T P U T , T i P E 6 = I N P U T . T A P F 8 = O I J T P U T l С . . . . C»c»

c<c»

c»c»

C«

PROGRAM GILOM ŚLUZY DO OBLICZAMI «AR7OSCI KONCENTRACJI WZGLĘDNEJ C/CO WSKAŹNIKA INIEKOVANE6O 00 WODY GPIJNTOWEJ « SPOSÓB C U G L Y 00 CHWILI Т = П * P-CIE X J « A S I E ł.lpf) CZASIE T] WYSTĘPUJE NAOAL IMEKCJ».KTOR* WYNIKA i ŁINtOWEGO ROZKŁADU KONCENTRACJI POD DNEM STAkU)

С Tl = CZAS 00 ZAKOŃCZENIA INIEKCJI * MIES.

С TJ« CZ«S J*KI UPŁYNIE 00 OSTATNIEJ PROGNOZY • MIFS.

С ОТ* ODSTĘP CZASOWY (ОТ=Тг/500 MAX>

С A»D/V ( CONST. . METRY ) С AP«OY/V (CONST..METRYI

С в» STOSUNEK PRZEPUSZCZALNOŚCI 00 POROWATOŚCI w VI OOHE

С L » ILOSC « E S C I NA JAKIE JEST PODZIELONY K A Z O Y ODCINEK PROFILU «' L«S I С КККг ILOSC ODCINKÓW NA PROFILU ( 00 10)

С NN'ILOSC «ARTOSCI CZASO T WCZYTYWANEGO w TF«JI.

С N= ILOSC CZESCI NA JAKIE JEST PODZIELONY CALV OKBES CZASU T2.

С XX(К) DŁUGOŚCI POSZCZEGÓLNYCH ODCINKÓW ,« DANYM PB0F1LU С W ( К ) ODPOWIADAJĄCE IM PRĘDKOŚCI FILTRACJI W M/KIES.

С X0= OLUGOSC STAWU POZA TAMA W METRACH

С YS SZEROKOŚĆ STAWU w METRACH /PODANIE VB=o OZNACZA NIE UWZGLEONIA>

С NIE DYSPERSJI POPRZECZNEJ.A YB<0 OOPOWTAOA NIESKOŃCZONOŚCI/

С YE 0DLE6LOSC -OSI Y OD POCZĄTKU UKL*OU WSP, W METRACH С V0> SREONIA P4EDK0SC FILTRACJI NA ODCINKU >0 W И/^ICS.

С TFtJ) KACIERZ WARTOŚCI CZASU T DLA KTÓRYCH CHCEPY "IEC PODANY С ROZKŁAD KONCENTRACJI w PRZESTRZENI.

(••••••••••••••••••••••••••••••••••••••••••••••••••«••••••••••«•«••«•••••••••••«

С Y M J ) MACIERZ WSPÓŁCZYNNIKÓW DO UWZGLEOMAKl* DYSPERSJI POPHZECZ.

С ZX<M) HACIERZ WARTOŚCI ODLEGŁOŚCI К N* KA7DYM ODCINKU.

С YY(N.H) MACIERZ WARTOS. С/СО Ы CZASIE T ORAZ 4 P-CIE X (T«DT«N, С Ж>0Я»М,Л=^Х1К-])«0К»М» BEZ UWZGLĘDNIANIA DYSPERSJI POPRZECZ.

С rVIJ.N) MACIERZ YY Z UWZGLEONIENIEM OYSPERSJI POPRZECZNEJ С PTT IN) MACIERZ WARTOŚCI CZASU T.

С XT(I) MACIERZ WARTOŚCI ODLEGŁOŚCI WZOLUZ C*LESO PROFILU.

С VTUiJl MACIERZ WARTOŚCI KONCENTRACJI W Р-кТАСН XT DLA USTALONEGO С CZASU T PODANEGO W TFIJ).

EXTERNAL FUNC.FEFE.FIFI."OFO.FUFU 001 COMMON/AAB/ AP.VP.T

COMMON /BAB/ A.X,V.JJ.L«XV.V0.TO _ 00?

COMMON /ABA/ YY(SOO.S( 003 COMMON /6LPRT/ H»L?.L3.L*.LS 004 DIMENSION VEC(SO0).TF(5I.PTT(SO0),YV(5O0.5I.XX(10).VV(10!.ZX(SI 095 DIMENSION YN(S00.SI,«T(5CI). W Y N O 0 0 ) .YTtSO.SI 006 DIMENSION YP(a00WTT<Z0G>iYR(8l .YA(SOO) 007 2 F0RMAT<5Xf36HK0NCENTRACJA SKAZEN PRZY ТАИJf(X»0.)//) 006 t, FORMAT<3X,2H7«.EJ1.*.2X.SMC/CO«,SI2X.E11.«.>> 009

S FORMAT(21X.?HX»,SI2X,£ll.*t/) 010

7 F0HMAT(3IS.7E4.2) 01]

9 F O R H » T < 8 X . 2 H T J . E 1 1 . 4 , S X . 5 H C / C 0 « . £ 1 1 . < . I 012 12 F0RHAT(10X.3KYE«.F6.1/l

13 F0RMAT<IHl.I7X.S(3X.aHT».EI0.3.3X)>

14 FORMATИОХ.гнА». El1.4/10Х,3HT1».E11.4/) 014 16 F0RMATC6E11.4) 01S

20

IT FORMAT<SX,2MX«.E10.3.S<2X.5HC/CO«.U1..4I> 014 2? F0R>«AT<lH1.21X.2HV=.EU.4/> 017 S2 FORM»T<5X«4HX«0.>«X.5(2X.5HC/|.;=.E1I.4>> 01R 90 FORMAT 12E11.4) ' 014 49 FORH»T110X.4HK/N=iF6.1.6t-M/008E/> 020 77 RE*0<6»7>KK«.L« NN.».OT.T].T2.X0.V0iR 021 IF(KKK.EO.O) со то зов огг REA0(6«90><(XXSK>.VV(K>)tK*l»KKK> 023 RE»0(6.16> CTFCJI .J*lt№»> 024 19 RE«0(6.16>AP.YB.YE

IFMP.EO.O.) GO TO 18

KRitEie.n.) A . T I 025

«RITE 18.94) в 0?*

МЯ1ТЕ1в<12) YE LI'

L3«L*«

L5«

0?7огв огчоэо

031 N»IFIX<T?.'OT) 043

тэ»тг-п «3603%

NK'tFIX (ТЭ/ОТ) 017 м»ГГ1Х(Т1/0Т) ОЗв YF«O.

00 35 J 4 . N VP'VT»DT»J

31 IF1YE.E0.0.IGO TO 341 C»Lt OGlOtO.tYE.FUFU.YF)

во то 3S

зг Y*CJ)>I.O GO TO 35

ЭЭ C»LL OG10(-YE»Ye-YE.FUFU.YF)

IS CONTINUE i 00 IS.9 J0-1.M 039]

VEC(JO>«OT»JO 0401 1S9 VEC(N*JO>«1.0 061]

XY«-XO 043 WRITE! 8.2! 046, Y0«0. { 00 ISO IJ«J.NK П45 IFIXO.EO.O.IGO TO IS*

TO»OT»1I 047' CALL Q6l0(XY.0..F0F0.Y0) 048;

|F(YO.LE.O,001> YO.0.0 049 158 YP(III*Y0 050 -

VEC(M>II)'TTtIII 051 VEC(HM»It)«Yf><U> 052 .140 CONTINUE . _04 3.

00 51 J«1»NN 054 Z»TF«J> O S S 1V.JR.N 057 41 YR(J)»VEC<IV) • 04»

ПО 1*0 !I«1<NK ПЯЧ 140 *RlTEie.9)TT(lIl.rP(II) 060

21

C»LL PSTi5l.T(0>VfC<M.?.O.F4NC.T?>o..l.»0.> 0*1 DO 20 1=1.L • »6F K»nx«FLO»T(l) 0*3 Z X I I ) = * **»

X T U I s * 065 Jil Oftfc I T«OT»J 0 6 7 PU<J)=T 0<>8 I F I T . O T . T l ) GO TO (ПО ОЬ9 V»T 070 C*LL QG10C0.»U.FIFl,Y«l 071 I F i Y X . L E . 0 . 0 0 1 1 YM=0.0 073 I F I Y i l . S E . 0 . 9 7 1 YW=1.0 0 7 3 Y Y I J . l ) = r n 074 J»J«1 07S GO TO 1 0 7 6 10» CONTINUE 077 YX«0. • 0T«

00 160 l Q i l t N K 0 7 9

и=т-(Т|»(г.о*ю«1.о)/г.о»от) овс

IFIU.LE-O.I GO TO 1*0 0»l y*»YPUO>»FIFHU>«OT 08?

YX«YX»Y« ОвЭ 1*0 CONTINUE 00*

1F(Y<.LT.0.001) YX=0.0 O8S YYU.I)*YX Овб

«•T-Tl 047 UU«T 0Я8 C»LL 0G101O.UU.F1F1.YO) Of4 IF(YO.LE.O.OOI1YO*0.0 . 090 IF1Y0.HE.0.971 Y0»1.0 091 Y7(J.l>*YY(J.t)«Y0 092 IF(T.GE.TZ) 00 TO ?0 093 JsJ.l 09*

GO TO 1 095 JO CONTINUE 096

00 71 I=1.L 00 12 J»liN

12 Yy(J.Il=YYIJ.I)»Y»UI 71 CONTINUE

чр.1Те«8.гг> v ' 097

WRITEle.SI IZJC(I).I»I.LI 098 DO ?S Jl=l.N 099

?S 4RITE(R.f.)PTT(Jl). (TV(J1 .11.1x1.1.) 100 00 «.5 U l . L lfll 00 53 J4.HN Юг

?«TF(J( 103 JK»IFIXU/OT) 104 S3 YTI|.J)«YVIIK.I1 105 4S CONTINUE ' 106 00 14 U l . L 107 00 55 J I ' l . N 108 J1»J1 109 VEC<JI1»PTT(J1) HO KKaN«|.JI , , , SS VEC(KK)cVVIJlil) | | ; IS CONTINUE 113 LL»L«l • Ц 4 C*UI. Pf»TPLllK.*EC.N.LL.0.FUNC.y2.0.. 1..0.) 115

XN.XX1U U6

KL>KKK>1 It7

22

ЪО 300 LKs|«KL MB 1 1 * OX*X»(Kl/FL0ATIL> i го

YF»O.

00 86 J*1«N VP»VT=OT»J 1F(YB)81.82.83 8) IF(YE.EO.O.)GO TO 8*

CALL 06iniO..YE.FiJFO.rr) a* YA<J>»O.S»YF

GO TO 86 вг YMJI-I.O

GO TO 86

83 CALL OGIO(-YE.YB-YE.FUFU.YF) Y»(J)«YF

вб COMTINUE oo io I»I«L

«TIKDaXN 1?' DO Й0 JS*1.N 1?6 80 YV(J5.I)»0. I?⁹ 00 30 JU«ltN 4°

DO 40 J2-1.N 131 T»PTT(J2) 13?

IF|J?.LE.JU) SO TO 3 133 KU*a»ju-i 13»

J J . J I J 13S, U"T-0.5"0T««O 13b YZ»FEFEIU>«OT 137 GO TO 11 138 Э Y2«0. 139 11 CONTINUE 1*0 IF(YZ.oe.0.971 YZ»!.O 1*1

IFIYZ.LE.o.oon тг«о.о |*г YN(J2«I)»YZ 1*3 fcO CONTINUE 1**

00 50 J3«1>N 1*5 50 YV(J3.I)»YV(Jl.It»rN<J3«l) 1*6 30 CONTINUE 1*9 00 «0 J**).« 150 1^

60 YV(J4,II»YY(J<..D»r*(J4( 15?

00 65 IN»1»NN IS»

I f I 156 YT(Kl.iN)>YVIIK.I) 1ST 65 CONTINUE I5T 10 CONTINUE ISe WRITE(В.гг) V 1S9 НЯ1Ге<В>5> (ZX(l).l-).L) 160 00 AS Jl*l»N 161 85 WRITE<8<4| PTT(JL)>(YV(JL<|I«I>1»L) 16?

00 70 l«l.L 163 00 75 JL»1«N 16*

JK«JL 165 VEC(JK)«PTTIJL) 166 KP«N«I»JK 167 75 VEC(KPt«YVIJL<I> 1*и

7* CONTINUE 169

23

CALL PRTPLHK«VEC»l?«U-"«>»FVINC.T3.0..l..0.l 170 200 CONTINUE 171 LI«KKK«L 17г WBITE1в.131 I T F ( J ) > J ' l > h K I 173 'Л1ТЕ1в«521 (T41J) .J»1.NNI 17*

DO 110 K 2 « i » U 175 WVNM>>0. 176 KR«KZ»1 !77 WYMKRXXTCKZ) 178 MRITECetlTi XT(KZI>(TT(K2rI>tl'ltNN) 179 110 CONTINUE 180 00 120 J«1.NN 181 KB««LI»lł»J»l 182 oo i3o кг>1>(.: ie«

Kt-xe-кг IBs 130 wrN(KA>«rTCKZ>j) lee 12S COIKTIHUe 187 L J « U M 188 KtiKN*] 189 C«LL Pi>TPLT(9<«YMtLJaK^>O.FUNCiXN«0.«l.^>0.) 190 00 TO 19

IB CONTINUE

GO TO 77 191 300 CONT1NUC 192 END 193

FUNCTION F l f l l T T I FI01 COMMON /e«8/ *»«tV.JJ.t_,XT,VO.TO FI02 OIMENSICK ]XZI6> FI03 ОН» 1ХЛ/6Ч-0)/ FIO»

|XXI*I«O FI05 CM.L SYSTEMC(115«I«X) FI06 G«V»(X/*"TT»»ł«/W-TTI/(4.»»«TTI FIOT GG«*.»3.1*l*»»»<TT"3ł«V FI08 FIFl«X»E»Pt-G)/SORT(061 FI09 RETURN FU0 ENO FtU

FUNCTION FEFEITTI FEOI

o v» ^EO2

COfWON /»0»/ TTCS00»5l ГЕ03 DIMENSION IXX(61 " FE04 0»T» lXK/6*(-0l/ FE05 1ХХ|*|«в FE06 C«tL S»STEHCIH5.IX*» FE07 R«V»«X/¥-TT»»łX/V-TTł/«t.»»»TTI FE08 RR«*.»3.1*16»*»ITT«»3)«V FE09 FeFE*rYUJ,D«X«£*?l-«>/SO«T(RRI FE10

«TORN , FEU END . ' FE12

FUNCTION rUPUtS) COMM0N/AAB/ AP.VP.T DIMENSION IXXIfc) DATA ке/6'i-ot/

IXX(*>sO

CALI SYSTEMCM15<IXX>

W«5»5/<4.»AP»VP»T) .••<..»3.l416»AP»VP*T FUFU«£XP <-•»/SORT(»«) RETURN

END

FU01 FU03 FUO*

FUOS FU06 FU07 FUOS FU09 FU10

run

FIKCTIOW

COFIMON lute/ «,X,V,JJ.L.XY,VJ,T.

01 1£NSI3H 1ХХЁ1 0AT*

CALL SYSTćKCIllb';IXX)

RETURN

FO 1Г0.2 F 0 . 3 FOLI.

fO. 5 fOi.f>

FO. ? fO.9FOll fOłi

SUBROUTINE O010(XL.XU»FCt.Yi

Y».03333S67»<FCT(»»C)«FCT<»-C)ł C«.*32S317»B

Y»Y».07*72567*(FCT<»»C)»FCT(»-CI) C«.33<J70«i8«8

Y«Y».109S432«(FCT(A*CI*FCTIA-C)>

C».2166«77*B

Y"Y•.13*633*»1FCTIA»C)«FCT«A-CI1 C«.07*«3717»B

Y»B«(Y».1*77621»(FCT(A«CI»FCT(»-CI»»

RETURN ENO

SUBROUTINE FUNC(X.Y) Y " 0 .

RETURN EM)

OGIOOG10

esioOGIO

воюOOIO 0610 06100610 OGIO OGIO QG1CGIO OGIO0610

го10₃₀

*oso fie70

e90 100 110120 l*v13»

ISO

FU010 FU020 FU030

SUBROUTINE SCAL. IXSCAL.XMAX.KMJN) DUMMY SUBS TO «VOID ROUNOEO SCALES

RETURN END

SCO 10 SC020 SC030 SC040

SUBROUTINE P R T P O •(NÓ.A.N.M,NFUNC<FUI«C.«I AX.XLIN.YLAX.YLI») PLOTO0O0 PURPOSE - TO PLOT * GRAPH nllH ONE 1Ч0ЕРЕк0ЕЫГ VABIABLE «NO UP PLOTOOIO TO 20 OEPENOENT VARIABLES. «ITH THE ADDITIONAL AR1LI1Y TO PLOT * PLOT0020 C»LCUL»TEO CUHVE. THE: IftOEPEMDENJ VAHIMJLF IS PLOTTED ON » PLOT0030 HORIZONTAL AXIS. THE DEPENDENT ONES ON A VF»TICAL AXIS. WIDTH »LOT00*O IS IflO PRINT POSITIONS» HEIGHT IS SO. EVERT POINT OF EACH PLOT0050 DEPENDENT VARIABLE IS IhOICATEO BY ;> NU»BFP 11-01. WHILE THE PLOT0060 CALCULATED POINTS ARF OENOTEO BY ASTERISKS. PLOTOOTO PARAMETER USAGE. PLOTOOSO ER PLOT0090 NO A F1XEO POINT NUMBER. UP TO 1 DIGITS- PRINTED AS THE

A A VECTOR «HOSE FIRST N POSITIONS CONTAtN THE INDEPENDENT PLOT0I10 VARIABLE. ANO WHOSE NEXT M SETS OF 4 POSITIONS CONTAIN PLOT0120 THE DEPENDENT VARIABLES PLOT0130 N NUMBER OF OBSERVATIONS PLOTOlfcO M NUMRER OF VARIABLES IINOEPENOENT . OEPFNOEWT) PLOTS1S0 NFUNC GREATER THAN ZERO IF A CALCULATED CURVE IS TO BE PLOT0160 PKtNTED PLOT01T0 FUNC SUBROUTINE TO GENERATE CALCULATEO CURVE. IF ONE WAMTEO. PLOTO 180 ELSE IS A OUHHY. PROGRAM CALLING PLOT MUST HAVE AN PLOT0I90

•EXTERNAL FUNC*. SUBROUTINE CALLED BY CALL FUNC (X.Y). PLOT0200 WHERE X IS GIVEN TO SUBROUTINE »Nf> Y RFTURNEO. PLOT0210 XLAX.XLIN.YLAX.YLlN MAXIMUM ANO MINIMIIM VALUES OF THE PL0T02?0 INDEPENDENT ANO DEPENOENT VARIABLES TO ЯЕ USED IN THE PLOT0230 PLOT IF XLAX • XLIK» THE PROGRAM CALCULATES ITS OWN PLOT0240 MAXIMUM AND MINIMUM FO» THE IN0EPFNOENT VARIABLE. PLOT0S50 SIMILARLY FOW YLAJ< • YLIM PLOT0260 REQUIRED SUBROUTINES. FUNC (IF USED). AND SCAL. PLOTO?70 01 HENS I ON I0UT(10?)tXPR(in>A(l).CALCll02).IREP<?ll

OIMENSION PIPIID.PAPAtll COMMON /HLPRT/ Ll.L?.L3.L<>.LS

EQUIVALENCE 1IOUT!1>.XPRCI>) PLOT0330 DATA PIPI/10H /

DATA PAPA/IH*/

DATA IREP /lH2.|Hl.)H?.lH3.1Ha.lHSilH6.1H7tlMe.lH9<IH0>

1IHA.1HB.1HC.1HD.JHE.1HF.1HG.1HH.IHI.1HJ/

С CALС IS 1 LARGER THAN NEEDED IN CASE XMIN . IOO*XSCAL PLOT0290 С SHOULD 8E LARGE» THAN XMAX. THIS PREVENT* SLOPOVER INTO PLOT0300 С NEXT LOCATION. LOOK AROUND CARD '950 TO SEE «MAT I MEAN. PLOT0310 С CALC IS WHERE CALCULATED FUNCTION GOES. PVOT0320 1 FORMATOH1.60X.7H CMAHT .13) " PLOT0340 2 FORMAT U « .1PE12.3.?H. .1O1A1.1H*,1PE12.3)

7 FORMAT IIH .K.X.10I10H*. .).1HM PLOT0370 9 F0RMATIJH0.9X.1P11E10.3)

118 FORMAT U H .12X.2H- . 1 0 1 А Ы Н - 1 С PRINT CHART NO.

«RITE (8.1) NO ' ICOUNT < *

С IF NO EXTREMES OF * GIVEN, FIND THEM IF (XLAX - XLIN) 20.10.20

10 XMIN ' A l l ) XMAX * XMIN 00 15 J • l.N

IF 1A(J) . XMINI U t l Z d f 12 IF (MJI - XMAJt) 15.15.14 11 XMIN > » O ]

GO TO 15 1» XMAX « A(J>

~Г5~С0МТШЯ GO то гог

20 XMAX»XLAX XM1N«XLIN

С CALCULATE ЯЫ SCALE SIZE гог XSC»L«(XMAX-XHIN)/IOO.

С ROUTINE TO CALCULATE EXACT SCALE SUE AND END POINTS

PLOT03S0 PLOTO390 PLOT0*00 PLOTO*10

рсооо*го

Pi.OT01.30 PLOT04*0 PLOTO«S0 PLOT0*60 PLOT0*70 PLOT0480 PLOT0«90 PLOT0500 PLOTOS10 PLOT 0S<!0 PL0T0S30 PLOTOSAO PLOT0550 PLOT0560 PLOT0570 .PLOT0580

26

CALL SC»L {XSCAL.XM*X.XM]N> PLO?0590 С IF NO EXTREMES OF Y GIVEN. FIND THE» PLOT0600 IF (YLAX - YLIN) 110.112.110 PLOT0610 U ? L « N • 1 PLOT0620 YM1N > »(L) PLOTS630 YKAX'YMIN PLOT0640 LL« M»h PLOT0650 00 «0 J > l_.Lt PLOT0660 IF(A(J)-YMIN> 26.26.26 PLOT0670 26 IFIA(J)-YMAX) 40.40.30 PLOT0680 28 YHINsA(J) PLOT0690 GO TO 40 PLOTOTOO 30 YMAX»A(J) PL0T0710 40 CONTINUE PLOT0720 GO TO 201 PLOT0730 110 YMAX.YLAX PLOT0740 YHIN»YLIN PLOT07SO С GET SCALE SIZE ANO END POINTS PLOT0760 201 YSCAL • 1YMAX - YHIN ) / 50. PLOT0770 CALL SCAL (YSCAL.YHAX.TMIN) PLOT0780 С PRINT TOP SCALE PLOT0790 XPBfl) • XHIN PLOT0800 DO 200 JP'1.10 PLOTOeiO XPRIJP'l) » XP4IJPI • XSCAL • 10. PLOT0820 С TO MAKE SURE THAT ZERO REALLY PRINTS AS ZFPO. NOT * SMALI NUhSER PLOT0830 С CAUSED BY ROUNGING ERRORS. PLOT0840 IF <ABSIXPR<JP.|) - .5 • XSCALI) 240.240.700 PLOT08S0 2*0 XPRIJPM) • 0. PLOT0860 200 CONTINUE PL0T0870 WRITE <8.9) IXPH(JP).JP«1,1D PLOTC880 WRITE (8.7) PLOT0890 С IF CALCULATED CURVE WANTED GET VALUES BETWEEN W\H AND ХИАХ PLOT0900 IF (NFUNO 310.210.211 ' PLOT0910 211 F • XMIN PLOT0920 JP « I PLUT0930 213 CALL FUNC(F.CALC(JP>) PLOT0940 fF IF - «MAXI 212.210.210 PLOT09SO 212 F « F • XSCAL PLOT0960 JP • .JP « 1 PLOTOvfO GO TO 213 PLOT0980 210 CONTINUE PLOT0990 С START PRINT AT MAXIMUM Y PLOT1000 tPR • YfAX PLOT1010 С CLEAR PRINT LINE PLOT1020 230 00 SS JP » I.101 PLOT 1030 55 lOUT(JP) >PIPI(1) PL031040 С IF CALCULATED CURVE WANTED SET UP POINTS PLOT10S0 IF (NFUNO 214.214.215 PLOT1060 215 F « XMIN PLOT)070 С SCAN ALL VALUES .OF Y FOR X 6ETWEEN XHIN «UJJ ХИАХ PLOT 1080 JP » 1 PLOT1090 С IS POINT WITHIN HALF A SCALE OF PRINT-POSITION PLOTU00 220 IF (ABS(YPR-CALCIJP)) - .5 • YSCALI 216.217.?18 ' PL0T1I10 С IF EXACTLY BETWEEN PRINT POSITIONS ONLY PPIfJT IT ONCE PLOT 1120 21T IF <YPR - CALC U P ) ) 218.Z16.216 PLOT!130 С RELIEVE IT OR NOT THIS IS AN ASTERISK (NUMgEO TOO LARGE TO WRITE PLOT11*0 С AS ONE NUMBER) PLOT]ISO С EXTRA 6AIO LINES AT • AND - TEN PER CENT

218 IF (ABSIYPR-.1)-.S»YSCAL) 216,245.244 2*4 IF IA8S<YPH-.1>-.S»YSCALI 216,243.246

С IF EXACTLY BETWEEN PRINT POSITIONS PRINT IT 04LY ONCE 245 IF(YPP.-.1)246.216.216

243 IF(YPR*.1)316.216.246

С EXTRA GRID LINES AT REGION BOUNDARIES

27

?** IF <JP-L1> 2S0.216.2*7 2*7 IF IJP-L2)2S<>.216,2*8 г*8 IF IJP-L3I2SO,216.2*9 8*9 IF <JP-L4I25O.216.2495 2*«5 I F I J P - L S > 2 S 0 . 2 1 6 . 2 S O 216 lOUTIJP» • PAPA(l)

2S0 IF (F-ХИДХ) 219.21«.21«

г» JP » JP »i

F • F • XSCAL GO TO 220

С RUN DOWN СДСН SET OF DEPENDENT VARIABLES С IF NO POINTS WANTED

21» IF (N) 70.70.300 300 DO 221 J * 2»H

00 22? L • l.N

С CALCULATE SUBSCRIPT FOR » LL • «J - 1) • N » L

С IS IT WITHIN HALF A SCALE Of PRINT POSITION IF <A8S<YPR - AILLM - -S • VSCAL) 223,224.225 С IF EXACTLY HALFvAY BETWEEN, PRINT ONLY OHCF

22* IF (YPR - A<LL1> 225.223.223 С FINO HORIZONTAL POSITION

123 JP • (A(L) -XMIN) / XSCAL • 1.5 С IF OFF GRAPH, FORGET IT

IF U P - u 225,226,226

226 IF (JP • 1011 227.227.225

С THIS GIVES 1.2.Э ETC. FOR J-2.3.4 ETC 227 IOUT U P ) • I R E P U )

225 CONTINUE ггг CONTINUE 221 CONTINUE

70 ICOUNT « ICOUNT «1

С PRINT VALUE ON VERTICAL AXIS EVERY FIVE POSITIONS IF I ICOUNT - S> 120.119.120

120 WRITE (в.Ив) (I0UTIJP)«JP«l<101) GO TO 80

С *ДК£ ZERO PRINT AS ZERO, NOT SHALL NUMBER 119 IF (ABSIYPR I-.S • VSCAL1 232,232(233 232 F« 0.

GO TO 23*

233 F • YPR

234 WRITE <в,2) F, <I0U1<JP>.JP«l,10!l»F _ ICOUNT * 0

С IF REACKEO ^V'*IN« STOP вО IF IYPR - Yl-IN) в6,86.<.5 С ELSE OECREMENT Y

<i5 YPR » YPR - YSCAL GO TO 230 86 WRITE 18,7) С PRINT BOTTOM SCALE

XPR11) • XMIN 00 90 JP «1.10

XPBIJP.1>«XPR(JP)*XSC»L«1O.

IF (ABSUPRIJP'l) - .S*XSCALH 231.231.90 231 «PR U P » 1) « 0.

90 CONTINUE

WRITE 18,9) <XPRIJP>»JP«)»1I) WRITEIS,999>

999 FORMATdNO) RETURN END

PLOT H I * PLOT 11«»

PLOT 1201 P L O T 1 2 U PLOT12?ł PLOT123»

PLOT 12*»

PLOT12S»

PL0T126t PL0T127»

PLOT 1280 PLOT 1290 PLOT 1ЗОв PLOT1310 PLOT i зга PLOT1338 PL0T1349 PLOT 1350 PLOT1360 PLOT1370 PLOT1380 PLOT1390 PLOT1*00 PLOT1*20 PLOT1*30 PLOTI*«O PLOT1*50 PLOT 1460 PLOT1*70 PLOT1480 PLOT! (.90 PLOT1S00 PLOT1S10 PLOTI520 PLOT1S30 PLOT1S*0 PLOT1S50 PLOT1S60 PL0015TO

PLOTI590 PLOT 1600 PLOTI610 PLOTI620 7LOTI630 PLOT1640 PLOT16S0 PLOT1460

•Ч.ОТ1670

PLOT1680 PLOT1690

28