diff --git a/README.md b/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..bcb98f430b98026fd4f1b207347cf4f5752d027e
--- /dev/null
+++ b/README.md
@@ -0,0 +1,79 @@
+# supersvd
+Python tool to perform maximal covariance analysis and calculate empirical orthogonal functions
+
+Написанная на языке `Python 3` программа состоит из двух файлов: `main.py` и `supersvd.py`. 
+В файле `supersvd.py` находится алгоритм вычисления максимально скоррелированных мод, а в `main.py` — вспомогательный код, 
+который анализирует ключи запуска программы, делает чтение входных данных из файлов, а также записывает в выходные файлы результаты работы алгоритма.
+
+Функцию `supersvd` можно напрямую использовать из кода на `Python`, в этом случае не обязательно сохранять массивы в виде файлов на диске.
+
+## Описание
+
+Функция `supersvd` по двум заданным наборам пространственно–временных полей строит матрицу ковариации, а затем вычисляет её неполное сингулярное разложение. 
+
+Функция `supersvd` принимает на вход 2 обязательных аргумента (два поля, максимально скоррелированные моды которых мы ищем) и 2 опциональных параметра: 
+количество пар максимально скоррелированных мод (по умолчанию 3) и значение переключателя режима вычитания из поля его среднего по времени значения 
+(по умолчанию `True`, то есть из поля _вычитается_ его среднее по времени значение).
+
+Пусть $`X(t), Y(t)`$  — два меняющихся во времени поля, максимально скоррелированные моды которых мы ищем, причём $`\operatorname{dim}(X)=nT \times nX`$[^1] и $`\operatorname{dim}(Y)=nT \times nY`$, где $`nX`$ и $`nY`$ могут быть одним или несколькими измерениями массивов (в случае среднемесячных данных INMCM $`nX`$ и $`nY`$ — $`120 \times 180`$). 
+
+[^1]: Здесь и далее размерности массивов указаны в порядке, принятом в `C` и `Python`. В `Fortran` размерности массивов следует развернуть в обратном порядке.
+
+Функция `supersvd` вычисляет разложение вида:
+```math
+\begin{aligned}
+X(t) &= \overline{X} + XV_1 XC_1(t) + XV_2 XC_2(t) + \ldots + XV_k XC_k(t) + \ldots,\\
+Y(t) &= \overline{Y} + YV_1 YC_1(t) + YV_2 YC_2(t) + \ldots + YV_k YC_k(t) + \ldots,
+\end{aligned}
+```
+где 
+```math
+\overline{X}=\frac{1}{nT}\sum_{t=1}^{nT} X(t),\qquad
+\overline{Y}=\frac{1}{nT}\sum_{t=1}^{nT} Y(t),
+```
+а $`k`$ — количество пар максимально скоррелированных мод. В формуле выше каждое новое слагаемое получается максимизацией корреляции между $`XC_k(t)`$ и $`YC_k(t)`$, а
+$`XV_k, YV_k`$ — два семейства ортогональных пространственных мод.
+
+Моды $`XV_k, YV_k`$ являются левыми и правыми сингулярными векторами матрицы ковариации 
+```math
+C =\frac{1}{nT} \sum_{t=1}^{nT} (X(t) - \overline{X}) (Y(t) - \overline{Y})^{\mathsf T}.
+```
+Функция `supersvd` возвращает:
+  * массивы `x_coeff`, `y_coeff`  временных коэффициентов $`XC(t), YC(t)`$ разложения (размерности $`k \times nT`$);
+  * массив `x_vect`  левых сингулярных векторов $`XV`$ (размерности $`k \times nX`$);
+  * массив `y_vect` правых сингулярных векторов $`YV`$ (размерности $`k \times nY`$);
+  * массив `corrcoeff`, содержащий $`k`$ коэффициентов корреляции между $`XC_k(t)`$ и  $`YC_k(t)`$;
+  * массив `x_variance_fraction` (`y_variance_fraction`), содержащий доли дисперсии, приходящиеся на каждый из $`k`$ левых (правых) сингулярных векторов;
+  * массив `eigenvalue_fraction`, содержащий долю дисперсии матрицы ковариации, приходящуюся на $`k`$-ую пару сингулярных векторов;
+  * массив `eigenvalues` сингулярных значений матрицы ковариации $`C`$.  
+
+
+## Использование
+
+Функция `supersvd` может вызываться как и из другой `Python`-функции, принимая на вход массивы данных, так и из командной строки, принимая на вход бинарные файлы (`.STD`). Последняя возможность реализована в функции `main`.
+
+ Функция `main` принимает на вход 3 обязательных аргумента: 
+  * `-x` имя файла, содержащего первое из полей (например, `X.STD`);
+  * `-y` имя файла, содержащего второе из полей  (например, `Y.STD`)[^2];
+  * `-t`, `--time` длину временного интервала (например, в случае среднемесячных данных исторического эксперимента с INMCM это 165 лет).
+
+ [^2]: Если нужно посчитать ЭОФы, то в качестве первого и второго нужно задать одно и то же поле, то есть передать два раза имя одного файла.
+ 
+ Также функция `main` принимает 7 необязательных (опциональных) параметров:
+  * `--type` тип используемых данных — `real` (4 байта) или `double` (8 байт), значение по умолчанию — `real`;
+  * `-k` количество вычисляемых пар максимально скоррелированных мод, значение по умолчанию — 3;
+  * `-xv` имя файла, в который запишется массив `x_vect`;
+  * `-yv` имя файла, в который запишется массив `y_vect`;
+  * `-xc` имя файла, в который запишется массив `x_coeff`;
+  * `-yc` имя файла, в который запишется массив `y_coeff`;
+  * `-stat` имя файла (предпочтительно в формате `.CSV`), в который для каждого $`k`$ запишутся домноженные на 100% элементы массивов: `corrcoeff`, `x_variance_fraction`, `y_variance_fraction`, `eigenvalue_fraction`.  
+
+Функция `main` также может быть запущена с ключом `--dont-subtract-mean`: при этом из полей $`X`$, $`Y`$ _не будут вычитаться_ их средние по времени значения.
+
+Итак, чтобы вычислить с помощью функции `main` 4 максимально скоррелированные моды аномалий температуры и давления (типа `float`) и сохранить все возможные результаты, достаточно в командной строке выполнить:
+```bash
+python3 main.py -x ts.std -y ps.std -t 1147 -k 4 -xv tsv.std -yv psv.std \
+                -xc tsc.std -yc psc.std -stat c.csv
+```
+
+Информацию, сохраняемую в файл `c.csv`, функция `main` также выводит на экран.