автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Полулагранжева модель динамики атмосферы мезомасштабного разрешения с использованием метода конечных объемов

На правах рукописи

Шашкин Владимир Валерьевич

Полулагранжева модель динамики атмосферы мезомасштабного разрешения с использованием метода конечных объемов

Специальность 05.13.18 — математическое моделирование, численные методы и комплексы программ

2 8 НОЯ 2013

Автореферат

диссертации на соискание учёной степени кандидата физико-математических наук

Москва-2013

005540060

Работа выполнена в федеральном государственном бюджетном учреждении науки Институте вычислительной математики Российской академии наук

Научный руководитель: доктор физико-математических наук,

Толстых Михаил Андреевич

Официальные оппоненты: Головизнин Василий Михайлович,

доктор физико-математических наук, профессор, федеральное государственное образовательное учреждение высшего профессионального образования «Московский государственный университет им. М.В. Ломоносова», факультет вычислительной математики и кибернетики, кафедра вычислительных методов,

Холодов Ярослав Александрович, кандидат физико-математических наук, федеральное государственное образовательное учреждение высшего профессионального образования «Московский физико-технический институт (государственный университет)», доцент.

Ведущая организация: федеральное государственное бюджетное учреждение науки Институт вычислительной математики и математической геофизики СО РАН

Защита состоится «20» декабря 2013 г. в 1400 часов на заседании диссертационного совета Д 002.045.01 при федеральном государственном бюджетном учреждении науки Институте вычислительной математики Российской академии наук по адресу: 119333, г. Москва, ул. Губкина д. 8

С диссертацией можно ознакомиться в библиотеке федерального государственного бюджетного учреждения науки Инсппута вычислительной математики Российской академии наук.

Автореферат разослан |3 2013 года.

Ученый секретарь

диссертационного совета Д 002.045.01

доктор физико-математических наук

Общая характеристика работы

Актуальность темы. Долгосрочный прогноз погоды и моделирование изменений климата невозможны без вычислительных экспериментов с математическими моделями климатической системы Земли. Повышение точности прогноза и уменьшение неопределенности моделирования климата - главная цель развития подобных моделей. Естественный путь достижения этой цели - совершенствование отдельных компонентов моделей климатической системы, в том числе моделей общей циркуляции атмосферы. Большое значение для развитии моделей атмосферы имеет увеличение разрешения.

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

На данный момент в России для целей моделирования климата применяется модель общей циркуляции атмосферы Института вычислительной математики РАН (ИВМ РАН) с горизонтальным разрешением около 150-200км (Володин, 2010)1, а для целей долгосрочного прогноза - модель ПЛАВ (ПолуЛагранжева, основанная на уравнении Абсолютного Вихря)2 с разрешением 100-150 км и другие модели более грубого разрешения. Данные модели способны разрешать на сетке явления синоптического масштаба. В недалеком будущем возможности вычислительной техники позволят применять для долгосрочного прогноза погоды и моделирования климата модели общей циркуляции атмосферы, способные разрешать на сетке а-мезомасштабные процессы (т.е. модели с горизонтальным разрешением 30-50 км). В мире ведутся работы по созданию подобных моделей.

В большинстве моделей, применяемых для долгосрочного прогноза погоды и моделирования климата, используются эйлеровы численные методы. Шаг по времени в подобных методах ограничен условием устойчивости Куранта (или более жестким условием в случае методов высокого порядка аппроксимации). При горизонтальном разрешении 30-50 км допустимый шаг по времени, определяемый условием Куранта, оказывается неоправданно малым с точки зрения точности.

Шаг по времени не ограничивается условием Куранта при использовании полунеявного полулагранжева подхода к интегрированию уравнений динамики атмо-

1 Володин Е.М., Дианский H.A., Гусев A.B. Воспроизведение современного климата с помощью совместной модели общей циркуляции атмосферы и океана INMCM 4.0 // Изв. РАН. Физика атмосферы и океана. 2010. Т. 46. № 4. С. 448-466.

2Толстых М.А. Глобальная полулагранжева модель численного прогноза погоды. М.; Обнинск: ОАО ФОП, 2010. 111с.

3

сферы (Robert, 1985)3. На практике, полунеявный полулагранжев подход позволяет проводить расчеты с шагом по времени в 3-5 раз большим, чем при использовании эйлеровых методов, и, таким образом, экономить вычислительные ресурсы. Дополнительным преимуществом полунагранжева подхода является высокая вычислительная эффективность при расчете переноса большого количества составляющих атмосферы (малых газовых составляющих и аэрозолей). Данное свойство имеет большое значение, так как современные модели долгосрочного прогноза погоды и моделирования климата включают описание динамики порядка 100 таких составляющих.

При моделировании атмосферы с высоким пространственным разрешением преимущества полунеявного полулагранжева подхода в совокупности с возможностью эффективной реализации на вычислительных системах с 0(Ю4) процессорами (White, 2011)4, (Mtieler, 2013)5 делают его многообещающей альтернативой явным эйлеровым методам. Возможность реализации полулагранжевых алгоритмов на сетках с квазиравномерным разрешением на сфере («кубическая сфера», редуцированная ширагно-долгогаая сетка) позволяет уменьшать шаги сетки по горизонтали до 0(10) км.

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

Проблему сохранения массы можно решить, применив локально-консервативный конечно-объемный вариант полулагранжева подхода (Machenhauer, 2009)6. При этом важно применять подобные методы как для расчета переноса компонентов атмосферы, так и для дискретизации уравнения неразрывности (уравнения сохранения массы). В случае применения разных методов для дискретизации уравнения неразрывности и уравнений переноса компонентов атмосферы может нарушаться глобальное и локальное сохранение массы переносимых величин, что было показано в работе (Jockel, 2010)7 в случае эйлеровых конечно-объемных методов.

3 Robert A., Yee Т., Ritchie Н. A semi-Lagrangian and semi-implicit numerical integration scheme for multilevel atmospheric models // Mon. Wea. Rev. 1985. Vol. 113. P. 388 - 394.

4J.B. White Ш, Dongarra J.J. High-performance high-resolution tracer transport on a sphere // J. Comput Phys. 2011. V. 230. P. 6778 - 6799.

5Muller E., Scheichl R. Massively parallel solvers for elliptic PDEs in numerical weather- and climate prediction, submitted preprint 2013. p. 24. URL: http://arxiv.org/abs/1307.2036vl.

6Machenhauer В., Kaas E., Lauritzen P. H. Finite-Volume Methods in Meteorology // Computational Methods for the Atmosphere and the Oceans. ELSEVIER, 2009. Vol. 14.

7Jockel P., von Kuhlmann R., Lawrence M. G. et al. On a fundamental problem in implementing flux-form advection schemes for tracer transport in 3-dimensional genera] circulation and chemistry transport models // Q. J. Roy. Meteorol. Soc. 2001. Vol. 127. P. 1035-1052.

