Сайт Георгия Мурого

Программа на Fortran90 для чтения файла NetCDF

В предыдущей заметке я описал разработку программы для создания файла NetCDF версии 3 с данными в узлах регулярной сетки. Теперь разработаем программу для чтения такого файла.

Описывается чтение файла формата NetCDF версии 3 с известной структурой. Указания соответствуют операционной системе Linux Debian, начиная с версии wheeze.

Считываемый файл создан на основе программ, описанных в предыдущей заметке, но для приближения к реальным условиям несколько усложнён:

Перед началом работы

Ознакомьтесь с моей заметкой, посвящённой созданию файла NetCDF, и выполните следующие действия:

Постановка задачи

Имеется файл 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.

Чаще всего, относительные значения используются в двух случаях:

В нашем тестовом файле 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 и просмотрите её вывод). Выше приведён полный дамп этого файла, поэтому структуру не приводим.

На что обращаем внимание при анализе?

  1. Убеждаемся, что координатами у переменной гидрометеорологического параметра являются время, широта и долгота. Из порядка следования координат получаем размерности массива для считывания параметра. Смотрим, какой тип данных у переменной параметра.

  2. Убеждаемся, что переменная широты — одномерный массив, значения её — абсолютные, а единица измерения — «degrees_north». Смотрим, какой тип данных у переменной широты.

  3. Убеждаемся, что переменная долготы — одномерный массив, значения её — абсолютные, а единица измерения — «degrees_east». Смотрим, какой тип данных у переменной долготы.

  4. Убеждаемся, что переменная времени — одномерный массив, а единица её измерения — «day as %Y%m%d.%f». Смотрим, какой тип данных у переменной времени.

  5. Проверяем, какую территорию покрывает наша сетка координат (земной шар, Северное или Южное полушария, или произвольную территорию.)

  6. Проверяем, есть ли у гидрометеорологического параметра хотя бы один из атрибутов «scale_factor» или «add_offset».

  7. Проверяем, есть ли у переменной гидрометеорологического параметра атрибут «missing_value».

  8. Убеждаемся, что у переменной времени нет атрибута «time_bnds».

Для чего нужны эти проверки?

Пункты с первого по четвёртый определяют структуру данных, на которую рассчитана наша программа.

Кроме того, единица измерения времени позволяет понять, нужна ли специализированная библиотека подпрограмм для обработки времени.

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

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

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

Так, седьмой пункт показывает, сколько условных значений обозначают отсутствующие данные. В конвенции ClimateForecast (CF) условное значение одно — значение атрибута _FillValue. В конвенции COARDS их два — значения атрибутов _FillValue и missing_value. То есть, атрибут _FillValue необходимо обрабатывать всегда, а во втором случае необходимо дополнительно обрабатывать и missing_value.

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

Однако программы, формирующие массивы реанализа NOAA, часто добавляют информацию о границах периодов, которым соответствуют значения переменной времени (например, для среднемесячных данных — значения первой и последней секунд месяцев), и значение атрибута time_bnds — имя переменной с такими границами.

В результате, программа должна реагировать на атрибут time_bnds у переменной времени: либо копировать в NetCDF-файл переменную с границами периодов, либо атрибут time_bnds в новый файл не выводить.

После анализа нашего файла проектируем программу так:

  1. Все массивы в программе будут динамически размещаемыми.
  2. Переменная гидрометеорологического параметра kq — трехмерная, тип её данных — REAL*4, описана как (time, lat, lon). Из этого следует, что массив для считывания переменной kq необходимо объявлять как real (lon, lat, time).
  3. При обработке значений kq будем использовать «рабочий» массив — с более привычной нам структурой real (time, lat, lon).
  4. Переменные широты и долготы нашим требованиям соответствуют: они — одномерные массивы типа REAL*4, значения координат — абсолютные.
  5. Переменная времени нашим требованиям тоже соответствует. Она — одномерный массив типа REAL*8, шкала времени — абсолютная, дополнительных библиотек для обработки времени не требуется.
  6. У переменной kq имеются атрибуты scale_factor и add_offset, поэтому добавим обработку этих двух атрибутов.
  7. Поставленная задача не требует обработки атрибута time_bnds, поэтому обработку этого атрибута не предусматриваем.

Дополнительные пояснения:

Во-первых, легко заметить, что программы CDO и NCO ко всем переменным нашего тестового файла добавили дополнительные атрибуты, в том числе, что важно, к переменной kq — атрибут scale_factor.

Поэтому обработку атрибутов scale_factor и add_offset рекомендуем добавлять всегда.

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

Вспомогательные процедуры

Для упрощения программы используем вспомогательные процедуры:

Подпрограмму 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 года.