CC BY
50
Вычислительные технологии
Том 11, часть 1, Специальный выпуск, 2006
ПОЛУЛАГРАНЖЕВА МОДЕЛЬ ПРОГНОЗА ПОГОДЫ С ПЕРЕМЕННЫМ РАЗРЕШЕНИЕМ И ЕЕ ДАЛЬНЕЙШЕЕ РАЗВИТИЕ*
М. А. Толстых Институт вычислительной математики РАН и Гидрометцентр России, Москва e-mail: tolstykh@inm.ras.ru
Р.Ю. Фадеев Московский физико-технический институт, Россия e-mail: Lord.EED@gmail.com
A global semi-Lagrangian finite-difference variable-resolution numerical weather prediction model was developed in INM RAS and Hydrometheorologieal center of Russia for 3 days forecasts. Precipitation parameterization schemes used in the model are described. The results of meteoelements fields forecasts at upper atmosphere and some preliminary results of precipitation forecast tests are presented. Two-dimensional nonhydrostatic version of the dynamical core have been developed on the base of the given model. The results of the commonly used test describing stationary flow over the mountain are given.
Введение
В Институте вычислительной математики РАН и Гидрометцентре России разработана глобальная полулагранжева конечно-разностная модель прогноза погоды в версии с переменным разрешением (ПЛАВ-ПР), предназначенная для прогноза с заблаговременностью до трех суток. Модель с переменным разрешением создана на основе версии той же модели с постоянным разрешением ill. Подробное описание динамического блока приведено в [2, 3]. Отличительными особенностями данной модели являются применение компактных разностей четвертого порядка для аппроксимации неадвективных слагаемых и использование вертикальной компоненты абсолютного вихря и горизонтальной дивергенции в качестве прогностических переменных. Численные алгоритмы, примененные в модели, подробно описаны в [4]. В модели используется набор параметризаций процессов подсе-точного масштаба из французской оперативной модели ARPEGE/IFS [5], который адаптирован к применению переменного разрешения. Программный комплекс модели успешно реализован на параллельных вычислительных системах с общей и распределенной памятью при использовании технологий распараллеливания ОрепМР и MPI,
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 04-05-64638).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.
Разрешение модели по широте изменяется от 30 км в полосе широт 48-90° с.ш. до 70 км в Южном полушарии. Шаг сетки по долготе составляет 0.5625°, и за счет сходимости меридианов к полюсу горизонтальное разрешение над Россией составляет около 30 км. Модель имеет 28 неравномерно расположенных сигма-уровней по вертикали. Это первая отечественная модель прогноза погоды с переменным разрешением. Также впервые в отечественной практике достигнуто горизонтальное разрешение в 30 км над всей территорией России.
Модель проходила квазиоперативные испытания в Гидрометцентре России с декабря 2004 по август 2005 гг. В качестве начальных данных использовались результаты анализа системы усвоения данных (СУД-ИОИ) Гидрометцентра России [6] на основе версии данной модели с постоянным разрешением 0.9 град, по долготе, 0.72 град, по широте, которые сначала интерполировались по горизонтали на сетку модели с переменным разрешением, а затем — по вертикали для учета разницы орографий в двух версиях модели. Вследствие этого начальные данные содержали ошибку, связанную с интерполяцией. Оценки качества прогнозов некоторых полей метеоэлементов в свободной атмосфере (на сетке 2.5 х 2.5 град, по стандарту Всемирной метеорологической организации) по региону Европа для версий модели с постоянным и переменным разрешением приведены на рис. 1, 2. По региону Азия модель с переменным разрешением в указанный период показала несколько мень-
Рис. 1. Градиентная ошибка геопотенциала 500 гПа (снизу) и среднеквадратичная ошибка температуры на уровне 850 гПа (вверху) для прогнозов на срок 48 ч (слева) и 24 ч (справа) для полулагранжевой модели с постоянным (БЬМ) и переменным (8ЬМ-Уаг) разрешением для каждого из месяцев в период с декабря 2004 по август 2005 гг.
0,0
123456789 1 23456789
Рис. 2. Среднеквадратичная ошибка геопотенциала 500 гПа для прогнозов на срок 72 ч (слева) и 24 ч (справа) для полулагранжевой модели с постоянным (ЯЬМ) и переменным (ЯЬМ-Уаг) разрешением для каждого из месяцев в период с декабря 2004 по август 2005 гг.
шую успешность прогнозов. В целом оказалось, что полулагранжева модель с переменным разрешением обеспечивает примерно такое же качество прогноза, что и полулагранжева модель с постоянным разрешением. Испытания показали в том числе пути улучшения результатов, в частности усовершенствование процедуры подготовки начальных данных.
1. Прогноз осадков
Рассмотрим сначала схемы параметризации осадков. В модели применяются параметризации из модели Метео-Франс ARPEGE/IFS [5] с модифицированным (гибридным) замыканием для глубокой конвекции [7].
Вычисление крупномасштабных осадков основано на диагностически определенном влагосодержании с использованием функции распределения размера капель Маршала — Палмера и предписанной концентрации капель. Скорость падения капель зависит от их диаметра. Таяние либо замерзание выпадающих осадков происходит в параметризации в то же самое время, что и их выпадение. Все перенасыщение удаляется как выпадение воды или льда. Различие между жидкой и твердой фазами воды состоит в следующем:
— для термодинамических эффектов фазовых переходов налагается ступенчатый переход осадков из одной фазы в другую в тройной точке для согласованности с функциями насыщения;
— для скорости испарения и/или таяния-замерзания, которая зависит от скорости выпадения, при генерации осадков делается различие между фазами воды, которое зависит только от температуры.
Коэффициент испарения осадков в ненасыщенных нижележащих уровнях вычисляется по модифицированной формуле Кесслера.
Параметризация глубокой конвекции выполнена на основе схемы [8] с многочисленными усовершенствованиями, реализованными в Метео-Франс за почти 20 лет существования этой схемы. Скорость выпадения конвективных осадков определяется разностью между общей конвергенцией влаги и увлажнением окружения благодаря перемешиванию с ним облака. Под основанием облака возможно испарение осадков. Облачный конденсат преобразуется в осадки, когда эквивалентная толщина насыщенного облачного слоя
превосходит некую критическую величину. Пропорция снега в осадках на данном вертикальном уровне зависит от пропорции снега на предыдущем уровне и температуры. В отличие от оригинальной схемы, введено изменение скорости вовлечения с высотой, которая равна максимальному значению Епггх в основании облака и затухает по экспоненте к стандартному значению Епгг по мере подъема облака. Значения Епгг и Епггх зависят от интегральной плавучести облака таким образом, что в глубоких облаках (тропики) скорость вовлечения меньше, чем в неглубоких (средние широты). В схеме учитываются нисходящие потоки по краям облака, а также перераспределение момента вследствие конвекции [9].
В отличие от оригинальной схемы, где применяется замыкание либо типа Куо, либо на основе конвективной доступной потенциальной энергии (САРЕ), используется замыкание САРЕ, когда температура на нижнем модельном уровне не превышает примерно 29 °С, в остальных случаях — замыкание типа Куо [7].
В настоящий момент доступны оценки осадков за период с 01.07 по 24.09 2006 г. Сравнивались две версии региональной модели ММ5, работающие в опытном режиме в Гид-
Рис. 3. Различные характеристики качества прогноза факта осадков для модели ММ5, работающей в Гидрометцентре РФ (ММ5-1), модели ММ5, работающей в Гидрометеобюро Москвы и Московской обл. (ММ5-2), и полулагранжевой модели с переменным разрешением (ЭГАУ-УК) для сроков прогноза 12, 24, 36 и 48 ч.
рометцентре РФ (с горизонтальным разрешением над Центрально-Европейским регионом Россиеи 18 км) и Московском гидрометеобюро (с разрешением 15 км), а также глобальная полулагранжева модель с переменным разрешением ПЛАВ-ПР (разрешение над данным регионом 30 км). Обе версии ММ5 в качестве начальных данных использовали анализы МСЕР, модель ПЛАВ-ПР стартовала с анализа системы усвоения данных Гидрометцентра РФ, Сравнение должно рассматриваться только как демонстрация возможностей модели с переменным разрешением, так как оценки доступны пока за небольшой период.
Известно, что разрешение модели играет существенную роль для прогноза осадков. Приведенные на рис, 3 различные виды оценок прогноза факта осадков по Центрально-Европейской части России показывают, что модель с переменным разрешением ПЛАВ-ПР способна давать конкурентоспособный прогноз осадков при более грубом по сравнению с моделью ММ5 горизонтальном разрешении. Немаловажным фактором здесь является вычислительная эффективность модели — прогноз на одни сутки рассчитывается 38 мин на четырех процессорах Итаниум-2, в то время как ММ5 требует около 1 ч на девяти процессорах,
2. Негидростатический динамический блок
Увеличение пространственного разрешения моделей прогноза погоды позволяет получать более точный прогноз. Однако при повышении горизонтального разрешения до примерно 10 км приближение гидростатики, предполагающее малость вертикальных ускорений по сравнению с ускорением свободного падения, становится неточным, особенно вблизи гор. При использовании высокопроизводительных компьютеров открывается возможность оперативного применения региональных моделей с таким разрешением уже в наши дни, В ближайшем будущем такое разрешение станет возможным и для глобальных моделей. Таким образом, весьма актуальным становится переход на полные негидростатические уравнения гидротермодинамики атмосферы,
В настоящее время существует целый ряд негидростатических моделей численного прогноза погоды [10], Обзор негидростатических моделей прогноза погоды по состоянию на 2002 г. можно найти, например, в [11]. Часть пз нпх основана на полулагранжевом полунеявном подходе [12, 13], Это, например, региональные модели А1ас1ш [14] и НШ1..\.\1 [15], К классу эйлеровых относятся региональные модели XVНК [16] и !..\1 [17].
Модели, использующие полунеявный подход [13], предполагают задание фонового состояния атмосферы, с помощью которого уравнения гидротермодинамики разделяются на две части: линейную, интегрируемую по времени по неявной схеме, и нелинейную, трактуемую явно. Устойчивость алгоритма зависит в первую очередь от того, насколько малы отклонения, интегрируемые по явной схеме. Поэтому полунеявный подход требует весьма аккуратного подхода к выбору состояния атмосферы, относительно которого производится линеаризация уравнений (фонового состояния). В качестве фонового состояния в некоторых случаях выбирают решение гидростатических уравнений. Фоновое состояние может пересчитываться на каждом шаге по времени, что часто необходимо для уменьшения нормы нелинейных слагаемых, интегрируемых по явной схеме, а следовательно, повышения устойчивости.
В данной работе реализована двумерная версия негидростатического динамического блока модели атмосферы на основе сочетания полулагранжевой полунеявной гидростатической модели ПЛАВ [1] и негидростатического варианта модели НШ1..\.\1 [15, 18-20].
Выбор модели ШШ.ЛМ обусловлен, в первую очередь, применением в ней вычислительно эффективного полулагранжева полунеявного подхода к решению уравнений гидротермодинамики атмосферы. Такой же подход принят и в исходной гидростатической модели ПЛАВ, Использование приближения, при котором негидростатическая добавка к давлению на поверхности земли обращается в нуль, позволяет относительно легко модифицировать алгоритм, разработанный применительно к вертикальной гибридной р-коордипате, для а-коордпнаты, используемой в модели ПЛАВ, Алгоритм, применяемый в модели ШШ.ЛМ. предполагает равенство нулю трехмерной дивергенции, тем самым исключаются из рассмотрения метеорологически незначимые звуковые волны. Кроме того, на основе масштабного анализа (приводимого в [19]) упрощается уравнение для вертикальной скорости.
При реализации динамического блока модели ШШ.ЛМ оригинальный алгоритм был переработан, хотя основная его структура осталась неизменной. Изменения коснулись только той части алгоритма, где модификации были необходимы в силу принципиального различия двух моделей. Проведена работа по увеличению устойчивости и улучшению аппроксимационных свойств. Отметим ряд ключевых отличий оригинального алгоритма.
1. Для расчета горизонтальных производных применяются не только центральные разности, как это сделано в модели ШШ.ЛМ. но и компактные схемы четвертого порядка точности [4].
2. Численные эксперименты показали, что используемый в ШШ.ЛМ дискретный оператор интегрирования по вертикали (применяется формула средней точки) для расчета, в частности, геопотенциала характеризуется недостаточной точностью. Использование формулы трапеций при расчете дискретного интеграла по высоте позволяет существенно повысить точность расчета флуктуационной части геопотенциала.
3. Оригинальный алгоритм ШШ.ЛМ предполагает использование неравномерной смещенной сетки по горизонтали и равномерной по вертикали. В реализованном двумерном блоке, напротив, применяется равномерная несмещенная сетка по горизонтали и неравномерная по вертикали.
В качестве теста двумерного варианта динамического блока моделировалось обтекание горы стационарным набегающим потоком. В рассматриваемом тесте во всей области интегрирования профили температуры Т = 280 К и горизонтальной скорости и = 30 м/с выбирались постоянными. Начальное значение вертикальной скорости рассчитывалось с применением уравнения неразрывности. Орография задавалась следующей формулой:
н
Ых) = --:-:-—. (1)
1 у 1+(х/ахУ к ;
Максимальная высота горы Н составляла 100 м , а полуширина ах = 3 км. Давление на поверхности земли в начальный момент времени определялось выражением
к
о
dz RI\z)
ps = p°s exp -g / , (2)
где фоновое значение давления выбиралось равным pO = 101 325 Па. Шаг интегрирования по времени составлял 30 с. Шаг сетки по горизонтали dx = 528 м, по вертикали использовалось равномерное разбиение а-коордпнаты на 100 интервалов.
На рис. 4 показаны линии уровня и (вертикальная скорость в p-системе координат: и = dp/dt) после 250 шагов интегрирования по времени. Можно видеть, что в результате
О 23.5 47.0 70.5 94.0 0 23.5 47.0 70.5 94.0
-2 -1.5 -1.0 -0.5 0 0.5 1.0 1.5 2 -2 -1.5 -1.0 -0.5 0 0.5 1.0 1.5 2
а б
Рис. 4. Вертикальная скорость и в p-системе координат (а) и отклонение температуры от среднего значения (б) в расчетной области (по горизонтали — ж-координата, км, по вертикали — безразмерная а-координата).
процесса обтекания горы (на рисунках движение происходит слева направо) с ее тыльной стороны формируются распространяющиеся вправо волны. Подобные волны являются негидростатическими и могут быть получены численно только при решении негидростатических уравнений гидротермодинамики атмосферы. В целом можно констатировать качественное согласие с результатами других авторов (например, [21,22]). Полученная амплитуда волн также близка к результатам упомянутых работ. Однако в верхней части модельной атмосферы волны заметно слабее, что, возможно, вызвано неоптимальными параметрами модели. В настоящее время ведутся работы по проверке реализованного динамического блока на других тестах, а также настройке модельных параметров.
С помощью реализованного негидростатического динамического блока можно оценить диапазон изменения различных параметров атмосферы, что позволит в дальнейшем правильно выбрать как прогностические переменные, так и фоновое состояние атмосферы. Реализованный динамический блок послужит основой для разработки оригинального негидростатического трехмерного динамического блока модели ПЛАВ.
Список литературы
[1] Толстых М.А. Полулагранжева модель атмосферы с высоким разрешением для численного прогноза погоды // Метеорология и гидрология. 2001. № 4. С. 5-16.
[2] Tolstykh М.А. Variable resolution global semi-Lagrangian atmospheric model // Russ J. Num. An. к Math. Mod. 2003. Vol. 18, N 4. P. 347-361.
[3] Толстых М.А. Модель прогноза погоды с переменным разрешением // Вычисл. технологии. 2004. Т. 9. Спецвыпуск: Тр. Междунар. конф. "Вычислительно-информационные технологии для наук об окружающей среде" — CITES-2003. Ч. 1. С. 86-94.
[4] Tolstykh М. Vorticity-divergence semi-lagrangian shallow-water model on the sphere based on compact finite differences // J. Comput. Phys. 2002. Vol. 179. P. 180-200.
[5] Geleyn J.-F., Bâzile E., Bougeault P. et al. Atmospheric parameterization schemes in Meteo-Franee's ARPEGE N.W.P. model // Parameterization of Subgrid-seale Physical Processes. ECMWF Seminar Proc. Reading, UK, 1994. P. 385-402.
[6] Цырульников M.Д., Толстых M.А., Багров А.Н., Зарипов Р.Б. Развитие глобальной системы усвоения данных с переменным разрешением // Метеорология и гидрология. 2003. № 4. С. 5-23.
[7] Tolstykh M. The use of combined closure in convection parameterization scheme // WMO/WCRP/WGNE Research Activities in Atmodpheric and Oceanic Modelling / J. Côté (Ed.) WMO, Geneva. 2003. Rep. N 33. P. 04.25-04.26.
[8] Bougeault P. A simple parameterization of the large-scale effects of cumulus convection // Mon. Weather Rev. 1985. Vol. 113. P. 2108-2121.
[9] Gregory D., Kershaw R., Inness P.M. Parameterization of momentum transport by convection. Pt IP Tests in single-column and general circulation models // Quart. J. Roy. Meteor. Soc. 1997. Vol. 123. P. 1153-1183.
flOj Толстых M.А., Фролов А.В. Некоторые проблемы численного прогноза погоды // Изв. РАН. Физика атмосферы и океана. 2005. Т. 41. С. 315-327.
[11] Steppeler J., Hess R., Schattler U., Bonaventura L. Review of numerical methods for nonhydrostatic weather prediction models // Meteor. Atmos. Phys. 2003. Vol. 82, N 1-4. P. 287301.
[12] Staniforth A., Cote J. Semi-Lagrangian integration schemes for atmospheric models. A review // Mon. Wea. Rev. 1991. Vol. 119. P. 2206-2223.
[13] Benard P. On the use of a wider class of linear systems for the design of constant-coefficients semi-implicit time schemes in NWP // Mon. Wea. Rev. 2004. Vol. 132. P. 1319-1324.
[14] Geleyn J.-F. The ALADIN Project: mesoscale operational forecasting and research tool from a new generation // Proc. of the 24th ICAM, Bled (Slovenia), Sept. 9-13. 1996. P. 85-92.
[15] room R., Mannik A. Non-hydrostatic Adiabatic Kernel for HIRLAM. Pt IV: Semi-implicit semi-Lagrangian scheme. HIRLAM Technical Rep. 65. Feb. 2006. http://hirlam.org
[16] Michalakes J., Chen S., Dudhia J. et al. Development of a next-generation regional weather research and forecast model // Developments in Teracomputing: Proc. of the 9th ECMWF Workshop on the Use of High Performance Computing in Meteorology / W. Zwiehofer, N. Kreitz (Ed.). Singapore: World Scientific, 2001. P. 269-276.
^lTj Steppeler JT.^ Donis Cj.^ Schattler, TJ. et al. Meso-gamma scale forecasts using the nonhydrostatic model LM // Meteor. Atmos. Phys. 2003. Vol. 82, N 1-4. P. 75-96.
[18] room R. Nonhydrostatic adiabatic kernel for HIRLAM. Pt 1: Fundamentals of nonhydrostatic dynamics in pressure-related coordinates. HIRLAM Technical Rep. 48. Feb. 2001.
http://hirlam.org/
[19] Mannik A., Room R. Nonhydrostatic adiabatic kernel for HIRLAM. Pt II: Anelastic, Hybrid-coordinate, Explicit-Eulerian Model. HIRLAM Technical Rep. 49. May 2001.
http://hirlam.org/
184
M.A. Tojicthx, P.K). (PafleeB
[20] room R., Mannik A. Non-hydrostatic Adiabatic Kernel for HIRLAM. Pt III: Semi-implicit Eulerian Scheme. HIRLAM Technical Rep. 55, Aug. 2002. http://hirlam.org/
[21] Room R.. Mannik A. Response of different nonhydrostatic, pressure-coordinate models to orographic forcing //J. Atmos. Sei. 1999. Vol. 56, N 15. P. 2553-2570.
[22] Janjic Z.I., Gerrity J.P., Nickovic S. An alternative approach to nonhydrostatic modeling // Mon. Wea. Rev. 2001. Vol. 129, N. 5. P. 1164-1178.
Ilocmynu/ia e pedanv^wo 9 noMÖpM 2006 '¿.