Целью данной работы является создание модели динамики атмосферы (динамического блока модели атмосферы), способной разрешать на сетке а-мезомасштабные процессы и предназначенной для долгосрочного прогноза погоды и моделирования климата.

Для достижения поставленной цели была создана модификация динамического блока модели ПЛАВ, сочетающая известные преимущества этой модели (точность, устойчивость и высокая вычислительная эффективность при мезо-масштабном разрешении) с необходимыми для долгосрочного моделирования свойствами локального и глобального сохранения массы. Создание модификации динамического блока модели ПЛАВ, обладающей необходимыми свойствами, подразумевало реализацию конечно-объемного полулагранжева метода для дискретизации уравнений неразрывности и переноса и включало следующие этапы:

1. Создание и тестирование конечно-обьемного полулагранжева алгоритма численного решения уравнения переноса на сфере.

2. Разработка численной модели мелкой воды на сфере с конечно-объемной полу-лагранжевой дискретизацией уравнения неразрывности для исследования свойств подобной дискретизации в упрощенном, двумерном случае.

3. Реализация в модели ПЛАВ конечно-объемной полулагранжевой дискретизации уравнений неразрывности и переноса.

4. Проверка точности и устойчивости созданной модификации модели ПЛАВ в стандартных экспериментах по воспроизведению динамики атмосферы и сравнение с результатами зарубежных моделей.

Научная новизна:

1. Впервые в России реализован локально-консервативный полулагранжев алгоритм численного решения уравнения переноса на сфере.

2. Впервые в мире подобный алгоритм реализован на редуцированной широтно-долготной сетке.

3. Независимо от группы исследователей Метеослужбы Великобритании (Wood, 2013)8, впервые в мире реализован блок решения уравнений динамики атмосферы для модели общей циркуляции атмосферы, основанный на конечно-объемном полулагранжевом подходе в сочетании с полунеявной схемой интегрирования по времени, сохраняющий массу локально и глобально.

Практическая значимость. Динамический блок, представленный в настоящей диссертации, может использоваться в моделях общей циркуляции атмосферы для долгосрочного прогноза и моделирования изменений климата, в частности, в глобальной полулагранжевой модели атмосферы ПЛАВ. При этом ожидается: • увеличение точности прогноза за счет повышения пространственного разрешения;

8N. Wood, A. Staniforth, A. White et al. An inherently mass-conserving semi-implicit semi-Lagrangian discretisation of the deep atmosphere global nonhydrostatic equations // Q. J. Roy. Meteorol. Soc. 2013.

• улучшение воспроизведения распределений компонентов атмосферы за счет использования локально-консервативного алгоритма расчета переноса.

Апробация работы. Результаты, вошедшие в диссертацию, докладывались на: 52-ой научной конференции МФТИ (2009), конференции по методам решения дифференциальных уравнений на сфере (Partial Differential Equations on the Sphere, Потсдам, Германия, 2010) и конференции по прогнозу погоды и климата на суперкомпьютерах нового поколения (Weather and Climate prediction on Next Generation Supercomputers: Numerical and Computational Aspects, Эксетер, Великобритания 2012), международной школе-конференции по вычислительно-информационным технологиям для окружающей среды CITES-2011 (Томск, 2011) и CITES-2013 (Петрозаводск 2013), на семинаре ФГБУ Гидрометцентр России (2012).

Полностью диссертация докладывалась на семинаре «Математическое моделирование геофизических процессов» в ФГБУН ИВМ РАН (2013).

Диссертационная работа была выполнена при подцеряже Российского фонда фундаментальных исследований (гранты 10-05-01066, 12-05-31441 и 13-0500868), федеральной целевой программы «Кадры» Минобрнауки РФ (контракт №14.132.21.1378, соглашения №№ 8344, 8350, 8326).

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

Публикации. Основные результаты по теме диссертации изложены в 6 печатных изданиях, 4 из которых изданы в журналах, в том числе 3 изданы в журналах, рекомендованных ВАК, 2 - в тезисах докладов.

Объем и структура работы. Диссертация состоит из введения, трех глав, заключения и приложения. Полный объем диссертации 108 страниц текста с 30 рисунками и 12 таблицами. Список литературы содержит 114 наименований.

Содержание работы

Во введении обосновывается актуальность исследований, проводимых в рамках данной диссертационной работы, приводится обзор научной литературы по изучаемой проблеме, формулируется цель, ставятся задачи работы, сформулированы научная новизна и практическая значимость представляемой работы.

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

Полулагранжев подход к решению уравнения переноса основан на его лагранжевой форме:

¿И' V

утверждающей, что концентрация переносимой величины с/ в лагранжевой частице (частице, перемещающейся вместе с воздухом) не меняется при = 0 (т.е. при отсутствии источников/стоков переносимой величины).

Полулагранжев подход к вычислению концентрации д^1 на шаге по времени п+1 в некоторой точке вычислительной сетки Хф заключается в рассмотрении движения лагранжевой частицы Х{1) такой, что Х(1п+1) = Хф. Кинематическое уравнение для данной частицы

где V - вектор скорости ветра, интегрируется на шаг по времени назад с помощью итерационного алгоритма. Таким образом, вычисляются координаты точки Х(Ьп), в которой находилась выбранная частица в момент времени 1 = 1'\ Точка Х*]1: = Х(Г) называется исходной точкой траекторий с конечной точкой в узле сетки Хф.

В силу уравнения переноса концентрация переносимой величины в лагранжевой частице сохраняется, следовательно, г/™].1 = где - значение концентрации в исходной точке Хф на шаге по времени п. Как правило, исходная точка траектории не совпадает ни с одним из узлов вычислительной сетки, поэтому для нахождения значения концентрации применяется интерполяция.

Недостатком традиционного полулагранжева подхода является отсутствие сохранения массы переносимой величины. Данная проблема может быть решена при использовании конечно-объемной формулировки полулагранжева метода, основанной на интегральной форме уравнения переноса:

4 / № = 0, (3)

dtJsv{t)

где 5У{Ь) - произвольный лагранжев объем (объем, перемещающийся вместе с воздухом), / = - плотность переносимой величины, р - плотность воздуха.

7

Рис. 1: Слева - схематичное изображение исходного объема и ячейки вычислительной сетки. Справа - аппроксимация исходного объема многогранником.

Физический смысл уравнения переноса в интегральной форме - сохранение массы переносимой величины в лагранжевом объеме.

В конечно-объемных полулагранжевых методах рассматривается лагранжев объем 5У{£) такой, что 6У(Ьп+1) совпадает с некоторой ячейкой вычислительной сетки АУук- Плотность переносимой величины Д^1, осредненная по этой ячейке на шаге по времени п+1, задается выражением

