Программа на Fortran90 для чтения файла NetCDF
В предыдущей заметке я описал разработку программы для создания файла NetCDF версии 3 с данными в узлах регулярной сетки. Теперь разработаем программу для чтения такого файла.
Описывается чтение файла формата NetCDF версии 3 с известной структурой. Указания соответствуют операционной системе Linux Debian, начиная с версии wheeze.
Считываемый файл создан на основе программ, описанных в предыдущей заметке, но для приближения к реальным условиям несколько усложнён:
- вместо абсолютных значений в файле содержатся относительные;
- вместо данных по всему Северному полушарию файл содержит данные по его части (прямоугольной трапеции).
Перед началом работы
Ознакомьтесь с моей заметкой, посвящённой созданию файла NetCDF, и выполните следующие действия:
- убедитесь, что на вашем компьютере установлены перечисленные в этой заметке программы;
- ознакомьтесь с разделом 1.2 («Reading a NetCDF Dataset with Known Names») документации на API;
- скомпилируйте программу netcdf_create.
Постановка задачи
Имеется файл NetCDF со значениями гидрометеорологического параметра kq в узлах регулярной сетки за несколько месяцев.
Необходимо вывести временной ряд среднемесячных значений kq в точке, координаты которой указываются в командной строке программы.
Формирование тестовых данных
Дамп NetCDF-файла с тестовыми данными (kq_region.cdl) приведён ниже:
netcdf kq_region {
dimensions:
lon = 9 ;
lat = 5 ;
time = UNLIMITED ; // (3 currently)
variables:
float lon(lon) ;
lon:standard_name = "longitude" ;
lon:long_name = "longitude" ;
lon:units = "degrees_east" ;
lon:axis = "X" ;
float lat(lat) ;
lat:standard_name = "latitude" ;
lat:long_name = "latitude" ;
lat:units = "degrees_north" ;
lat:axis = "Y" ;
double time(time) ;
time:standard_name = "time" ;
time:long_name = "time" ;
time:units = "day as %Y%m%d.%f" ;
time:calendar = "proleptic_gregorian" ;
time:axis = "T" ;
float kq(time, lat, lon) ;
kq:standard_name = "kq" ;
kq:long_name = "Sea Level KQ" ;
kq:units = "counts" ;
kq:add_offset = 40.f ;
kq:scale_factor = 1.f ;
kq:_FillValue = -999.9f ;
kq:missing_value = -999.9f ;
// global attributes:
:CDI = "Climate Data Interface version ?? (http://mpimet.mpg.de/cdi)" ;
:Conventions = "CF-1.6" ;
:history = "Wed Oct 24 12:38:08 2018: cdo sellonlatbox,60,100,40,60 kq2018.nc kq_region.nc\n",
"Wed Oct 24 12:38:08 2018: ncatted -O -a add_offset,kq,a,f,40. kq2018.nc\n",
"Wed Oct 24 12:38:08 2018: cdo mergetime jan.nc feb.nc mar.nc kq2018.nc\n",
"Wed Oct 24 12:38:08 2018 : ./netcdf_create jan.dat jan.nc 2018 01 15" ;
:source = "objective analyses of hydrometeorological fields" ;
:institution = "DataSource" ;
:title = "Database of SuperNII" ;
:references = "None" ;
:NCO = "\"4.6.3\"" ;
:CDO = "Climate Data Operators version 1.7.2 (http://mpimet.mpg.de/cdo)" ;
data:
lon = 60, 65, 70, 75, 80, 85, 90, 95, 100 ;
lat = 60, 55, 50, 45, 40 ;
time = 20180115, 20180215, 20180315 ;
kq =
21.06, 21.065, 21.07, 21.075, 21.08, 21.085, 21.09, 21.095, 21.1,
16.06, 16.065, 16.07, 16.075, 16.08, 16.085, 16.09, 16.095, 16.1,
11.06, 11.065, 11.07, 11.075, 11.08, 11.085, 11.09, 11.095, 11.1,
6.06, 6.065, 6.07, 6.075, 6.08, 6.085, 6.09, 6.095, 6.1,
1.06, 1.065, 1.07, 1.075, 1.08, 1.085, 1.09, 1.095, 1.1,
22.06, 22.065, 22.07, 22.075, 22.08, 22.085, 22.09, 22.095, 22.1,
17.06, 17.065, 17.07, 17.075, 17.08, 17.085, 17.09, 17.095, 17.1,
12.06, 12.065, 12.07, 12.075, 12.08, 12.085, 12.09, 12.095, 12.1,
7.06, 7.065, 7.07, 7.075, 7.08, 7.085, 7.09, 7.095, 7.1,
2.06, 2.065, 2.07, 2.075, 2.08, 2.085, 2.09, 2.095, 2.1,
23.06, 23.065, 23.07, 23.075, 23.08, 23.085, 23.09, 23.095, 23.1,
18.06, 18.065, 18.07, 18.075, 18.08, 18.085, 18.09, 18.095, 18.1,
13.06, 13.065, 13.07, 13.075, 13.08, 13.085, 13.09, 13.095, 13.1,
8.06, 8.065, 8.07, 8.075, 8.08, 8.085, 8.09, 8.095, 8.1,
3.06, 3.065, 3.07, 3.075, 3.08, 3.085, 3.09, 3.095, 3.1 ;
}
Тестовый NetCDF-файл можно сформировать из дампа следующей командой:
$ ncgen -o kq_region.nc kq_region.cdl
После создания тестового файла можно пропустить текст до конца раздела и перейти к следующему («Проектирование программы»). Тем не менее, опишем, как мы тестовый файл формировали.
Прежде всего, определимся с относительностью значений. Это означает, что в файле NetCDF хранятся не сами значения гидрометеорологического параметра, а информация для их вычисления по формуле y = scale_factor * x + add_offset.
Чаще всего, относительные значения используются в двух случаях:
- абсолютные значения настолько велики или настолько малы, что при оперировании с ними возможна потеря точности и искажение данных;
- в практической работе используются как относительное, так и абсолютное значения (например, вместо абсолютной температуры в NetCDF-файле удобно хранить температуру в градусах Цельсия и атрибут add_offset, равный 273,15).
В нашем тестовом файле scale_factor равен 1 (поэтому этот атрибут можно не приводить), а add_offset — 40.
Чтобы данные за разные месяцы различались, к значению в каждой точке добавим число, вводимое в командной строке.
Исходный код модифицированной команды testdata:
program testdata2
implicit none
character ( len = 80 ) :: outfile
character ( len = 5 ) :: chaddvalue
integer :: ilat, ilon
real, dimension ( 0:18, 0:71 ) :: array
real, parameter:: offset = 40.
real :: addvalue
call get_command_argument ( 1, outfile )
call get_command_argument ( 2, chaddvalue )
read (chaddvalue, '(f5.0)') addvalue
array ( 18, : ) = 90
do ilat = 17, 0, -1
do ilon = 0, 71
array (ilat, ilon ) = float (ilat * 5) + float (ilon * 5) / 1000.
end do
end do
open (1, file = outfile)
do ilat = 18, 0, -1
write (1, fmt = 800 ) ( array ( ilat, ilon) - offset + addvalue, ilon=0, 71 )
end do
close (1)
800 format (72f8.3)
stop 'Файл сформирован'
end program testdata2
Компилируем:
$ gfortran -Wall -o testdata2 testdata2.f90
Формируем файлы с помесячными данными:
$ ./testdata2 jan.dat 1.
$ ./testdata2 feb.dat 2.
$ ./testdata2 mar.dat 3.
На этом этапе вам необходимо решить, каким образом переменной kq NetCDF-файла присваивать атрибут «add_offset». Есть два варианта:
Первый — изменить программу netcdf_create, добавив в блок записи атрибутов гидрометеорологического параметра следующие строки:
status = nf90_put_att ( ncid = ncid, varid = hmparam_varid, &
name = 'add_offset', values = 40. )
call check_err ( status )
Второй — программу netcdf_create не изменять, а атрибут add_offset добавить в готовый NetCDF-файл.
Формируем файлы netcdf:
$ ./netcdf_create jan.dat jan.nc 2018 1 15
$ ./netcdf_create feb.dat feb.nc 2018 2 15
$ ./netcdf_create mar.dat mar.nc 2018 3 15
$ cdo mergetime jan.nc feb.nc mar.nc kq2018.nc
Если атрибут add_offset добавлялся не в программе netcdf_create, убедитесь, что на вашем компьютере установлен пакет NCO, и выполните команду:
$ ncatted -O -a add_offset,kq,a,f,40. kq2018.nc
Делаем выборку по району от 40° до 60° с. ш., от 60° до 100° в.д.:
$ cdo sellonlatbox,60,100,40,60 kq2018.nc kq_region.nc
Проверяем правильность получившегося файла. У переменной kq наименьшее значение должно быть в юго-западной точке района в январе (40 + 60/1000 + 1 = 41,060), максимальное — северо-восточной точке в марте (60 + 100/1000 + 3 = 63,100).
Выводим статистику:
$ cdo info kq_region.nc
Результат (лишние пробелы удалены):
-1 : Date Time Level Gridsize Miss : Minimum Mean Maximum : Parameter ID
1 : 2018-01-15 00:00:00 0 45 0 : 41.060 51.080 61.100 : -1
2 : 2018-02-15 00:00:00 0 45 0 : 42.060 52.080 62.100 : -1
3 : 2018-03-15 00:00:00 0 45 0 : 43.060 53.080 63.100 : -1
Проектирование программы
Необходимо проанализировать структуру файла kq_region.nc (выполните команду ncdump -c kq_region.nc и просмотрите её вывод). Выше приведён полный дамп этого файла, поэтому структуру не приводим.
На что обращаем внимание при анализе?
-
Убеждаемся, что координатами у переменной гидрометеорологического параметра являются время, широта и долгота. Из порядка следования координат получаем размерности массива для считывания параметра. Смотрим, какой тип данных у переменной параметра.
-
Убеждаемся, что переменная широты — одномерный массив, значения её — абсолютные, а единица измерения — «degrees_north». Смотрим, какой тип данных у переменной широты.
-
Убеждаемся, что переменная долготы — одномерный массив, значения её — абсолютные, а единица измерения — «degrees_east». Смотрим, какой тип данных у переменной долготы.
-
Убеждаемся, что переменная времени — одномерный массив, а единица её измерения — «day as %Y%m%d.%f». Смотрим, какой тип данных у переменной времени.
-
Проверяем, какую территорию покрывает наша сетка координат (земной шар, Северное или Южное полушария, или произвольную территорию.)
-
Проверяем, есть ли у гидрометеорологического параметра хотя бы один из атрибутов «scale_factor» или «add_offset».
-
Проверяем, есть ли у переменной гидрометеорологического параметра атрибут «missing_value».
-
Убеждаемся, что у переменной времени нет атрибута «time_bnds».
Для чего нужны эти проверки?
Пункты с первого по четвёртый определяют структуру данных, на которую рассчитана наша программа.
Кроме того, единица измерения времени позволяет понять, нужна ли специализированная библиотека подпрограмм для обработки времени.
Пятый пункт говорит о размерностях массивов в нашей программе: если территория — земной шар или полушария, можно использовать массивы фиксированного размера. Если территория — произвольная, целесообразно использовать динамически размещаемые массивы.
Шестой пункт показывает, требуется ли перевод относительных значений в абсолютные.
Последние два пункта связаны с особенностями программного обеспечения, сформировавшего обрабатываемый NetCDF-файл (оно может добавлять дополнительные атрибуты).
Так, седьмой пункт показывает, сколько условных значений обозначают отсутствующие данные. В конвенции ClimateForecast (CF) условное значение одно — значение атрибута _FillValue. В конвенции COARDS их два — значения атрибутов _FillValue и missing_value. То есть, атрибут _FillValue необходимо обрабатывать всегда, а во втором случае необходимо дополнительно обрабатывать и missing_value.
Восьмой пункт важен, если ваша программа считывает один файл формата NetCDF и формирует другой. Обычно у гидрометеорологических параметров в обоих файлах — одни и те же координаты, поэтому переменные координат (значения, размерности и атрибуты) удобно копировать из одного файла в другой без дополнительной обработки. В том числе, и переменную времени.
Однако программы, формирующие массивы реанализа NOAA, часто добавляют информацию о границах периодов, которым соответствуют значения переменной времени (например, для среднемесячных данных — значения первой и последней секунд месяцев), и значение атрибута time_bnds — имя переменной с такими границами.
В результате, программа должна реагировать на атрибут time_bnds у переменной времени: либо копировать в NetCDF-файл переменную с границами периодов, либо атрибут time_bnds в новый файл не выводить.
После анализа нашего файла проектируем программу так:
- Все массивы в программе будут динамически размещаемыми.
- Переменная гидрометеорологического параметра kq — трехмерная, тип её данных — REAL*4, описана как (time, lat, lon). Из этого следует, что массив для считывания переменной kq необходимо объявлять как real (lon, lat, time).
- При обработке значений kq будем использовать «рабочий» массив — с более привычной нам структурой real (time, lat, lon).
- Переменные широты и долготы нашим требованиям соответствуют: они — одномерные массивы типа REAL*4, значения координат — абсолютные.
- Переменная времени нашим требованиям тоже соответствует. Она — одномерный массив типа REAL*8, шкала времени — абсолютная, дополнительных библиотек для обработки времени не требуется.
- У переменной kq имеются атрибуты scale_factor и add_offset, поэтому добавим обработку этих двух атрибутов.
- Поставленная задача не требует обработки атрибута time_bnds, поэтому обработку этого атрибута не предусматриваем.
Дополнительные пояснения:
Во-первых, легко заметить, что программы CDO и NCO ко всем переменным нашего тестового файла добавили дополнительные атрибуты, в том числе, что важно, к переменной kq — атрибут scale_factor.
Поэтому обработку атрибутов scale_factor и add_offset рекомендуем добавлять всегда.
Во-вторых, для нашей задачи (вывод временного ряда в одной точке), строго говоря, достаточно массива «для считывания». Мы добавляем «рабочий» массив, чтобы приблизить учебную задачу к реальным условиям.
Вспомогательные процедуры
Для упрощения программы используем вспомогательные процедуры:
- для проверки результата операции с NetCDF-файлом (check_err);
- для определения позиции точки в массиве (location);
- для перевода переменной времени в дату (time2date).
Подпрограмму check_err возьмём из предыдущей заметки.
Функция location возвращает номер элемента массива, равного заданному значению. Номер элемента определяется перебором.
integer function location ( array, size, value )
implicit none
real, dimension ( size ), intent ( in ) :: array
integer, intent (in ) :: size
real, intent ( in ) :: value
integer :: i
location = -999
do i = 1, size
if ( array ( i ) == value ) then
location = i
exit
end if
end do
return
end function location
Подпрограмма time2date:
subroutine time2date ( value, year, month, day)
implicit none
real ( kind = 8 ), intent (in) :: value
integer, intent (out) :: year, month, day
real ( kind = 8 ) :: dt
year = int ( value / 1d4 )
dt = value - year * 1d4
month = int ( dt / 1d2 )
dt = dt - month * 1d2
day = int ( dt )
return
end subroutine time2date
Алгоритм работы программы
Пошагово опишем алгоритм с выписками из исходного текста (однотипные фрагменты кода не приводятся).
Вначале из командной строки считываем имя NetCDF-файла и координаты требуемой точки:
call get_command_argument ( 1, infile )
Широта и долгота преобразуются из символьной переменной в вещественные числа:
read ( chlon, '(f3.0)') lonreq
Открываем файл NetCDF:
status = nf90_open ( path = trim (infile), mode = nf90_nowrite, ncid = ncid )
call check_err ( status )
Считываем размерности переменных координат:
status = nf90_inq_dimid ( ncid = ncid, name ='lat', dimid = lat_dimid )
call check_err ( status )
Получаем идентификаторы переменных координат:
status = nf90_inq_varid ( ncid = ncid, name ='lat', varid = lat_varid )
call check_err ( status )
Получаем идентификатор переменной гидрометпараметра:
status = nf90_inq_varid ( ncid = ncid, name ='kq', varid = param_varid )
call check_err ( status )
Определяем размерности переменной гидрометпараметра и убеждаемся, что размерения следуют в правильном порядке (единственная проверка, без которой нельзя обойтись):
status = nf90_inquire_variable ( ncid = ncid, varid = param_varid, dimids = &
param_dimids )
call check_err ( status )
if ( param_dimids ( 1 ) /= lon_dimid .or. param_dimids ( 2 ) /= lat_dimid .or.&
param_dimids ( 3 ) /= time_dimid ) stop 'Неожиданная структура NetCDF-файла!'
Определяем размеры массивов координат:
status = nf90_inquire_dimension ( ncid = ncid, dimid = lat_dimid, len = numlat )
call check_err ( status )
На следующем этапе должны были считываться атрибуты координат (например, стандартное и полное имена переменных, единицы измерения и т. п.), но в нашей задаче эти атрибуты не используются.
Считываем атрибуты гидрометеорологического параметра, которые важны для оперирования данными (_FillValue, scale_factor, add_offset).
status = nf90_get_att ( ncid = ncid, varid = param_varid, name = '_FillValue',&
values = fillvalue )
if ( status /= nf90_noerr ) fillvalue = nf90_fill_real
Данная запись означает, что если у переменной атрибут _FillValue отсутствует, программа использует значение этого атрибута по умолчанию. Поскольку переменная kq имеет тип REAL*4, использовано значение по умолчанию для этого типа.
status = nf90_get_att ( ncid = ncid, varid = param_varid, name = &
'scale_factor', values = scalefactor )
if ( status /= nf90_noerr ) scalefactor = 1.
Данная запись означает, что если атрибут scale_factor не указан, программа использует его значение по умолчанию (при котором преобразование данных не требуется). То же указано для атрибута add_offset.
Размещаем в памяти динамические массивы для координат и «считываемого» массива гидрометпараметра:
allocate ( lat ( numlat ), lon ( numlon ), time ( numtime ), &
param_in ( numlon, numlat, numtime ) )
Строго говоря, можно было одновременно разместить в памяти все остальные массивы («рабочий» массив гидрометпараметра, массивы для номеров годов, месяцев, дат), но я привык размещать массивы непосредственно перед их использованием.
Считываем все переменные (и координаты, и гидрометеорологический параметр):
status = nf90_get_var ( ncid = ncid, varid = lat_varid, values = lat )
call check_err ( status )
Закрываем файл NetCDF:
status = nf90_close ( ncid )
call check_err ( status )
Нам удобнее, когда координаты у массива гидрометеорологического параметра идут в привычном порядке (время, широта, долгота), поэтому размещаем в памяти «рабочий» массив и наполняем его данными из «считываемого» массива:
do itime = 1, numtime
do ilat = 1, numlat
do ilon = 1, numlon
param ( itime, ilat, ilon ) = param_in ( ilon, ilat, itime )
end do
end do
end do
Если данные необходимо преобразовать из относительных значений в абсолютные (когда значение атрибута scale_factor не равно единице, а add_offset — нулю), выполняем это преобразование. Код для учёта scale_factor:
if ( scalefactor /= 1.0 ) then
do itime = 1, numtime
do ilat = 1, numlat
do ilon = 1, numlon
if ( param ( itime, ilat, ilon ) /= fillvalue ) &
param ( itime, ilat, ilon ) = param ( itime, ilat, ilon ) * scalefactor
end do
end do
end do
end if
Код для учёта add_offset — аналогичен.
Теперь у нас есть массив с реальными значениями гидрометеорологического параметра за несколько сроков.
По координатам требуемой точки определяем её позиции в массивах координат:
lonreqpos = location ( lon, numlon, lonreq )
Если точка не попадает ни в один из узлов нашей сетки координат, выводим сообщение об ошибке и прекращаем работу программы:
if ( latreqpos < 0 .or. lonreqpos < 0 ) stop 'Точка - вне района с данными!'
Преобразуем значения переменной времени в даты (в нашем примере годы, месяцы и дни размещаются в отдельных массивах):
do itime = 1, numtime
call time2date ( time ( itime ), datey ( itime ), datem ( itime ), dated &
( itime ))
end do
Наконец, выводим на экран временной ряд в нужной точке
do itime = 1, numtime
print '(i4,2(x,i2),f8.3)', datey ( itime ), datem (itime ), dated (itime), &
param ( itime, latreqpos, lonreqpos )
end do
Осталось освободить память от динамических массивов и завершить программу.
Полный исходный текст программы
program netcdf_load
use netcdf
implicit none
!
! Переменные, касающиеся программы в целом
!
character ( len = 80 ) :: infile
character ( len = 3 ) :: chlon
character ( len = 2 ) :: chlat
integer :: ncid, status
!
! Переменные, относящиеся к широте
real, dimension (:), allocatable :: lat
real :: latreq
integer :: lat_dimid, lat_varid, numlat
integer :: ilat, latreqpos
!
! Переменные, относящиеся к долготе
real, dimension (:), allocatable :: lon
real :: lonreq
integer :: lon_dimid, lon_varid, numlon
integer :: ilon, lonreqpos
!
! Переменные, относящиеся к времени
real ( kind = 8 ), dimension (:), allocatable :: time
integer, dimension ( : ), allocatable :: datey, datem, dated
integer :: time_dimid, time_varid, numtime
integer :: itime
!
! Переменные, относящиеся к гидрометпараметру
real, dimension ( :, :, : ), allocatable :: param_in, param
real :: scalefactor, offset, fillvalue
integer, dimension ( 3 ) :: param_dimids
integer :: param_varid
!
! Считывание и обработка командной строки
call get_command_argument ( 1, infile )
call get_command_argument ( 2, chlat )
call get_command_argument ( 3, chlon )
read ( chlat, '(f2.0)') latreq
read ( chlon, '(f3.0)') lonreq
!
! Открываем файл NetCDF
status = nf90_open ( path = trim (infile), mode = nf90_nowrite, ncid = ncid )
call check_err ( status )
!
! Получаем размерности переменных координат
status = nf90_inq_dimid ( ncid = ncid, name ='lat', dimid = lat_dimid )
call check_err ( status )
status = nf90_inq_dimid ( ncid = ncid, name ='lon', dimid = lon_dimid )
call check_err ( status )
status = nf90_inq_dimid ( ncid = ncid, name ='time', dimid = time_dimid )
call check_err ( status )
!
! Получаем идентификаторы переменных координат
status = nf90_inq_varid ( ncid = ncid, name ='lat', varid = lat_varid )
call check_err ( status )
status = nf90_inq_varid ( ncid = ncid, name ='lon', varid = lon_varid )
call check_err ( status )
status = nf90_inq_varid ( ncid = ncid, name ='time', varid = time_varid )
call check_err ( status )
!
! Получаем идентификатор переменной гидрометпараметра
status = nf90_inq_varid ( ncid = ncid, name ='kq', varid = param_varid )
call check_err ( status )
!
! Получаем и проверяем размерности переменной гидрометпараметра
status = nf90_inquire_variable ( ncid = ncid, varid = param_varid, dimids = &
param_dimids )
call check_err ( status )
if ( param_dimids ( 1 ) /= lon_dimid .or. param_dimids ( 2 ) /= lat_dimid .or.&
param_dimids ( 3 ) /= time_dimid ) stop 'Неожиданная структура NetCDF-файла!'
!
! Определяем размеры массивов координат
status = nf90_inquire_dimension ( ncid = ncid, dimid = lat_dimid, len = numlat )
call check_err ( status )
status = nf90_inquire_dimension ( ncid = ncid, dimid = lon_dimid, len = numlon )
call check_err ( status )
status = nf90_inquire_dimension ( ncid = ncid, dimid = time_dimid, len = numtime )
call check_err ( status )
!
! Считываем необходимые для оперирования с данными атрибуты гидрометпараметра
status = nf90_get_att ( ncid = ncid, varid = param_varid, name = '_FillValue',&
values = fillvalue )
if ( status /= nf90_noerr ) fillvalue = nf90_fill_real
status = nf90_get_att ( ncid = ncid, varid = param_varid, name = &
'scale_factor', values = scalefactor )
if ( status /= nf90_noerr ) scalefactor = 1.
status = nf90_get_att ( ncid = ncid, varid = param_varid, name = 'add_offset',&
values = offset)
if ( status /= nf90_noerr ) offset = 0.
!
! Размещаем массивы для считывания в памяти
allocate ( lat ( numlat ), lon ( numlon ), time ( numtime ), &
param_in ( numlon, numlat, numtime ) )
!
! Считываем все переменные
status = nf90_get_var ( ncid = ncid, varid = lat_varid, values = lat )
call check_err ( status )
status = nf90_get_var ( ncid = ncid, varid = lon_varid, values = lon )
call check_err ( status )
status = nf90_get_var ( ncid = ncid, varid = time_varid, values = time )
call check_err ( status )
status = nf90_get_var ( ncid = ncid, varid = param_varid, values = param_in )
call check_err ( status )
!
! Закрываем файл NetCDF
status = nf90_close ( ncid )
call check_err ( status )
!
! Переводим массив параметра в привычный вид
allocate ( param ( numtime, numlat, numlon ) )
do itime = 1, numtime
do ilat = 1, numlat
do ilon = 1, numlon
param ( itime, ilat, ilon ) = param_in ( ilon, ilat, itime )
end do
end do
end do
!
! Учитываем scale_factor
if ( scalefactor /= 1.0 ) then
do itime = 1, numtime
do ilat = 1, numlat
do ilon = 1, numlon
if ( param ( itime, ilat, ilon ) /= fillvalue ) &
param ( itime, ilat, ilon ) = param ( itime, ilat, ilon ) * scalefactor
end do
end do
end do
end if
!
! Учитываем add_offset
if ( offset /= 0. ) then
do itime = 1, numtime
do ilat = 1, numlat
do ilon = 1, numlon
if ( param ( itime, ilat, ilon ) /= fillvalue ) &
param ( itime, ilat, ilon ) = param ( itime, ilat, ilon ) + offset
end do
end do
end do
end if
!
! Определяем позиции точки в массивах по её координатам
latreqpos = location ( lat, numlat, latreq )
lonreqpos = location ( lon, numlon, lonreq )
if ( latreqpos < 0 .or. lonreqpos < 0 ) stop 'Точка - вне района с данными!'
!
! Переводим переменные времени в даты
allocate ( datey ( numtime ), datem ( numtime ), dated ( numtime ) )
do itime = 1, numtime
call time2date ( time ( itime ), datey ( itime), datem ( itime ), dated &
( itime ))
end do
!
! Выводим временной ряд в нужной точке
do itime = 1, numtime
print '(i4,2(x,i2),f8.3)', datey ( itime ), datem (itime ), dated (itime), &
param ( itime, latreqpos, lonreqpos )
end do
!
! Завершение работы программы
deallocate ( lat, lon, time, datey, dated, datem, param_in, param )
stop 'Программа доработала до конца'
!
contains
!
subroutine check_err ( stat )
use netcdf
implicit none
integer :: stat
if (stat /= NF90_NOERR) then
print *, nf90_strerror ( stat )
stop
endif
return
end subroutine check_err
!
integer function location ( array, size, value )
implicit none
real, dimension ( size ), intent ( in ) :: array
integer, intent (in ) :: size
real, intent ( in ) :: value
integer :: i
location = -999
do i = 1, size
if ( array ( i ) == value ) then
location = i
exit
end if
end do
return
end function location
!
subroutine time2date ( value, year, month, day)
implicit none
real ( kind = 8 ), intent (in) :: value
integer, intent (out) :: year, month, day
real ( kind = 8 ) :: dt
year = int ( value / 1d4 )
dt = value - year * 1d4
month = int ( dt / 1d2 )
dt = dt - month * 1d2
day = int ( dt )
return
end subroutine time2date
!
end program netcdf_load
Компиляция и отладка
Вначале готовим данные для отладки.
Выбираем точку (допустим, 50° с. ш. 90° в. д.):
$ cdo sellonlatbox,90,90,50,50 kq_region.nc test.nc
Выводим временной ряд в этой точке:
$ cdo outputtab,year,month,day,value test.nc
Получаем:
# year month day value
2018 1 15 51.09
2018 2 15 52.09
2018 3 15 53.09
Этот же результат должна вывести и наша программа.
Компилируем её:
$ gfortran netcdf_load.f90 -o netcdf_load -l netcdff -I/usr/include
Запускаем:
$ ./netcdf_load kq_region.nc 50 90
Вывод программы:
2018 1 15 51.090
2018 2 15 52.090
2018 3 15 53.090
STOP Программа доработала до конца
То есть, программа работала правильно.
Заключение
В заметке описана разработка программы, пригодной для учебных целей. Однако приведены приёмы программирования, используемые в практической работе с данными (использование динамических массивов и преобразование данных).
Программа максимально упрощена — исключены проверки входных данных, освобождение памяти от динамически размещаемых массивов вынесено в конец программы, и т. п.
Надеюсь, этот текст будет полезен осваивающим программирование для формата NetCDF.
Автор: Георгий Мурый, site@moury.ru, 26 ноября 2018 года.