metgrid/src/scan_input.F

   1 module scan_input
   2
   3     use netcdf
   4
   5     type input_handle_type
   6         integer :: ncid
   7         integer :: num_vars = 0
   8         integer :: current_var = 0
   9         integer, dimension(:), pointer :: varids => null()
  10         integer :: unlimited_dimid
  11     end type input_handle_type
  12
  13     type input_field_type
  14         character (len=64) :: name
  15         logical :: isTimeDependent = .false.
  16         integer :: varid = -1
  17         integer :: xtype = -1
  18         integer :: ndims = -1
  19         character (len=64), dimension(:), pointer :: dimnames
  20         integer, dimension(:), pointer :: dimlens
  21         integer, dimension(:), pointer :: dimids
  22         type (input_handle_type), pointer :: file_handle
  23
  24         !  Members to store field data
  25         real :: array0r
  26         real, dimension(:), pointer :: array1r => null()
  27         real, dimension(:,:), pointer :: array2r => null()
  28         real, dimension(:,:,:), pointer :: array3r => null()
  29         real, dimension(:,:,:,:), pointer :: array4r => null()
  30         double precision :: array0d
  31         double precision, dimension(:), pointer :: array1d => null()
  32         double precision, dimension(:,:), pointer :: array2d => null()
  33         double precision, dimension(:,:,:), pointer :: array3d => null()
  34         double precision, dimension(:,:,:,:), pointer :: array4d => null()
  35         integer :: array0i
  36         integer, dimension(:), pointer :: array1i => null()
  37         integer, dimension(:,:), pointer :: array2i => null()
  38         integer, dimension(:,:,:), pointer :: array3i => null()
  39     end type input_field_type
  40
  41     integer, parameter :: FIELD_TYPE_UNSUPPORTED = -1, &
  42                           FIELD_TYPE_REAL        =  1, &
  43                           FIELD_TYPE_DOUBLE      =  2, &
  44                           FIELD_TYPE_INTEGER     =  3, &
  45                           FIELD_TYPE_CHARACTER   =  4
  46
  47
  48     contains
  49
  50
  51     integer function scan_input_open(filename, handle, nRecords) result(stat)
  52
  53         implicit none
  54
  55         character (len=*), intent(in) :: filename
  56         type (input_handle_type), intent(out) :: handle
  57         integer, intent(out), optional :: nRecords
  58
  59         integer :: i
  60
  61         stat = 0
  62
  63         stat = nf90_open(trim(filename), NF90_NOWRITE, handle % ncid)
  64         if (stat /= NF90_NOERR) then
  65             stat = 1
  66             return
  67         end if
  68
  69         stat = nf90_inquire(handle % ncid, nVariables=handle % num_vars)
  70         if (stat /= NF90_NOERR) then
  71             stat = 1
  72             return
  73         end if
  74
  75         allocate(handle % varids(handle % num_vars))
  76
  77 #ifdef HAVE_NF90_INQ_VARIDS
  78         stat = nf90_inq_varids(handle % ncid, handle % num_vars, handle % varids)
  79         if (stat /= NF90_NOERR) then
  80             stat = 1
  81             return
  82         end if
  83 #else
  84         ! Newer versions of the netCDF4 library (perhaps newer than 4.2.0?)
  85         ! provide a function to return a list of all variable IDs in a file; if
  86         ! we are using an older version of the netCDF library, we can apparently
  87         ! assume that the variable IDs are numbered 1 through nVars.
  88         ! See http://www.unidata.ucar.edu/software/netcdf/docs/tutorial_ncids.html
  89         do i=1,handle % num_vars
  90             handle % varids(i) = i
  91         end do
  92 #endif
  93
  94         stat = nf90_inquire(handle % ncid, unlimitedDimId=handle % unlimited_dimid)
  95         if (stat /= NF90_NOERR) then
  96             stat = 1
  97             return
  98         end if
  99
 100         if (present(nRecords)) then
 101             stat = nf90_inquire_dimension(handle % ncid, handle % unlimited_dimid, len=nRecords)
 102             if (stat /= NF90_NOERR) then
 103                 stat = 1
 104                 return
 105             end if
 106
 107             ! In case we have an input file that no time-varying records but
 108             ! does have time-invariant fields, set nRecords = 1 so that we can
 109             ! at least extract these fields
 110             if ((nRecords == 0) .and. (handle % num_vars > 0)) then
 111                 nRecords = 1
 112             end if
 113         end if
 114
 115         handle % current_var = 1
 116
 117     end function scan_input_open
 118
 119
 120     integer function scan_input_close(handle) result(stat)
 121
 122         implicit none
 123
 124         type (input_handle_type), intent(inout) :: handle
 125
 126
 127         stat = 0
 128
 129         stat = nf90_close(handle % ncid)
 130         if (stat /= NF90_NOERR) then
 131             stat = 1
 132         end if
 133
 134         if (associated(handle % varids)) then
 135             deallocate(handle % varids)
 136         end if
 137         handle % current_var = 0
 138
 139     end function scan_input_close
 140
 141
 142     integer function scan_input_rewind(handle) result(stat)
 143
 144         implicit none
 145
 146         type (input_handle_type), intent(inout) :: handle
 147
 148
 149         stat = 0
 150
 151         handle % current_var = 1
 152
 153     end function scan_input_rewind
 154
 155
 156     integer function scan_input_next_field(handle, field) result(stat)
 157
 158         implicit none
 159
 160         type (input_handle_type), intent(inout), target :: handle
 161         type (input_field_type), intent(out) :: field
 162
 163         integer :: idim
 164
 165
 166         stat = 0
 167
 168         if (handle % current_var < 1 .or. handle % current_var > handle % num_vars) then
 169             stat = 1
 170             return
 171         end if
 172
 173         field % varid = handle % varids(handle % current_var)
 174         stat = nf90_inquire_variable(handle % ncid, field % varid, &
 175                                      name=field % name, &
 176                                      xtype=field % xtype, &
 177                                      ndims=field % ndims)
 178         if (stat /= NF90_NOERR) then
 179             stat = 1
 180             return
 181         end if
 182
 183         if (field % xtype == NF90_FLOAT) then
 184             field % xtype = FIELD_TYPE_REAL
 185         else if (field % xtype == NF90_DOUBLE) then
 186             field % xtype = FIELD_TYPE_DOUBLE
 187         else if (field % xtype == NF90_INT) then
 188             field % xtype = FIELD_TYPE_INTEGER
 189         else if (field % xtype == NF90_CHAR) then
 190             field % xtype = FIELD_TYPE_CHARACTER
 191         else
 192             field % xtype = FIELD_TYPE_UNSUPPORTED
 193         end if
 194
 195         allocate(field % dimids(field % ndims))
 196
 197         stat = nf90_inquire_variable(handle % ncid, field % varid, &
 198                                      dimids=field % dimids)
 199         if (stat /= NF90_NOERR) then
 200             stat = 1
 201             deallocate(field % dimids)
 202             return
 203         end if
 204
 205         allocate(field % dimlens(field % ndims))
 206         allocate(field % dimnames(field % ndims))
 207
 208         do idim=1,field % ndims
 209             stat = nf90_inquire_dimension(handle % ncid, field % dimids(idim), &
 210                                           name=field % dimnames(idim), &
 211                                           len=field % dimlens(idim))
 212             if (field % dimids(idim) == handle % unlimited_dimid) then
 213                 field % isTimeDependent = .true.
 214             end if
 215         end do
 216
 217         field % file_handle => handle
 218
 219         handle % current_var = handle % current_var + 1
 220
 221     end function scan_input_next_field
 222
 223
 224     integer function scan_input_for_field(handle, fieldname, field) result(stat)
 225
 226         implicit none
 227
 228         type (input_handle_type), intent(inout), target :: handle
 229         character (len=*), intent(in) :: fieldname
 230         type (input_field_type), intent(out) :: field
 231
 232         integer :: idim
 233
 234         stat = 0
 235
 236         stat = nf90_inq_varid(handle % ncid, trim(fieldname), field % varid)
 237         if (stat /= NF90_NOERR) then
 238             stat = 1
 239             return
 240         end if
 241
 242         stat = nf90_inquire_variable(handle % ncid, field % varid, &
 243                                      name=field % name, &
 244                                      xtype=field % xtype, &
 245                                      ndims=field % ndims)
 246         if (stat /= NF90_NOERR) then
 247             stat = 1
 248             return
 249         end if
 250
 251         if (field % xtype == NF90_FLOAT) then
 252             field % xtype = FIELD_TYPE_REAL
 253         else if (field % xtype == NF90_DOUBLE) then
 254             field % xtype = FIELD_TYPE_DOUBLE
 255         else if (field % xtype == NF90_INT) then
 256             field % xtype = FIELD_TYPE_INTEGER
 257         else if (field % xtype == NF90_CHAR) then
 258             field % xtype = FIELD_TYPE_CHARACTER
 259         else
 260             field % xtype = FIELD_TYPE_UNSUPPORTED
 261         end if
 262
 263         allocate(field % dimids(field % ndims))
 264
 265         stat = nf90_inquire_variable(handle % ncid, field % varid, &
 266                                      dimids=field % dimids)
 267         if (stat /= NF90_NOERR) then
 268             stat = 1
 269             deallocate(field % dimids)
 270             return
 271         end if
 272
 273         allocate(field % dimlens(field % ndims))
 274         allocate(field % dimnames(field % ndims))
 275
 276         do idim=1,field % ndims
 277             stat = nf90_inquire_dimension(handle % ncid, field % dimids(idim), &
 278                                           name=field % dimnames(idim), &
 279                                           len=field % dimlens(idim))
 280             if (field % dimids(idim) == handle % unlimited_dimid) then
 281                 field % isTimeDependent = .true.
 282             end if
 283         end do
 284
 285         field % file_handle => handle
 286
 287     end function scan_input_for_field
 288
 289
 290     integer function scan_input_read_field(field, frame) result(stat)
 291
 292         implicit none
 293
 294         type (input_field_type), intent(inout) :: field
 295         integer, intent(in), optional :: frame
 296
 297         integer :: local_frame
 298         integer, dimension(5) :: start, count
 299         real, dimension(1) :: temp1r
 300         double precision, dimension(1) :: temp1d
 301         integer, dimension(1) :: temp1i
 302
 303
 304         stat = 0
 305
 306         local_frame = 1
 307         if (present(frame)) then
 308             local_frame = frame
 309         end if
 310
 311         if (field % xtype == FIELD_TYPE_REAL) then
 312             if (field % ndims == 0 .or. (field % ndims == 1 .and. field % isTimeDependent)) then
 313                 if (field % isTimeDependent) then
 314                     start(1) = local_frame
 315                     count(1) = 1
 316                     stat = nf90_get_var(field % file_handle % ncid, field % varid, temp1r, &
 317                                         start=start(1:1), count=count(1:1))
 318                     field % array0r = temp1r(1)
 319                 else
 320                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array0r)
 321                 end if
 322             else if (field % ndims == 1 .or. (field % ndims == 2 .and. field % isTimeDependent)) then
 323                 if (field % isTimeDependent) then
 324                     start(1) = 1
 325                     count(1) = field % dimlens(1)
 326                     start(2) = local_frame
 327                     count(2) = 1
 328                     allocate(field % array1r(count(1)))
 329                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array1r, &
 330                                         start=start(1:2), count=count(1:2))
 331                 else
 332                     allocate(field % array1r(field%dimlens(1)))
 333                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array1r)
 334                 end if
 335             else if (field % ndims == 2 .or. (field % ndims == 3 .and. field % isTimeDependent)) then
 336                 if (field % isTimeDependent) then
 337                     start(1) = 1
 338                     count(1) = field % dimlens(1)
 339                     start(2) = 1
 340                     count(2) = field % dimlens(2)
 341                     start(3) = local_frame
 342                     count(3) = 1
 343                     allocate(field % array2r(count(1),count(2)))
 344                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array2r, &
 345                                         start=start(1:3), count=count(1:3))
 346                 else
 347                     allocate(field % array2r(field%dimlens(1),field%dimlens(2)))
 348                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array2r)
 349                 end if
 350             else if (field % ndims == 3 .or. (field % ndims == 4 .and. field % isTimeDependent)) then
 351                 if (field % isTimeDependent) then
 352                     start(1) = 1
 353                     count(1) = field % dimlens(1)
 354                     start(2) = 1
 355                     count(2) = field % dimlens(2)
 356                     start(3) = 1
 357                     count(3) = field % dimlens(3)
 358                     start(4) = local_frame
 359                     count(4) = 1
 360                     allocate(field % array3r(count(1),count(2),count(3)))
 361                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array3r, &
 362                                         start=start(1:4), count=count(1:4))
 363                 else
 364                     allocate(field % array3r(field%dimlens(1),field%dimlens(2),field%dimlens(3)))
 365                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array3r)
 366                 end if
 367             else if (field % ndims == 4 .or. (field % ndims == 5 .and. field % isTimeDependent)) then
 368                 if (field % isTimeDependent) then
 369                     start(1) = 1
 370                     count(1) = field % dimlens(1)
 371                     start(2) = 1
 372                     count(2) = field % dimlens(2)
 373                     start(3) = 1
 374                     count(3) = field % dimlens(3)
 375                     start(4) = 1
 376                     count(4) = field % dimlens(4)
 377                     start(5) = local_frame
 378                     count(5) = 1
 379                     allocate(field % array4r(count(1),count(2),count(3),count(4)))
 380                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array4r, &
 381                                         start=start(1:5), count=count(1:5))
 382                 else
 383                     allocate(field % array4r(field%dimlens(1),field%dimlens(2),field%dimlens(3),field%dimlens(4)))
 384                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array4r)
 385                 end if
 386             end if
 387         else if (field % xtype == FIELD_TYPE_DOUBLE) then
 388             if (field % ndims == 0 .or. (field % ndims == 1 .and. field % isTimeDependent)) then
 389                 if (field % isTimeDependent) then
 390                     start(1) = local_frame
 391                     count(1) = 1
 392                     stat = nf90_get_var(field % file_handle % ncid, field % varid, temp1d, &
 393                                         start=start(1:1), count=count(1:1))
 394                     field % array0d = temp1d(1)
 395                 else
 396                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array0d)
 397                 end if
 398             else if (field % ndims == 1 .or. (field % ndims == 2 .and. field % isTimeDependent)) then
 399                 if (field % isTimeDependent) then
 400                     start(1) = 1
 401                     count(1) = field % dimlens(1)
 402                     start(2) = local_frame
 403                     count(2) = 1
 404                     allocate(field % array1d(count(1)))
 405                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array1d, &
 406                                         start=start(1:2), count=count(1:2))
 407                 else
 408                     allocate(field % array1d(field%dimlens(1)))
 409                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array1d)
 410                 end if
 411             else if (field % ndims == 2 .or. (field % ndims == 3 .and. field % isTimeDependent)) then
 412                 if (field % isTimeDependent) then
 413                     start(1) = 1
 414                     count(1) = field % dimlens(1)
 415                     start(2) = 1
 416                     count(2) = field % dimlens(2)
 417                     start(3) = local_frame
 418                     count(3) = 1
 419                     allocate(field % array2d(count(1),count(2)))
 420                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array2d, &
 421                                         start=start(1:3), count=count(1:3))
 422                 else
 423                     allocate(field % array2d(field%dimlens(1),field%dimlens(2)))
 424                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array2d)
 425                 end if
 426             else if (field % ndims == 3 .or. (field % ndims == 4 .and. field % isTimeDependent)) then
 427                 if (field % isTimeDependent) then
 428                     start(1) = 1
 429                     count(1) = field % dimlens(1)
 430                     start(2) = 1
 431                     count(2) = field % dimlens(2)
 432                     start(3) = 1
 433                     count(3) = field % dimlens(3)
 434                     start(4) = local_frame
 435                     count(4) = 1
 436                     allocate(field % array3d(count(1),count(2),count(3)))
 437                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array3d, &
 438                                         start=start(1:4), count=count(1:4))
 439                 else
 440                     allocate(field % array3d(field%dimlens(1),field%dimlens(2),field%dimlens(3)))
 441                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array3d)
 442                 end if
 443             else if (field % ndims == 4 .or. (field % ndims == 5 .and. field % isTimeDependent)) then
 444                 if (field % isTimeDependent) then
 445                     start(1) = 1
 446                     count(1) = field % dimlens(1)
 447                     start(2) = 1
 448                     count(2) = field % dimlens(2)
 449                     start(3) = 1
 450                     count(3) = field % dimlens(3)
 451                     start(4) = 1
 452                     count(4) = field % dimlens(4)
 453                     start(5) = local_frame
 454                     count(5) = 1
 455                     allocate(field % array4d(count(1),count(2),count(3),count(4)))
 456                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array4d, &
 457                                         start=start(1:5), count=count(1:5))
 458                 else
 459                     allocate(field % array4d(field%dimlens(1),field%dimlens(2),field%dimlens(3),field%dimlens(4)))
 460                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array4d)
 461                 end if
 462             end if
 463
 464         else if (field % xtype == FIELD_TYPE_INTEGER) then
 465             if (field % ndims == 0 .or. (field % ndims == 1 .and. field % isTimeDependent)) then
 466                 if (field % isTimeDependent) then
 467                     start(1) = local_frame
 468                     count(1) = 1
 469                     stat = nf90_get_var(field % file_handle % ncid, field % varid, temp1i, &
 470                                         start=start(1:1), count=count(1:1))
 471                     field % array0i = temp1i(1)
 472                 else
 473                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array0i)
 474                 end if
 475             else if (field % ndims == 1 .or. (field % ndims == 2 .and. field % isTimeDependent)) then
 476                 if (field % isTimeDependent) then
 477                     start(1) = 1
 478                     count(1) = field % dimlens(1)
 479                     start(2) = local_frame
 480                     count(2) = 1
 481                     allocate(field % array1i(count(1)))
 482                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array1i, &
 483                                         start=start(1:2), count=count(1:2))
 484                 else
 485                     allocate(field % array1i(field%dimlens(1)))
 486                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array1i)
 487                 end if
 488             else if (field % ndims == 2 .or. (field % ndims == 3 .and. field % isTimeDependent)) then
 489                 if (field % isTimeDependent) then
 490                     start(1) = 1
 491                     count(1) = field % dimlens(1)
 492                     start(2) = 1
 493                     count(2) = field % dimlens(2)
 494                     start(3) = local_frame
 495                     count(3) = 1
 496                     allocate(field % array2i(count(1),count(2)))
 497                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array2i, &
 498                                         start=start(1:3), count=count(1:3))
 499                 else
 500                     allocate(field % array2i(field%dimlens(1),field%dimlens(2)))
 501                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array2i)
 502                 end if
 503             else if (field % ndims == 3 .or. (field % ndims == 4 .and. field % isTimeDependent)) then
 504                 if (field % isTimeDependent) then
 505                     start(1) = 1
 506                     count(1) = field % dimlens(1)
 507                     start(2) = 1
 508                     count(2) = field % dimlens(2)
 509                     start(3) = 1
 510                     count(3) = field % dimlens(3)
 511                     start(4) = local_frame
 512                     count(4) = 1
 513                     allocate(field % array3i(count(1),count(2),count(3)))
 514                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array3i, &
 515                                         start=start(1:4), count=count(1:4))
 516                 else
 517                     allocate(field % array3i(field%dimlens(1),field%dimlens(2),field%dimlens(3)))
 518                     stat = nf90_get_var(field % file_handle % ncid, field % varid, field % array3i)
 519                 end if
 520             end if
 521
 522         else if (field % xtype == FIELD_TYPE_CHARACTER) then
 523             write(0,*) ' '
 524             write(0,*) 'Processing of character fields is not supported; skipping read of field '//trim(field % name)
 525             write(0,*) ' '
 526
 527         else
 528             write(0,*) ' '
 529             write(0,*) 'Unsupported type; skipping read of field '//trim(field % name)
 530             write(0,*) ' '
 531         end if
 532
 533         if (stat /= NF90_NOERR) then
 534             write(0,*) ' '
 535             write(0,*) 'NetCDF error: reading '//trim(field % name)//' returned ', stat
 536             write(0,*) ' '
 537             stat = 1
 538         else
 539             stat = 0
 540         end if
 541
 542     end function scan_input_read_field
 543
 544
 545     integer function scan_input_free_field(field) result(stat)
 546
 547         implicit none
 548
 549         type (input_field_type), intent(inout) :: field
 550
 551
 552         stat = 0
 553
 554         if (associated(field % dimids)) then
 555             deallocate(field % dimids)
 556         end if
 557         if (associated(field % dimlens)) then
 558             deallocate(field % dimlens)
 559         end if
 560         if (associated(field % dimnames)) then
 561             deallocate(field % dimnames)
 562         end if
 563
 564         if (associated(field % array1r)) then
 565             deallocate(field % array1r)
 566         end if
 567         if (associated(field % array2r)) then
 568             deallocate(field % array2r)
 569         end if
 570         if (associated(field % array3r)) then
 571             deallocate(field % array3r)
 572         end if
 573         if (associated(field % array4r)) then
 574             deallocate(field % array4r)
 575         end if
 576
 577         if (associated(field % array1d)) then
 578             deallocate(field % array1d)
 579         end if
 580         if (associated(field % array2d)) then
 581             deallocate(field % array2d)
 582         end if
 583         if (associated(field % array3d)) then
 584             deallocate(field % array3d)
 585         end if
 586         if (associated(field % array4d)) then
 587             deallocate(field % array4d)
 588         end if
 589
 590         if (associated(field % array1i)) then
 591             deallocate(field % array1i)
 592         end if
 593         if (associated(field % array2i)) then
 594             deallocate(field % array2i)
 595         end if
 596         if (associated(field % array3i)) then
 597             deallocate(field % array3i)
 598         end if
 599
 600         nullify(field % file_handle)
 601
 602     end function scan_input_free_field
 603
 604 end module scan_input