где АУ1*к = 5У(Г1) - исходный объем, т.е. объем, совпадающий с рассматриваемым лагранжевым объемом на шаге по времени п. Физический смысл интеграла по исходному объему в правой части уравнения (4) - масса переносимой величины, заключенная в исходном объеме.

Вершины углов исходного объема Д У£к - это исходные точки траекторий лагранжевых частиц с конечными точками в вершинах углов ячейки вычислительной сетки Л Уф- Координаты вершин углов исходного объема вычисляются так же, как координаты исходных точек в традиционных полулагранжевых методах. Исходный объем ДУ?к полагается равным многограннику, полученному соединением вершин его углов отрезками (см. Рис. 1 лев.).

В разделе 1.1 рассматривается применение конечно-объемного полулагран-жева метода к решению уравнения переноса в одномерном случае. Для нахождения плотности переносимой величины Д1+1, осредненной по некоторой ячейке вычислительной сетки Д"Ц = \х1 — ц + |Дх], необходимо интегрирование /™ по исходному объему /\У* (который полагается известным). Функция плотности /п задана набором средних значений /к по ячейкам сетки. Для вычисления /г™+1 с точностью выше первого порядка плотность в ячейках вычислительной сетки аппроксимируется кусочно-параболической функцией Ьк(х)=ак+Ькх+Скх2, х&Ук,

(4)

k=l,...,N (Colella, 1984)9, коэффициенты которой определяются выражениями:

[ hk{x)dx = &x~fl (5)

Jvk

h(xk±^Ax)=fn(xk±^Ax), (6)

где значения плотности на границах ячейки /"(:г;,:±|Д.т) вычисляются с помощью интерполяции. Условие (5) означает равенство средних значений функции плотности / и ее приближения h в ячейке Vk и необходимо для сохранения массы.

При численном решении уравнения переноса важно сохранение монотонности решения, т.е. не должны появляться новые экстремумы концентрации переносимой величины, а существующие максимумы не должны увеличиваться (минимумы - уменьшаться). Использование кусочно-параболической функции распределения плотности h не гарантирует сохранение монотонности, поэтому к функции h применяется монотонный фильтр (Barth, 1989)10.

В случае решения двумерного или трехмерного уравнения переноса (разделы 1.2 и 1.3) производить интегрирование по точной геометрии исходного объема сложно и вычислительно неэффективно. В работе (Nair, 2002)11 на примере двумерного уравнения переноса на сфере было показано, что аппроксимация исходных объемов многоугольниками со сторонами, параллельными координатным осям, существенно упрощает вычисления и лишь незначительно влияет на точность. В трехмерном случае логично аппроксимировать исходные объемы многогранниками с гранями, параллельными координатным плоскостям (см. Рис. 1, прав.).

Выбор аппроксимирующего многогранника (многоугольника) можно осуществить различными способами. В работе (Nair, 2002)12 для двумерного случая был предложен такой выбор аппроксимирующего многоугольника, который позволяет вычислить массу, заключенную в исходном объеме, посредством двух последовательных одномерных интегрирований вместо одного двумерного (т.н. каскадный подход). Каскадный подход позволяет в многомерном случае не строить многомерные функции, аппроксимирующие распределение плотности, а использовать одномерные и, таким образом, уменьшить вычислительное время. Кроме того, к одномерным функциям плотности легче применять монотонные фильтры.

В разделе 1.2 диссертации описывается алгоритм численного решения двумерного уравнения переноса на сфере. Автор диссертации обобщает двумерный

'Colella P., Woodward P. The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations // J. Comput. Phys. 1984. Vol. 54. P. 174 - 201.

10Barth T. J., Jespersen D. C. The Design and Application of Upwind Schemes on Unstructured Meshes // 27th Aerospace Sciences Meeting. 1989.

"Nair R. D., Machenhauer В. The mass-conservative cell-integrated semi-Lagrangian advection scheme on the sphere // Mon. Wea. Rev. 2002. Vol. 130. P. 649 —667.

12Nair R., Scroggs J., Semazzi F. Efficient Conservative Global Transport Schemes for Climate and Atmospheric Chemistry Models // Mon. Wea. Rev. 2002. Vol. 130. P. 2059 - 2073.

каскадный подход (Nair, 2002) на случай редуцированной широтно-долготной сетки, в которой разрешение по долготе увеличивается при приближении к полюсам. В разделе 1.3 описывается разработанная автором диссертации аппроксимация исходных объемов в трехмерном пространстве, допускающая каскадный подход к вычислению масс (т.е. использование трех одномерных интегрирований вместо одного трехмерного). Данная аппроксимация является основой конечно-объемного полупагранжева алгоритма численного решения трехмерного уравнения переноса на сфере на редуцированной широтно-долготной сетке, описанного в разделе 1.3.

Созданный алгоритм численного решения уравнения переноса испытывался на задаче «твердое вращение с колебаниями по вертикали» (Jablonowski, 2008)13, которая является стандартным тестом для алгоритмов численного решения уравнения переноса в трехмерном случае. Поле горизонтальной скорости описывает вращение атмосферы, как целого, вокруг географической оси (второй вариант - вокруг оси, наклоненной под углом а = тг/4 к географической). Движения по вертикали - колебательные, с амплитудой 1000 м. Начальное распределение - непрерывное, локальное с горизонтальным радиусом около 2000 км, вертикальным - 1000 м или 2000 м. Точное решение после одного оборота совпадает с начальным распределением.

Ошибки численного решения созданного алгоритма в этом эксперименте оказались в 3 раза ниже, чем в случае неконсервативного полупагранжева алгоритма, основанного на трикубической интерполяции. Кроме того, созданный алгоритм сохраняет массу локально и глобально. При использовании монотонного фильтра (Barth, 1989) 1\- и ^-нормы ошибок численного решения разработанного алгоритма увеличивались до 2,5 раз, однако, оставались меньше соответствующих норм ошибок неконсервативного полупагранжева алгоритма. При численном решении на редуцированной сетке с количеством точек на 20% меньше, чем в регулярной сетке с таким же разрешением на экваторе, нормы ошибок численного решения увеличивались незначительно.

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

I3Jablonowski С., Lauritzen P. H., Taylor М. et al. Idealized test cases for the dynamical cores of Atmospheric General Circulation Models // A proposal for the NCAR ASP 2008 summer colloquium. 2008. URL: http://vnvw.cgd.ucar.edu/cms/peI/asp200S/idealized_testoases.pdf.

Автор диссертации разработал и реализовал в полунеявной полулагран-жевой модели мелкой воды (Tolstykh, 2002)14 конечно-объемную полулагранжеву дискретизацию уравнения неразрывности, что позволило совместить устойчивость при больших шагах по времени с локальным и глобальным сохранением массы.

В разделе 2.1 приведены уравнения мелкой воды на вращающейся сфере:

= _уФ-УФ5 (7)

§ = -ФЯ. ^ (8)

Используются следующие обозначения: V = (u,v) - вектор скорости ветра, Г2 - вектор угловой скорости вращения Земли, г - радиус-вектор точки, Ф - геопотенциал - высота слоя воздуха, умноженная на ускорение свободного падения д, Ф8 - высота поверхности земли, умноженная на д, (...)я - обозначает взятие горизонтальной проекции векторной величины в скобках (с учетом кривизны земли), D = V ■ V -дивергенция поля ветра.

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

Конечно-объемная дискретизация уравнения неразрывности основана на его интегральной форме:

4 [ ®dS = 0, (9)

где SV(t) - произвольный плоский лагранжев объем. Интегральная форма уравнения неразрывности линеаризуется относительно постоянного фонового значения геопотенциала <T>"'f\

Ф'dS=-f ФrefDdS, (10)

dtJsv(t) J6V{t)

Ф' - отклонение геопотенциала от фонового значения. Дискретизация этого уравнения по полунеявной схеме приводит к выражению:

Ф'"+15,(ДУ)+^АгФге/£'п+15,(ДУ)= [ [ф'п-^Д4Фге/Бп1сг5, (11) 2 Jav L 2 J

где AV - некоторая ячейка вычислительной сетки, S(AV) - ее площадь, AV* -исходный объем, соответствующий AV. Интегралы по исходным объемам AV* вычисляются с помощью конечно-объемного полулагранжева алгоритма численного решения уравнения переноса из раздела 1.2.

l4Tolstykh М. A. Vorticity-Divergence Semi-Lagrangian Shallow-Water Model of the Sphere Based on Compact Finite Differences //J. Comput. Phys. 2002. Vol. 179. P. 180-200.

Система дискретных уравнений модели разрешается относительно неизвестной переменной Ф'п+1. В результате получается уравнение Гельмгольца, которое решается с помощью метода конечных объемов.

В разделе 2.3 описываются результаты численных экспериментов с моделью мелкой воды на сфере (Tolstykh, 2002) с разработанной автором диссертации конечно-объемной дискретизацией уравнения неразрывности (11). Были проведены эксперименты «квази-геострофический балланс», «обтекание изолированной горы потоком», «квазиреальные данные» из набора (Williamson, 1992)15, а также эксперимент «баротропная неустойчивость» (Galewski, 2002)16. Данный набор экспериментов является мировым стандартом для тестирования моделей мелкой воды на сфере.

Численные эксперименты показали, что конечно-объемная полупагранжева дискретизация уравнения неразрывности не оказывает отрицательного влияния на точность численного решения уравнений мелкой воды. Несомненным преимуществом конечно-объемной дискретизации по сравнению со стандартной полулагран-жевой дискретизацией, использовавшейся в модели (Tolstykh, 2002) ранее, является локальное и глобальное сохранение массы. Сравнение с результатами, доступными в литературе, говорит о том, что точность разработанной модели при эквивалентном пространственном разрешении не уступает точности зарубежных моделей.

При численном решении на редуцированной сетке с количеством точек на 20% меньше, чем в регулярной сетке с таким же разрешением на экваторе, нормы ошибок численного решения увеличивались незначительно. Например, на Рис. 2 показано, что в тесте «Квазиреальные данные» структура и амплитуда поля ошибки численного решения на регулярной сетке и редуцированной сетке практически не отличаются.

В третьей главе приводится описание разработанной автором диссертации полулагранжевой модели динамики атмосферы (динамического блока модели атмосферы) с использованием метода конечных объемов. Разработанная модель является модификацией модели ПЛАВ (Толстых, 2010) с использованием конечно-о&ьемной дискретизации уравнений неразрывности и переноса влажности. Применение конечно-объемных дискретизаций позволило достичь локального и глобального сохранения массы в полулагранжевой модели.

В разделе 3.1 приводятся уравнения модели - примитивные уравнения динамики атмосферы в а = р/р,-системе координат по вертикали, где р - давление,

15Williamson D., Drake J., Hack J. et al. A standard test set for numerical approximations to the shallow water equations in spherical geometry II J. Comput. Phys. 1992. Vol. 102. P. 211 — 224.

"Galewsky J., Scott R., Polvani L. An initial value problem for testing numerical models of the global shallow water equations // Tellus A. 2004. Vol. 56. P. 429 - 440.

-150 -100 -50

50 100 150 300

Рис. 2: Задача «квазиреальные данные», регион вокруг северного полюса (до 70° с.ш.). Контуры - поле высоты слоя воздуха (м) после 5 дней моделирования, заливка - ошибка поля высоты (м). Слева - численное решение на регулярной широтно-долготной сетке, справа - на редуцированной широтно-долготной сетке (на 20% меньше узлов).

a ps - приземное давление:

(f+2Q х f) н = (12)

(13)

J^ = -RdTv, (15)

di = Fq- 06)

В уравнениях (12-16): Ф - геопотенциал ст-поверхности (высота, умноженная на ускорение свободного падения д), V - оператор горизонтального градиента (сг = const), R,i - газовая постоянная сухого воздуха, Tv = T(l + (S—l)q) - виртуальная температура, Т - температура, q - удельная влажность, <5 = Rq/Rd - отношение газовых постоянных влажного и сухого воздуха, CIxi - теплоемкость сухого воздуха при постоянном давлении, Т - постоянное фоновое значение температуры, а -вертикальная скорость в сг-системе координат, остальные обозначения аналогичны уравнениям мелкой воды. Члены Fy, Ft, Fq уравнений (12-16) отвечают за вклад процессов подсеточного масштаба (испарение и конденсация влаги, конвекция, пограничный слой и др.).

Для формулировки конечно-объемной полулагранжевой дискретизации уравнения неразрывности (14) это уравнение переписывается в интегральной форме:

^ [ р.с№ = 0, (17)

atJsv(t)

где 5У(£) - произвольный лагранжев объем.

В разделе 3.2 описывается полунеявная полулагранжева дискретизация уравнений (12-15) в стандартной версии модели ПЛАВ. Для дискретизации по времени применяется двухслойная полунеявная полулагранжева схема с устойчивой экстраполяцией по времени ЗЕТТЬБ (НоЛа1, 2002)17. Для дискретизации уравнений по пространству используется метод конечных разностей. Система дискретных уравнений модели сводится к решению двумерного уравнения Гельмгольца относительно неизвестной дивергенции Ип+1 на каждом уровне по вертикали.

В разделе 3.3 описывается разработанная автором диссертации модель динамики атмосферы (версия динамического блока модели ПЛАВ с конечно-объемной дискретизацией уравнений неразрывности и переноса влаги). Уравнение неразрывности в интегральной форме (17) линеаризуется относительно фонового значения приземного давления рте^

[ \у.(ргс}У)+Ргг}~]йУ, (18)

cttJ5V(t) Jsv(t)1 ост \

где р'$ - отклонение приземного давления от фонового значения.

Уравнение (18) дискретизуется по времени по схеме БЕЛЬВ следующим образом:

е (г Вп-п1 Г г Я/т^+^ет \

+ 2 ( [^■{Рге/Уп)+рге1—\ АУ+¡^ ]¿V 1, (19)

где ДУ - некоторая ячейка вычислительной сетки, а АУ* - соответствующий ей исходный

объем, (...)(п+1)« = 2(...)" - (...)"_1.

Для вычисления интегралов по исходным объемам используется конечно-объемный полулагранжев алгоритм численного решения уравнения переноса, представленный в главе 1.

Уравнение переноса влаги дискретизуется в интегральной форме так же, как в алгоритме численного решения уравнения переноса (глава 1). Источники влаги (Рд) учитываются явно, результирующее дискретное уравнение выглядит

17Hoital M. The development and testing of a new two-time-level semi-Lagrangian scheme (SETTLS) in the ECMWF forecast model // Q. J. Roy. Meteor. Soc. 2002. Vol. 128. P. 1671-1688.

следующим образом:

-fn+1AV= f (fn + AtpsFn)dVJ = psq (20)

JAV-

Раздел 3.4 посвящен вопросам реализации разработанной модели динамики атмосферы на параллельных вычислительных системах. Для параллельной реализации применяется метод одномерной декомпозиции расчетной области по широте, т.е. каждый процессор производит вычисления в определенной полосе широт. В настоящий момент модель реализована на вычислительных системах с общей памятью и использует 16 процессоров (максимальное число процессоров на общей памяти) с эффективностью 70%.

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

В разделе 3.5 описываются результаты численных экспериментов с разработанной моделью динамики атмосферы. Были проведены эксперименты по моделированию адиабатической атмосферы: «Гидростатическое обтекание изолированной горы потоком» и «Бароклинная неустойчивость», а также неадиабатические эксперименты: «Моделирование циркуляции атмосферы на длительный промежуток времени при наличии внешнего воздействия» и «Моделирование атмосферной циркуляции на 72 часа» (численный прогноз погоды). Эксперименты «Гидростатическое обтекание изолированной горы потоком», «Бароклинная неустойчивость» и «Моделирование циркуляции атмосферы на длительный промежуток времени» входят в стандартизованный набор экспериментов международного проекта по сравнению динамических блоков моделей общей циркуляции атмосферы DCMIP (http://www.earthsystemcog.oig/projects/dcmip/).

Начальные условия эксперимента «Бароклинная неустойчивость» (Jablonowski, 2006)18 представляют собой атмосферу в геострофическом равновесии с небольшим возмущением зональной скорости ветра. Начальное состояние воспроизводит такие особенности реальной циркуляции, как струйные течения в средних широтах и тропопаузу. Возмущение зональной скорости нарушает геострофический баланс в районе максимума скорости струйного течения в средних широтах Северного полушария и вызывает развитие бароклинной волны.

'"Jablonowski С., Williamson D. L. A baroclinic instability test case for atmospheric model dynamical cores // Q. J. Roy. Meteorol. Soc. 2006. Vol. 132. P. 2943 - 2975.

Эксперименты с созданной моделью проводились на регулярных широтно-долшгных сетках с разрешением 100 км, 63 км, 50 км и 33 км на экваторе. При самом высоком разрешении модель была устойчива при тех же числах Куранта, что и при самом грубом, шумы и другие негативные эффекты увеличения разрешения не возникли. Была обнаружена сходимость численных решений на более грубых сетках к решению на сетке с самьм высоким разрешением.

Развитие бароклинной волны по результатам эксперимента с разработанной моделью согласуется с результатами зарубежных моделей общей циркуляции атмосферы, опубликованными в литературе и на странице проекта DCMIP в сети интернет. Нормы отклонения численного решения разработанной модели от численного решения спектральной полулагранжевой модели высокого разрешения, общепринятого в качестве эталонного результата, меньше уровня неопределенности численного решения данной задачи, вычисленного в (Jablonowski, 2006). В экспериментах с разработанной моделью струйное течение в Южном полушарии остается невозмущенным, в отличие от ряда моделей (ICON-DWD, OLAM), использующих гексагональные вычислительные сетки (сетки из треугольных или шестиугольных элементов). В подобных моделях структура сетки вызывает ложные нарушения геострофического равновесия и развитие волны №5 (т.н. «отпечаток сетки»).

Полная энергия атмосферы в экспериментах по моделированию бароклинной неустойчивости уменьшается, т.е. разработанная модель диссипативна. Отток энергии около 0,5 Вт/м2, обнаруженный в эксперименте на сетке с разрешением 100 км на экваторе, оказался в 2,5 раза больше, чем в случае модели (Gassmann, 2013)19 с эквивалентным горизонтальным разрешением, динамический блок которой сохраняет полную энергию с машинной точностью (отток энергии в данной модели возникает за счет введения явной диссипации, необходимой для устойчивости). Однако, величина 0, 5 Вт/м2 значительно меньше потери энергии, наблюдаемой в большинстве современных моделей климатической системы Земли, вследствие неконсервативности и несогласованности некоторых параметризаций подсеточных процессов (Lucarini, 2011)20.

Цель эксперимента по моделированию циркуляции атмосферы на длительный промежуток времени (Held, 1994)21 - проверка точности воспроизведения характеристик атмосферной циркуляции, осредненных за длительный промежуток времени. В этом эксперименте изотермическая атмосфера, находящаяся в состоянии покоя с небольшими возмущениями, приводится в движение с помощью внешнего

I9Gassmann A. A global hexagonal C-grid non-hydrostatic dynamical core (ICON-1AP) designed for energetic consistency // Quart. J. Roy. Meteorol. Soc. 2012. Vol. 139. P. 152-175.

20Lucarini V., F. Ragone. Enei^etics of climate models: net energy balance and meridional enthalpy transport // Reviews of Geophysics. 2011. Vol. 49. 29 p.

2'Held I. M., Suarez M. J. A proposal for the intercomparison of the dynamical cores of atmospheric general circulation models // Bull. Am. Meteorol. Soc. 1994. Vol. 75. P. 1825 — 1830.

воздействия, создающего разницу температур полюс - экватор. Производится интегрирование уравнений динамики атмосферы на 1200 дней. Характеристики возникшей циркуляции осредняются за последние 1000 дней интегрирования (первые 200 дней отводятся на развитие и установление движения атмосферы).

На Рис. 3 показаны поля зонального ветра и температуры, осредненные по долготе и времени (за последние 1000 дней интегрирования). Положение и величина максимумов скорости ветра с хорошей точностью совпадают с результатами других моделей, доступными в литературе, в том числе с результатами моделей Немецкой метеослужбы и Европейского центра среднесрочных прогнозов погоды (Jablonowski, 1998)22. Распределение температуры на Рис. 3 хорошо согласуется с моделями ICON-LAP (Gassmann, 2013) и спектральной модели (Held, 1994), сохраняющими полную энергию (без учета диффузии).

Рис. 3: Эксперимент «Моделирование циркуляции атмосферы на длительный промежуток времени». Слева - зонально осредненное поле температуры, справа - зонально осредненное поле компоненты и скорости ветра, прерывистые линии соответствуют отрицательным значениям (восточный ветер). Поля осреднены с 201-ого по 1200-й день интегрирования.

В данном эксперименте большое значение имеет воспроизведение вихревых потоков тепла и импульса, так как бароклинный вихревой перенос является главным механизмом перераспределения этих величин от полюса к экватору. Интенсивность вихревого переноса можно оценить по полю вихревой кинетической энергии. Характеристики вихревого переноса в эксперименте с разработанной моделью (потоки тепла и импульса, вихревая кинетическая энергия), изображенные на Рис. 4, хорошо согласуются с результатами зарубежных моделей, в том числе с результатами моделей, сохраняющих полную энергию (Gassmann, 2013), (Held, 1994).

22Jablonowski С. Test of the Dynamics of two global Weather Prediction Models of the German Weather Service: The Held-Suarez Test: Master's thesis: Metorological Institute of the University of Bonn, Germany. 1998. URL: http://www-personal.umich.edu/~cjablono/comparison.html.

К'=(и'г+у'2)/2, (тг/зг)

ц | 1 • 1 . . 1 , , 1. , • и!/ 118

-Л А

■ И1 1\Н 1 '

У/ \ч 1:

: 1 --1-4- К( м

уТ, (К»т/я)

90Э бОЭ ЗОв

30N 60И ЭОИ и'У, (тг/аг)

1 1 г~

905 60Э ЗОБ О ЗОЫ 60« 90К

Рис. 4: Эксперимент «Долгопериодное интегрирование». (1) - Зонально осредненная вихревая кинетическая энергия К'. (2) - Зонально осредненный вихревой поток тепла г/Т". (3) - Зонально осредненный вихревой поток зонального импульса и'ь'. Все поля осреднены по времени с 201-ого по 1200-й день интегрирования. Прерывистые линии соответствуют отрицательным значениям.

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

Точность прогноза погоды по большинству показателей соответствовала стандартной версии модели ПЛАВ, применяемой в качестве основной модели глобального прогноза погоды в Гидрометцентре РФ. Отмечено некоторое усиление ложного орографического резонанса (паразитный эффект в численном решении, возникающий во всех полулагранжевых моделях) по сравнению со стандартной

версией модели. Разработанная модель, в отличие от стандартной версии модели ПЛАВ, сохраняет массу атмосферы и водяного пара с точностью машинной арифметики. Эффекты несохранения массы не проявляются при моделировании на короткие промежутки времени (такие, как 72 часа), но в более длительных экспериментах (сезон, год) неизбежно приведут к ошибкам моделирования.

В заключении приведены основные результаты работы, которые заключаются в следующем:

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

2. Разработан динамический блок модели общей циркуляции атмосферы (модель динамики атмосферы) мезомасштабного разрешения с использованием конечно-объемного полулагранжева подхода. Динамический блок может использоваться в моделях общей циркуляции атмосферы для целей долгосрочного прогноза и моделирования изменений климата, в частности, в глобальной полулагранжевой модели атмосферы ПЛАВ. Разработанный динамический блок воспроизводит атмосферную циркуляцию с точностью, соответствующей современному мировому уровню, и устойчив при высоком разрешении.

3. Показана принципиальная возможность применения численных методов, разработанных в диссертации, для моделирования атмосферы с разрешением около 10 км (за счет возможности работы на редуцированной широтно-долготной сетке).

Публикации автора по теме диссертации

1. Шашкин В.В. Локально-консервативный полулагранжев алгоритм решения уравнения переноса на сфере // Труды 52-ой научной конференции МФТИ. Современные проблемы фундаментальных и прикладных наук. Часть 8: Проблемы современной физики. М.: МФТИ, 2009, С. 269-271.

2. Шашкин В.В. Толстых М.А. Полулагранжева модель мелкой воды на сфере на редуцированной сетке, сохраняющая массу // Избранные труды школы-конференции СИТЕС-2011. Томск: ЦНТИ, 2011, С. 57-60.

3. Tolstykh М., Shashkin V. Vorticity-divergence mass-conserving semi-Lagrangian shallow-water model using the reduced grid on the sphere // J. Comput. Phys. 2012. Vol. 231. P. 4205 - 4233. doi: 10.1016/j.jcp.2012.02.016.

4. Шашкин В.В. Локально-консервативный полулагранжев алгоритм численного решения уравнения переноса на сфере в трехмерном случае в z-системе координат по вертикали // Труды Гидрометцентра России. 2012. Вып. 348. С. 64-82. URL: http://method.hydromet.ru/publ/tr/tr348/shashkin.pdf.

5. Shashkin V., Tolstykh М. Inherently mass-conservative version of the semi-Lagrangian SL-AV atmospheric model dynamical core // Geosci. Model. Dev. Discuss. 2013. Vol. 6. P. 4809—4832. doi:10.5194/gmdd-6-4809-2013.

6. Lauritzen P., Ullrich P., Jablonowski C., Bosler P. A., Calhoun D., Conley A. J., Enomoto Т., Dong L., Dubey S., Guba O., Hansen А. В., Kaas E., Kent J., Lamarque J.-F., Prather M. J., Reinert D., Shashkin V. V., Skamarock W. C., S0rensen В., Taylor M. A., Tolstykh M. A. A standard test case suite for two-dimensional linear transport on the sphere: results from a collection of state-of-the-art schemes // Geosci. Model. Dev. Discuss. 2013. Vol. 6. P. 4983 — 5076. doi: 10.5194/gmdd-6-4983-2013.

Заказ № 56-а/11/2013 Подписано в печать 14.11.2013 Тираж 100 экз. Усл. п.л. 1,0

ООО "Цифровичок", тел. (495) 649-83-30 www.cfr.ru; е-тай:zak@cfr.ru

в

Федеральное государственное бюджетное учреждение науки Институт вычислительной

математики Российской академии наук

04201450933

На правах рукописи УДК 551.5

Шашкин Владимир Валерьевич

Полулагранжева модель динамики атмосферы мезомасштабного разрешения с использованием метода конечных объемов

Специальность 05.13.18 — «Математическое моделирование, численные методы и комплексы программ»

Диссертация на соискание учёной степени кандидата физико-математических наук

Научный руководитель: д.ф-м.н. Толстых М.А.

Москва - 2013

Содержание

Введение..................................................................................4

1 Конечно-объемные полулагранжев алгоритм численного решения уравнения переноса на сфере..................................................................................21

1.1 Конечно-объемный полулагранжев подход к решению одномерного уравнения переноса ........................................................................................24

1.2 Конечно-объемный полулагранжев алгоритм численного решения двумерного уравнения переноса на сфере на редуцированной широтно-долготной сетке .... 27

1.2.1 Описание алгоритма................................................................27

1.2.2 Результаты численного решения двумерного уравнения переноса на сфере . 31

1.3 Конечно-объемный полулагранжев алгоритм численного решения уравнения переноса на сфере в трехмерном случае....................................................37

1.3.1 Описание алгоритма................................................................37

1.3.2 Результаты численного решения трехмерного уравнения переноса на сфере 42

1.4 Основные результаты главы ..............................................................45

2 Полулагранжева полунеявная модель мелкой воды на сфере с использованием метода конечных объемов......................................................................47

2.1 Уравнения мелкой воды на вращающейся сфере........................................48

2.2 Полулагранжева полунеявная дискретизация уравнений мелкой воды................49

2.3 Численные эксперименты..................................................................51

2.3.1 Геострофическое равновесие......................................................51

2.3.2 Зональный поток над изолированной горой ....................................55

2.3.3 Квази-реальные данные............................................................55

2.3.4 Баротропная неустойчивость......................................................59

2.4 Основные результаты главы ..............................................................64

3 Полунеявная полулагранжева модель динамики атмосферы с конечно-объемной дискретизацией уравнений неразрывности и переноса водяного пара................65

3.1 Формулировка уравнений модели атмосферы ПЛАВ ..................................66

3.2 Дискретизация уравнений динамики атмосферы в модели ПЛАВ....................68

3.3 Версия динамического блока модели ПЛАВ с конечно-объемной полулагранжевой дискретизацией уравнения неразрывности и переноса водяного пара................70

3.4 Реализация модели на параллельных вычислительных системах......................72

3.5 Численные эксперименты..................................................................74

3.5.1 Эксперимент «Орографически возбужденные волны Россби»................76

3.5.2 Эксперимент «Бароклинная неустойчивость» ..................................78

3.5.3 Эксперимент «Моделирование циркуляции атмосферы на длительный промежуток времени»..................................................................85

3.5.4 Эксперимент «Моделирование атмосферной циркуляции на 72 часа» .... 89

3.6 Основные результаты главы ..............................................................91

Заключение..........................................................................................95

А Конечно-объемная аппроксимация горизонтальной дивергенции и оператора Лапласа ..............................................................................................98

Литература..........................................................................................99

Введение

Долгосрочный прогноз погоды, изучение глобальных климатических процессов и прогнозирование изменений климата невозможны без вычислительных экспериментов с математическими моделями климатической системы Земли. Модель климатической системы Земли современного уровня включает в себя модели общей циркуляции атмосферы и океана, модели динамики биосферы, морского и наземного льда, модели процессов на поверхности суши и в почве, а также модели взаимодействия на границах раздела всех перечисленных сред. Модели атмосферы и океана занимают центральное место в этом комплексе, так как процессы, формирующие климат, определяются, в основном, циркуляцией системы атмосфера, океан, почва, в то время как остальные компоненты климатической системы (биосфера, лед, поверхность суши, почва) влияют на эти процессы опосредованно, через взаимодействие с атмосферой и океаном.

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

Общая циркуляция атмосферы описывается системой уравнений, представляющей собой дифференциальную формулировку законов сохранения. Одна из возможных форм системы уравнений динамики атмосферы:

• Уравнения движения (закон сохранения импульса)

— + / - к х V = -УРФ +

(1)

• Уравнения притока тепла (первое начало термодинамики)

й% _ ИдТу Ы _

(И Срй[ 1 + (5 - 1)д] р

= Рт-

(2)

• Уравнение неразрывности (закон сохранения массы)

(3)

• Уравнение переноса водяного пара (сохранение массы водяного пара)

Данная система дополняется уравнениями гидростатики (используется гидростатическое приближение):

др = рдФ (5)

и уравнением состояния идеального газа

р = PRdTv. (6)

Здесь используются следующие обозначения V — (и, v) - горизонтальная скорость ветра, /

- параметр Кориолиса, р - давление, Vp - горизонтальный оператор набла при р = const, Ф

- геопотенциал (высота умноженная на g), ш - вертикальная скорость в р-системе координат, Tv = (1 4- [5 — 1 ]q)T - виртуальная температура, Т - температура, q - удельная влажность, 5 -отношение теплоемкостей сухого и влажного воздуха, Rd - газовая постоянная сухого воздуха, Cpd - теплоемкость сухого воздуха при постоянном давлении, F$, FT, Fq - источники и стоки физических величин вследствие процессов подсеточного масштаба.

Численные модели общей циркуляции атмосферы, как правило, состоят из двух блоков: динамическое ядро отвечает за решение системы уравнений (1-6), в блоке параметризованного описания физических процессов подсеточного масштаба рассчитываются источники (стоки) Ft, Fq. Параметризованное описание процессов подсеточного масштаба (конвекция, осадки, испарение, облака и др.) является одним из главных источников неопределенности моделирования и ошибок прогноза. При повышении разрешения некоторые процессы подсеточного масштаба начинают разрешаться на модельной сетке и становятся доступными для прямого моделирования. Естественным результатом является увеличение точности воспроизведения этих процессов.

На данный момент, перспективные версии ведущих моделей общей циркуляции, применяемых для климатических исследований и долгосрочного прогноза погоды (например, [1-4] и других) имеют разрешение 25-50 км, т.е. разрешают на сетке процессы масштабом 100-200 км (здесь учтено, что модель может эффективно разрешать процессы масштабом не менее 4 шагов сетки). Таким образом, современные модели климата имеют возможность прямого моделирования атмосферных процессов а-мезомасштаба (200 - 2000 км), таких, как фронты и ураганы [5], и вплотную подошли к моделированию явлений /3-мезомасштаба (20 - 200 км), таких, как бризы и тропические циклоны. Наибольшего успеха достигла японская группа разработчиков модели NICAM [3], достигшей горизонтального разрешения 3,5 км и способной разрешать на сетке крупные облачные скопления.

В России для целей моделирования климата применяется модель общей циркуляции атмосферы ИВМ РАН с горизонтальным разрешением около 150-200км [6], а для целей долгосрочного прогноза - модель ПЛАВ [7] с разрешением 100-150 км. Данные модели способны разрешать на сетке явления синоптического масштаба.

Повышение горизонтального разрешения моделей требует решения определенного количества проблем физического, математического и вычислительного плана. С точки зрения физики, основные проблемы - это переход от формулировки уравнений динамики атмосферы в гидростатическом приближении (1-6) к негидростатическим уравнениям, а также настройка блока описания процессов подсеточного масштаба с учетом того, что часть атмосферных процессов

начинает разрешаться на модельной сетке. Математические вопросы касаются выбора вычислительных сеток и методов дискретизации уравнений динамики атмосферы и отчасти взаимосвязаны с вопросами вычислительной эффективности и параллельной реализации моделей.

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

Одним из главных вопросов построения динамических блоков моделей общей циркуляции атмосферы является выбор вычислительной сетки на сфере. В настоящее время, наиболее распространены регулярные широтно-долготные сетки (см. Рисунок 1), состоящие из точек, лежащих на пересечениях меридианов и параллелей, взятых с постоянными шагами по широте и долготе. Популярность подобных сеток объясняется простой прямоугольной структурой в сферических координатах и, как следствие, простотой реализации горизонтальных дискретизаций. Недостатком регулярных сеток является сходимость меридианов к полюсам, вследствие чего, шаг сетки по долготе, выраженный в единицах длины, прогрессивно уменьшается в высоких широтах. Например, при разрешении по долготе около 10 км на экваторе, шаг по долготе около полюса будет в 100 раз меньше и составит 100 м, а при шаге на экваторе около 1 км, к полюсу шаг уменьшиться уже в 1000 раз (и составит всего лишь 1 м).

Рисунок 1: Сетки на сфере, слева направо: регулярная широтно-долготная сетка, редуцированная широтно-долготная сетка, треугольная икосаэдральная сетка, сетка «кубическая сфера». Рисунки взяты из [8]

Уменьшение шага регулярной сетки у полюсов сильно ограничивает шаг по времени в явных конечно-разностных моделях (вследствие условия Куранта около полюсов). При любом разумно малом, с точки зрения вычислительных затрат, шаге по времени подобные модели неустойчивы без применения полярных фильтров, удаляющих самые быстрые волны, уменьшая при этом точность решения. Кроме того, с точки зрения точности, увеличение разрешения около полюсов является излишним и приводит к неоправданным вычислительным затратам.

Интерес к моделированию атмосферы на сетках с квазиравномерным разрешением на сфере возник в 1960-х годах из-за нехватки вычислительных ресурсов для моделирования на регулярной сетке с достаточно малыми шагами по времени. К середине 1970-х годов, вследствие роста возможностей ЭВМ и совершенствования численных методов, этот интерес постепенно угас. В последнее время вопрос о моделировании атмосферы на сетках с квази-равномерным разрешением на сфере снова приобрел актуальность, так как модели атмосферы вплотную подошли к горизонтальному разрешению (порядка 1-10 км), при котором расчеты на регулярной сетке, из-за уменьшения шага по долготе у полюса, теряют смысл. Существует большое количество сеток с квазиравномерным разрешением на сфере, наиболее полный обзор которых приведен в [9], здесь же рассмотрены только основные варианты.

Самой простой сеткой с квази-равномерным разрешением является редуцированная широтно-долготная сетка [10], в которой, точки лежат на параллелях, взятых с постоянным шагом по широте (как в регулярной сетке), а шаг по долготе в угловых единицах растет при приближении к полюсам (см. Рисунок 1). Благодаря этому, в единицах длины, шаг по долготе уменьшается не так быстро или не уменьшается вовсе. Редуцированная сетка, в силу конструкции, в основном сохраняет все преимущества регулярной.

Редуцированные сетки активно применяются в спектральных моделях [11,12], так как, в силу свойств сферических функций, позволяют сократить вычислительные затраты (за счет уменьшения количества точек) без потери точности. Применение редуцированных сеток в конечно-разностных моделях было сочтено бесперспективным [9] из-за неудовлетворительной точности решения в высоких широтах. Данный вывод был поставлен под сомнение в работах [13], [14], показавших, что, при использовании специальным образом построенной редуцированной сетки [15], рост ошибок численного решения уравнения переноса и уравнений мелкой воды на сфере, при увеличении редукции до 20% (уменьшение количества узлов в редуцированной сетке относительно регулярной, с аналогичным разрешением на экваторе) остается незначительным.

Недостатком редуцированных широтно-долготных сеток является недостаточная равномерность разрешения (скачкообразное изменение разрешения от широты к широте) и анизотропность распределения точек, что отрицательно влияет на распространение волн. Более изотропное и равномерное распределение точек достигается в геодезических или икосаэдральных сетках (см. Рисунок 1). Подобные сетки основаны на делении треугольных граней икосаэдра, вписанного в сферу, на более мелкие треугольники и последующем гномоническом проецировании полученных треугольников на поверхность сферы. Существует так же вариант икосаэдральных сеток, основанный на делении граней икосаэдра на шестиугольники. Любая икосаэдральная сетка, основанная на шестиугольниках, содержит 12 пятиугольников.

Дискретизации уравнений динамики атмосферы на икосаэдральных сетках посвящено множество работ, из которых стоит выделить [16], в которой представлена глобальная оперативная модель численного прогноза погоды немецкой метеослужбы, и два варианта модели ICON (на треугольниках [2] и шестиугольниках [17]). Недостатком икосаэдральных сеток является их неортогональность, затрудняющая построение конечно-разностных и конечно-объемных дискретизаций высокого порядка. Кроме того, в силу структуры сетки (неортогональность, несиммет-

ричность, отличие соседних треугольников/шестиугольников по форме и ориентации, отсутствие выделенных направлений), могут возникать нефизические решения дискретных уравнений (вычислительные моды), например, ложная циркуляция масштаба ячейки сетки, нарушения reo строфического баланса и даже формирование ложных крупномасштабных циклонов. Подобные проблемы обозначаются понятием «отпечаток сетки» (от англ.- «grid imprinting»), большинство из них было решено в работах [17], [18].

На сетках типа «кубическая сфера» (см. Рисунок 1), построенных путем гномонической или конформной проекции регулярно распределенных точек сетки с граней вписанного куба на сферу, распределение точек менее изотропное и равномерное, зато более структурированное, чем на икосаэдральных сетках. Лучшая структурированность сетки упрощает дискретизацию уравнений. В качестве примеров моделей на сетке «кубическая сфера» можно привести [1,19,20]. В численных решениях уравнений динамики атмосферы на «кубической сфере» могут содержаться отпечатки сетки (хотя и в меньшей степени, чем на икосаэдральных сетках). Причинами появления отпечатков сетки является резкое изменение формы ячейки при переходе через ребро, а так же сложность построения точной дискретизации уравнений около ребер и вершин куба. Как правило, проблема исчезает с увеличением разрешения.

Еще одним вариантом вычислительной сетки, который удалось реализовать в модели атмосферы, является сетка Инь-Янь (см. Рисунок 2), представляющая объединение области регулярной широтно-долготной сетки от 45° ю.ш. до 45° с.ш. и от 0° до 270° в.д. с аналогичной областью, повернутой на 90°. В каждой из областей сетки Инь-Янь уравнения динамики атмосферы дискретизуются так же удобно и просто, как на регулярной сетке, однако, постановка граничных условий и дискретизация уравнений в областях стыковки и перекрытия двух областей вызывает трудности. Сетка Инь-Янь используется в перспективной модели метеослужбы Канады [21].

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

Вопрос выбора метода го