metgrid/src/storage_module.F

   1 module storage_module
   2
   3    use datatype_module
   4    use minheap_module
   5    use misc_definitions_module
   6    use module_debug
   7    use parallel_module
   8
   9    ! Maximum umber of words to keep in memory at a time
  10    ! THIS MUST BE AT LEAST AS LARGE AS THE SIZE OF THE LARGEST ARRAY TO BE STORED
  11    integer, parameter :: MEMSIZE_MAX = 1E9
  12
  13    ! Name (when formatted as i9.9) of next file to be used as array storage
  14    integer :: next_filenumber = 1
  15
  16    ! Time counter used by policy for evicting arrays to Fortran units
  17    integer :: global_time = 0
  18
  19    ! Current memory usage of module
  20    integer :: memsize = 0
  21
  22    ! Primary head and tail pointers
  23    type (head_node), pointer :: head => null()
  24    type (head_node), pointer :: tail => null()
  25
  26    ! Pointer for get_next_output_fieldname
  27    type (head_node), pointer :: next_output_field  => null()
  28
  29    contains
  30
  31    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  32    ! Name: storage_init
  33    !
  34    ! Purpose: Initialize the storage module.
  35    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  36    subroutine storage_init()
  37
  38       implicit none
  39
  40       call init_heap()
  41
  42    end subroutine storage_init
  43
  44
  45    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  46    ! Name: reset_next_field
  47    !
  48    ! Purpose: Sets the next field to the first available field
  49    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  50    subroutine reset_next_field()
  51
  52       implicit none
  53
  54       next_output_field => head
  55
  56    end subroutine reset_next_field
  57
  58
  59    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  60    ! Name: storage_put_field
  61    !
  62    ! Purpose: Stores an fg_input type. Upon return, IT MUST NOT BE ASSUMED that
  63    !      store_me contains valid data, since all such data may have been written
  64    !      to a Fortran unit
  65    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  66    subroutine storage_put_field(store_me)
  67
  68       implicit none
  69
  70       ! Arguments
  71       type (fg_input), intent(in) :: store_me
  72
  73       ! Local variables
  74       integer :: funit
  75       logical :: is_used
  76       character (len=64) :: fname
  77       type (head_node), pointer :: name_cursor
  78       type (data_node), pointer :: data_cursor
  79       type (data_node), pointer :: newnode
  80       type (data_node), pointer :: evictnode
  81
  82       !CWH Initialize local pointer variables
  83       nullify(evictnode)       !MGD initialization for evictnode should not be necessary
  84
  85       ! We'll first see if there is already a list for this fieldname
  86       name_cursor => head
  87       do while (associated(name_cursor))
  88          if (primary_cmp(name_cursor%fg_data, store_me) == EQUAL) exit
  89          name_cursor => name_cursor%next
  90       end do
  91
  92       ! If not, create a new node in the primary list
  93       if (.not. associated(name_cursor)) then
  94          allocate(name_cursor)
  95          call dup(store_me, name_cursor%fg_data)
  96          nullify(name_cursor%fg_data%r_arr)
  97          nullify(name_cursor%fg_data%valid_mask)
  98          nullify(name_cursor%fg_data%modified_mask)
  99          nullify(name_cursor%fieldlist_head)
 100          nullify(name_cursor%fieldlist_tail)
 101          nullify(name_cursor%prev)
 102          name_cursor%next => head
 103          if (.not. associated(head)) tail => name_cursor
 104          head => name_cursor
 105       else
 106          if ((name_cursor%fg_data%header%time_dependent .and. .not. store_me%header%time_dependent) .or. &
 107              (.not. name_cursor%fg_data%header%time_dependent .and. store_me%header%time_dependent)) then
 108             call mprintf(.true.,ERROR,'Cannot combine time-independent data with '// &
 109                          'time-dependent data for field %s',s1=store_me%header%field)
 110          end if
 111       end if
 112
 113       ! At this point, name_cursor points to a valid head node for fieldname
 114       data_cursor => name_cursor%fieldlist_head
 115       do while ( associated(data_cursor) )
 116          if ((secondary_cmp(store_me, data_cursor%fg_data) == LESS) .or. &
 117              (secondary_cmp(store_me, data_cursor%fg_data) == EQUAL)) exit
 118          data_cursor => data_cursor%next
 119       end do
 120
 121       if (associated(data_cursor)) then
 122          if (secondary_cmp(store_me, data_cursor%fg_data) == EQUAL) then
 123             if (data_cursor%filenumber > 0) then
 124 ! BUG: Might need to deal with freeing up a file
 125 call mprintf(.true.,WARN,'WE NEED TO FREE THE FILE ASSOCIATED WITH DATA_CURSOR')
 126 call mprintf(.true.,WARN,'PLEASE REPORT THIS BUG TO THE DEVELOPER!')
 127             end if
 128             data_cursor%fg_data%r_arr => store_me%r_arr
 129             data_cursor%fg_data%valid_mask => store_me%valid_mask
 130             data_cursor%fg_data%modified_mask => store_me%modified_mask
 131             return
 132          end if
 133       end if
 134
 135       allocate(newnode)
 136       call dup(store_me, newnode%fg_data)
 137
 138       newnode%field_shape = shape(newnode%fg_data%r_arr)
 139       memsize = memsize + size(newnode%fg_data%r_arr)
 140       newnode%last_used = global_time
 141       global_time = global_time + 1
 142       newnode%filenumber = 0
 143       call add_to_heap(newnode)
 144
 145       do while (memsize > MEMSIZE_MAX)
 146          call get_min(evictnode)
 147          evictnode%filenumber = next_filenumber
 148          next_filenumber = next_filenumber + 1
 149          do funit=10,100
 150             inquire(unit=funit, opened=is_used)
 151             if (.not. is_used) exit
 152          end do
 153          memsize = memsize - size(evictnode%fg_data%r_arr)
 154          write(fname,'(i9.9,a2,i3.3)') evictnode%filenumber,'.p',my_proc_id
 155          open(funit,file=trim(fname),form='unformatted',status='unknown')
 156          write(funit) evictnode%fg_data%r_arr
 157          close(funit)
 158          deallocate(evictnode%fg_data%r_arr)
 159       end do
 160
 161       ! Inserting node at the tail of list
 162       if (.not. associated(data_cursor)) then
 163          newnode%prev => name_cursor%fieldlist_tail
 164          nullify(newnode%next)
 165
 166          ! List is actually empty
 167          if (.not. associated(name_cursor%fieldlist_head)) then
 168             name_cursor%fieldlist_head => newnode
 169             name_cursor%fieldlist_tail => newnode
 170          else
 171             name_cursor%fieldlist_tail%next => newnode
 172             name_cursor%fieldlist_tail => newnode
 173          end if
 174
 175       ! Inserting node at the head of list
 176       else if ((secondary_cmp(name_cursor%fieldlist_head%fg_data, newnode%fg_data) == GREATER) .or. &
 177                (secondary_cmp(name_cursor%fieldlist_head%fg_data, newnode%fg_data) == EQUAL)) then
 178          nullify(newnode%prev)
 179          newnode%next => name_cursor%fieldlist_head
 180          name_cursor%fieldlist_head%prev => newnode
 181          name_cursor%fieldlist_head => newnode
 182
 183       ! Inserting somewhere in the middle of the list
 184       else
 185          newnode%prev => data_cursor%prev
 186          newnode%next => data_cursor
 187          data_cursor%prev%next => newnode
 188          data_cursor%prev => newnode
 189       end if
 190
 191    end subroutine storage_put_field
 192
 193
 194    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 195    ! Name: storage_get_field
 196    !
 197    ! Purpose: Retrieves an fg_input type from storage; if the fg_input type whose
 198    !    header matches the header of get_me does not exist, istatus = 1 upon
 199    !    return; if the requested fg_input type is found, istatus = 0
 200    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 201    subroutine storage_get_field(get_me, istatus)
 202
 203       implicit none
 204
 205       ! Arguments
 206       type (fg_input), intent(inout) :: get_me
 207       integer, intent(out) :: istatus
 208
 209       ! Local variables
 210       integer :: funit
 211       logical :: is_used
 212       character (len=64) :: fname
 213       type (head_node), pointer :: name_cursor
 214       type (data_node), pointer :: data_cursor
 215       type (data_node), pointer :: evictnode
 216
 217       !CWH Initialize local pointer variables
 218       nullify(evictnode)     !MGD initialization for evictnodeshould not be necessary
 219
 220       global_time = global_time + 1
 221
 222       istatus = 1
 223
 224       ! We'll first see if there is already a list for this fieldname
 225       name_cursor => head
 226       do while (associated(name_cursor))
 227          if (primary_cmp(name_cursor%fg_data, get_me) == EQUAL) exit
 228          name_cursor => name_cursor%next
 229       end do
 230
 231       if (.not. associated(name_cursor)) return
 232
 233       ! At this point, name_cursor points to a valid head node for fieldname
 234       data_cursor => name_cursor%fieldlist_head
 235       do while ( associated(data_cursor) )
 236          if (secondary_cmp(get_me, data_cursor%fg_data) == EQUAL) then
 237             call dup(data_cursor%fg_data, get_me)
 238
 239             ! Before deciding whether we need to write an array to disk, first consider
 240             !   that reading the requested array will use memory
 241             if (data_cursor%filenumber > 0) then
 242                memsize = memsize + data_cursor%field_shape(1)*data_cursor%field_shape(2)
 243             end if
 244
 245             ! If we exceed our memory limit, we need to evict
 246             do while (memsize > MEMSIZE_MAX)
 247                call get_min(evictnode)
 248                evictnode%filenumber = next_filenumber
 249                next_filenumber = next_filenumber + 1
 250                do funit=10,100
 251                   inquire(unit=funit, opened=is_used)
 252                   if (.not. is_used) exit
 253                end do
 254                memsize = memsize - size(evictnode%fg_data%r_arr)
 255                write(fname,'(i9.9,a2,i3.3)') evictnode%filenumber,'.p',my_proc_id
 256                open(funit,file=trim(fname),form='unformatted',status='unknown')
 257                write(funit) evictnode%fg_data%r_arr
 258                close(funit)
 259                deallocate(evictnode%fg_data%r_arr)
 260             end do
 261
 262             ! Get requested array
 263             if (data_cursor%filenumber > 0) then
 264                data_cursor%last_used = global_time
 265                global_time = global_time + 1
 266                call add_to_heap(data_cursor)
 267                write(fname,'(i9.9,a2,i3.3)') data_cursor%filenumber,'.p',my_proc_id
 268                do funit=10,100
 269                   inquire(unit=funit, opened=is_used)
 270                   if (.not. is_used) exit
 271                end do
 272                open(funit,file=trim(fname),form='unformatted',status='old')
 273                allocate(data_cursor%fg_data%r_arr(data_cursor%field_shape(1),data_cursor%field_shape(2)))
 274                read(funit) data_cursor%fg_data%r_arr
 275                get_me%r_arr => data_cursor%fg_data%r_arr
 276                close(funit,status='delete')
 277                data_cursor%filenumber = 0
 278             else
 279                get_me%r_arr => data_cursor%fg_data%r_arr
 280
 281                call remove_index(data_cursor%heap_index)
 282                data_cursor%last_used = global_time
 283                global_time = global_time + 1
 284                call add_to_heap(data_cursor)
 285             end if
 286
 287             istatus = 0
 288             return
 289          end if
 290          data_cursor => data_cursor%next
 291       end do
 292
 293    end subroutine storage_get_field
 294
 295
 296    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 297    ! Name: storage_query_field
 298    !
 299    ! Purpose:
 300    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 301    subroutine storage_query_field(get_me, istatus)
 302
 303       implicit none
 304
 305       ! Arguments
 306       type (fg_input), intent(inout) :: get_me
 307       integer, intent(out) :: istatus
 308
 309       ! Local variables
 310       type (head_node), pointer :: name_cursor
 311       type (data_node), pointer :: data_cursor
 312
 313       istatus = 1
 314
 315       ! We'll first see if there is already a list for this fieldname
 316       name_cursor => head
 317       do while (associated(name_cursor))
 318          if (primary_cmp(name_cursor%fg_data, get_me) == EQUAL) exit
 319          name_cursor => name_cursor%next
 320       end do
 321
 322       if (.not. associated(name_cursor)) return
 323
 324       ! At this point, name_cursor points to a valid head node for fieldname
 325       data_cursor => name_cursor%fieldlist_head
 326       do while ( associated(data_cursor) )
 327          if (secondary_cmp(get_me, data_cursor%fg_data) == EQUAL) then
 328             get_me%r_arr => data_cursor%fg_data%r_arr
 329             get_me%valid_mask => data_cursor%fg_data%valid_mask
 330             get_me%modified_mask => data_cursor%fg_data%modified_mask
 331             istatus = 0
 332             return
 333          end if
 334          data_cursor => data_cursor%next
 335       end do
 336
 337    end subroutine storage_query_field
 338
 339
 340    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 341    ! Name: get_next_output_fieldname
 342    !
 343    ! Purpose:
 344    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 345    subroutine get_next_output_fieldname(nest_num, field_name, ndims, &
 346                                         min_level, max_level, &
 347                                         istagger, mem_order, dim_names, units, description, &
 348                                         sr_x, sr_y, derived_from, &
 349                                         istatus)
 350
 351       implicit none
 352
 353       ! Arguments
 354       integer, intent(in) :: nest_num
 355       integer, intent(out) :: ndims, min_level, max_level, istagger, istatus
 356       integer, intent(out) :: sr_x, sr_y
 357       character (len=128), intent(out) :: field_name, mem_order, units, description, derived_from
 358       character (len=128), dimension(3), intent(out) :: dim_names
 359
 360 #include "wrf_io_flags.h"
 361 #include "wrf_status_codes.h"
 362
 363       ! Local variables
 364       type (data_node), pointer :: data_cursor
 365
 366       istatus = 1
 367
 368       derived_from = ''
 369
 370       if (.not. associated(next_output_field)) return
 371
 372       min_level = 1
 373       max_level = 0
 374       ndims = 2
 375
 376       do while (max_level == 0 .and. associated(next_output_field))
 377
 378          data_cursor => next_output_field%fieldlist_head
 379          if (associated(data_cursor)) then
 380             if (.not. is_mask_field(data_cursor%fg_data)) then
 381                do while ( associated(data_cursor) )
 382                   istatus = 0
 383                   max_level = max_level + 1
 384                   data_cursor => data_cursor%next
 385                end do
 386             end if
 387          end if
 388
 389          if (max_level == 0) next_output_field => next_output_field%next
 390       end do
 391
 392       if (max_level > 0 .and. associated(next_output_field)) then
 393
 394          if (max_level > 1) ndims = 3
 395          if (ndims == 2) then
 396             mem_order = 'XY '
 397             dim_names(3) = ' '
 398          else
 399             mem_order = 'XYZ'
 400             if (is_time_dependent(next_output_field%fg_data)) then
 401                dim_names(3) = ' '
 402                dim_names(3)(1:32) = next_output_field%fg_data%header%vertical_coord
 403             else
 404                write(dim_names(3),'(a11,i4.4)') 'z-dimension', max_level
 405             end if
 406          end if
 407          field_name = get_fieldname(next_output_field%fg_data)
 408          istagger = get_staggering(next_output_field%fg_data)
 409          if (istagger == M .or. istagger == HH .or. istagger == VV) then
 410             dim_names(1) = 'west_east'
 411             dim_names(2) = 'south_north'
 412          else if (istagger == U) then
 413             dim_names(1) = 'west_east_stag'
 414             dim_names(2) = 'south_north'
 415          else if (istagger == V) then
 416             dim_names(1) = 'west_east'
 417             dim_names(2) = 'south_north_stag'
 418          else if (istagger == CORNER) then
 419             dim_names(1) = 'west_east_stag'
 420             dim_names(2) = 'south_north_stag'
 421          else
 422             dim_names(1) = 'i-dimension'
 423             dim_names(2) = 'j-dimension'
 424          end if
 425          units = get_units(next_output_field%fg_data)
 426          description = get_description(next_output_field%fg_data)
 427          call get_subgrid_dim_name(nest_num, field_name, dim_names(1:2), &
 428                                    sr_x, sr_y, istatus)
 429
 430          next_output_field => next_output_field%next
 431       end if
 432
 433    end subroutine get_next_output_fieldname
 434
 435
 436    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 437    ! Name: get_subgrid_dim_name
 438    !
 439    ! Purpose:
 440    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 441    subroutine get_subgrid_dim_name(nest_num, field_name, dimnames, &
 442                                    sub_x, sub_y, istatus)
 443
 444       use gridinfo_module
 445
 446       implicit none
 447
 448       ! Arguments
 449       integer, intent(in) :: nest_num
 450       integer, intent(out) :: sub_x, sub_y, istatus
 451       character(len=128), intent(in) :: field_name
 452       character(len=128), dimension(2), intent(inout) :: dimnames
 453
 454       ! Local variables
 455       integer :: idx, nlen
 456
 457       sub_x = next_output_field%fg_data%header%sr_x
 458       sub_y = next_output_field%fg_data%header%sr_y
 459
 460       if (sub_x > 1) then
 461         dimnames(1) = trim(dimnames(1))//"_subgrid"
 462       end if
 463       if (sub_y > 1) then
 464         dimnames(2) = trim(dimnames(2))//"_subgrid"
 465       end if
 466
 467       istatus = 0
 468
 469    end subroutine get_subgrid_dim_name
 470
 471
 472    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 473    ! Name: get_next_output_field
 474    !
 475    ! Purpose:
 476    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 477    subroutine get_next_output_field(field_name, r_array, &
 478                                     start_i, end_i, start_j, end_j, min_level, max_level, istatus)
 479
 480       implicit none
 481
 482       ! Arguments
 483       integer, intent(out) :: start_i, end_i, start_j, end_j, min_level, max_level, istatus
 484       real, pointer, dimension(:,:,:) :: r_array
 485       character (len=128), intent(out) :: field_name
 486
 487 #include "wrf_io_flags.h"
 488 #include "wrf_status_codes.h"
 489
 490       ! Local variables
 491       integer :: k
 492       type (data_node), pointer :: data_cursor
 493       type (fg_input) :: temp_field
 494
 495       istatus = 1
 496
 497       if (.not. associated(next_output_field)) return
 498
 499       min_level = 1
 500       max_level = 0
 501
 502       do while (max_level == 0 .and. associated(next_output_field))
 503
 504          data_cursor => next_output_field%fieldlist_head
 505          if (associated(data_cursor)) then
 506             if (.not. is_mask_field(data_cursor%fg_data)) then
 507                do while ( associated(data_cursor) )
 508                   istatus = 0
 509                   max_level = max_level + 1
 510                   data_cursor => data_cursor%next
 511                end do
 512             end if
 513          end if
 514
 515          if (max_level == 0) next_output_field => next_output_field%next
 516       end do
 517
 518       if (max_level > 0 .and. associated(next_output_field)) then
 519
 520          start_i = 1
 521          end_i = next_output_field%fieldlist_head%field_shape(1)
 522          start_j = 1
 523          end_j = next_output_field%fieldlist_head%field_shape(2)
 524
 525          allocate(r_array(next_output_field%fieldlist_head%field_shape(1), &
 526                           next_output_field%fieldlist_head%field_shape(2), &
 527                           max_level) )
 528
 529          k = 1
 530          data_cursor => next_output_field%fieldlist_head
 531          do while ( associated(data_cursor) )
 532             call dup(data_cursor%fg_data, temp_field)
 533             call storage_get_field(temp_field, istatus)
 534             r_array(:,:,k) = temp_field%r_arr
 535             k = k + 1
 536             data_cursor => data_cursor%next
 537          end do
 538
 539          field_name = get_fieldname(next_output_field%fg_data)
 540
 541          next_output_field => next_output_field%next
 542       end if
 543
 544    end subroutine get_next_output_field
 545
 546
 547    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 548    ! Name: storage_delete_field
 549    !
 550    ! Purpose: Deletes the stored fg_input type whose header matches delete_me
 551    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 552    subroutine storage_delete_field(delete_me)
 553
 554       implicit none
 555
 556       ! Arguments
 557       type (fg_input), intent(in) :: delete_me
 558
 559       ! Local variables
 560       integer :: funit
 561       logical :: is_used
 562       character (len=64) :: fname
 563       type (head_node), pointer :: name_cursor
 564       type (data_node), pointer :: data_cursor
 565
 566       ! We'll first see if there is a list for this fieldname
 567       name_cursor => head
 568       do while (associated(name_cursor))
 569          if (primary_cmp(name_cursor%fg_data, delete_me) == EQUAL) exit
 570          name_cursor => name_cursor%next
 571       end do
 572
 573       if (.not. associated(name_cursor)) return
 574
 575       ! At this point, name_cursor points to a valid head node for fieldname
 576       data_cursor => name_cursor%fieldlist_head
 577       do while ( associated(data_cursor) )
 578          if (secondary_cmp(delete_me, data_cursor%fg_data) == EQUAL) then
 579
 580             if (data_cursor%filenumber > 0) then
 581                do funit=10,100
 582                   inquire(unit=funit, opened=is_used)
 583                   if (.not. is_used) exit
 584                end do
 585                write(fname,'(i9.9,a2,i3.3)') data_cursor%filenumber,'.p',my_proc_id
 586                open(funit,file=trim(fname),form='unformatted',status='old')
 587                close(funit,status='delete')
 588             else
 589                call remove_index(data_cursor%heap_index)
 590                memsize = memsize - size(data_cursor%fg_data%r_arr)
 591                deallocate(data_cursor%fg_data%r_arr)
 592             end if
 593             if (associated(data_cursor%fg_data%valid_mask)) call bitarray_destroy(data_cursor%fg_data%valid_mask)
 594             nullify(data_cursor%fg_data%valid_mask)
 595             if (associated(data_cursor%fg_data%modified_mask)) call bitarray_destroy(data_cursor%fg_data%modified_mask)
 596             nullify(data_cursor%fg_data%modified_mask)
 597
 598             ! Only item in the list
 599             if (.not. associated(data_cursor%next) .and. &
 600                 .not. associated(data_cursor%prev)) then
 601                nullify(name_cursor%fieldlist_head)
 602                nullify(name_cursor%fieldlist_tail)
 603                deallocate(data_cursor)
 604 ! DO WE REMOVE THIS HEADER NODE AT THIS POINT?
 605                return
 606
 607             ! Head of the list
 608             else if (.not. associated(data_cursor%prev)) then
 609                name_cursor%fieldlist_head => data_cursor%next
 610                nullify(data_cursor%next%prev)
 611                deallocate(data_cursor)
 612                return
 613
 614             ! Tail of the list
 615             else if (.not. associated(data_cursor%next)) then
 616                name_cursor%fieldlist_tail => data_cursor%prev
 617                nullify(data_cursor%prev%next)
 618                deallocate(data_cursor)
 619                return
 620
 621             ! Middle of the list
 622             else
 623                data_cursor%prev%next => data_cursor%next
 624                data_cursor%next%prev => data_cursor%prev
 625                deallocate(data_cursor)
 626                return
 627
 628             end if
 629
 630          end if
 631          data_cursor => data_cursor%next
 632       end do
 633
 634    end subroutine storage_delete_field
 635
 636
 637    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 638    ! Name: storage_delete_all_td
 639    !
 640    ! Purpose: Deletes the stored time-dependent data
 641    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 642    subroutine storage_delete_all_td()
 643
 644       implicit none
 645
 646       ! Local variables
 647       integer :: funit
 648       logical :: is_used
 649       character (len=64) :: fname
 650       type (head_node), pointer :: name_cursor
 651       type (data_node), pointer :: data_cursor, next_cursor
 652
 653       ! We'll first see if there is a list for this fieldname
 654       name_cursor => head
 655       do while (associated(name_cursor))
 656
 657          data_cursor => name_cursor%fieldlist_head
 658          do while ( associated(data_cursor) )
 659             if ( is_time_dependent(data_cursor%fg_data) ) then
 660
 661                if (data_cursor%filenumber > 0) then
 662                   do funit=10,100
 663                      inquire(unit=funit, opened=is_used)
 664                      if (.not. is_used) exit
 665                   end do
 666                   write(fname,'(i9.9,a2,i3.3)') data_cursor%filenumber,'.p',my_proc_id
 667                   open(funit,file=trim(fname),form='unformatted',status='old')
 668                   close(funit,status='delete')
 669                else
 670                   call remove_index(data_cursor%heap_index)
 671                   memsize = memsize - size(data_cursor%fg_data%r_arr)
 672                   deallocate(data_cursor%fg_data%r_arr)
 673                end if
 674                if (associated(data_cursor%fg_data%valid_mask)) call bitarray_destroy(data_cursor%fg_data%valid_mask)
 675                nullify(data_cursor%fg_data%valid_mask)
 676                if (associated(data_cursor%fg_data%modified_mask)) call bitarray_destroy(data_cursor%fg_data%modified_mask)
 677                nullify(data_cursor%fg_data%modified_mask)
 678
 679                ! We should handle individual cases, that way we can deal with a list
 680                !   that has both time independent and time dependent nodes in it.
 681
 682                ! Only item in the list
 683                if (.not. associated(data_cursor%next) .and. &
 684                    .not. associated(data_cursor%prev)) then
 685                   next_cursor => null()
 686                   nullify(name_cursor%fieldlist_head)
 687                   nullify(name_cursor%fieldlist_tail)
 688                   deallocate(data_cursor)
 689 ! DO WE REMOVE THIS HEADER NODE AT THIS POINT?
 690
 691                ! Head of the list
 692                else if (.not. associated(data_cursor%prev)) then
 693                   name_cursor%fieldlist_head => data_cursor%next
 694                   next_cursor => data_cursor%next
 695                   nullify(data_cursor%next%prev)
 696                   deallocate(data_cursor)
 697
 698                ! Tail of the list
 699                else if (.not. associated(data_cursor%next)) then
 700 ! THIS CASE SHOULD PROBABLY NOT OCCUR
 701                   name_cursor%fieldlist_tail => data_cursor%prev
 702                   next_cursor => null()
 703                   nullify(data_cursor%prev%next)
 704                   deallocate(data_cursor)
 705
 706                ! Middle of the list
 707                else
 708 ! THIS CASE SHOULD PROBABLY NOT OCCUR
 709                   next_cursor => data_cursor%next
 710                   data_cursor%prev%next => data_cursor%next
 711                   data_cursor%next%prev => data_cursor%prev
 712                   deallocate(data_cursor)
 713
 714                end if
 715
 716             end if
 717             data_cursor => next_cursor
 718          end do
 719
 720          name_cursor => name_cursor%next
 721       end do
 722
 723    end subroutine storage_delete_all_td
 724
 725
 726    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 727    ! Name: storage_get_levels
 728    !
 729    ! Purpose: Returns a list of all levels for the field indicated in the_header.
 730    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 731    subroutine storage_get_levels(the_header, list)
 732
 733       implicit none
 734
 735       ! Arguments
 736       integer, pointer, dimension(:) :: list
 737       type (fg_input), intent(in) :: the_header
 738
 739       ! Local variables
 740       integer :: n
 741       type (head_node), pointer :: name_cursor
 742       type (data_node), pointer :: data_cursor
 743
 744       if (associated(list)) deallocate(list)
 745       nullify(list)
 746
 747       ! We'll first see if there is a list for this header
 748       name_cursor => head
 749       do while (associated(name_cursor))
 750          if (primary_cmp(name_cursor%fg_data, the_header) == EQUAL) exit
 751          name_cursor => name_cursor%next
 752       end do
 753
 754       if (.not. associated(name_cursor)) return
 755
 756       n = 0
 757       ! At this point, name_cursor points to a valid head node for fieldname
 758       data_cursor => name_cursor%fieldlist_head
 759       do while ( associated(data_cursor) )
 760          n = n + 1
 761          if (.not. associated(data_cursor%next)) exit
 762          data_cursor => data_cursor%next
 763       end do
 764
 765       if (n > 0) allocate(list(n))
 766
 767       n = 1
 768       do while ( associated(data_cursor) )
 769          list(n) = get_level(data_cursor%fg_data)
 770          n = n + 1
 771          data_cursor => data_cursor%prev
 772       end do
 773
 774    end subroutine storage_get_levels
 775
 776
 777    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 778    ! Name: storage_delete_all
 779    !
 780    ! Purpose: Deletes all data, both time-independent and time-dependent.
 781    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 782    subroutine storage_delete_all()
 783
 784       implicit none
 785
 786       ! Local variables
 787       integer :: funit
 788       logical :: is_used
 789       character (len=64) :: fname
 790       type (head_node), pointer :: name_cursor
 791       type (data_node), pointer :: data_cursor
 792
 793       ! We'll first see if there is already a list for this fieldname
 794       name_cursor => head
 795       do while (associated(name_cursor))
 796
 797          if (associated(name_cursor%fieldlist_head)) then
 798             data_cursor => name_cursor%fieldlist_head
 799             do while ( associated(data_cursor) )
 800                name_cursor%fieldlist_head => data_cursor%next
 801
 802                if (data_cursor%filenumber > 0) then
 803                   do funit=10,100
 804                      inquire(unit=funit, opened=is_used)
 805                      if (.not. is_used) exit
 806                   end do
 807                   write(fname,'(i9.9,a2,i3.3)') data_cursor%filenumber,'.p',my_proc_id
 808                   open(funit,file=trim(fname),form='unformatted',status='old')
 809                   close(funit,status='delete')
 810                else
 811                   call remove_index(data_cursor%heap_index)
 812                   memsize = memsize - size(data_cursor%fg_data%r_arr)
 813                   deallocate(data_cursor%fg_data%r_arr)
 814                end if
 815                if (associated(data_cursor%fg_data%valid_mask)) call bitarray_destroy(data_cursor%fg_data%valid_mask)
 816                nullify(data_cursor%fg_data%valid_mask)
 817                if (associated(data_cursor%fg_data%modified_mask)) call bitarray_destroy(data_cursor%fg_data%modified_mask)
 818                nullify(data_cursor%fg_data%modified_mask)
 819
 820                deallocate(data_cursor)
 821                data_cursor => name_cursor%fieldlist_head
 822             end do
 823          end if
 824
 825          head => name_cursor%next
 826          deallocate(name_cursor)
 827          name_cursor => head
 828       end do
 829
 830       nullify(head)
 831       nullify(tail)
 832
 833       call heap_destroy()
 834
 835    end subroutine storage_delete_all
 836
 837
 838    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 839    ! Name: storage_get_all_headers
 840    !
 841    ! Purpose:
 842    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 843    subroutine storage_get_all_headers(header_list)
 844
 845       implicit none
 846
 847       ! Arguments
 848       type (fg_input), pointer, dimension(:) :: header_list
 849
 850       ! Local variables
 851       integer :: nheaders
 852       type (head_node), pointer :: name_cursor
 853
 854       nullify(header_list)
 855
 856       ! First find out how many time-dependent headers there are
 857       name_cursor => head
 858       nheaders = 0
 859       do while (associated(name_cursor))
 860          if (associated(name_cursor%fieldlist_head)) then
 861             if (.not. is_mask_field(name_cursor%fieldlist_head%fg_data)) then
 862                nheaders = nheaders + 1
 863             end if
 864          end if
 865          name_cursor => name_cursor%next
 866       end do
 867
 868       allocate(header_list(nheaders))
 869
 870       name_cursor => head
 871       nheaders = 0
 872       do while (associated(name_cursor))
 873          if (associated(name_cursor%fieldlist_head)) then
 874             if (.not. is_mask_field(name_cursor%fieldlist_head%fg_data)) then
 875                nheaders = nheaders + 1
 876                call dup(name_cursor%fieldlist_head%fg_data, header_list(nheaders))
 877             end if
 878          end if
 879          name_cursor => name_cursor%next
 880       end do
 881
 882    end subroutine storage_get_all_headers
 883
 884
 885    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 886    ! Name: storage_get_all_td_headers
 887    !
 888    ! Purpose:
 889    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 890    subroutine storage_get_td_headers(header_list)
 891
 892       implicit none
 893
 894       ! Arguments
 895       type (fg_input), pointer, dimension(:) :: header_list
 896
 897       ! Local variables
 898       integer :: nheaders
 899       type (head_node), pointer :: name_cursor
 900
 901       nullify(header_list)
 902
 903       ! First find out how many time-dependent headers there are
 904       name_cursor => head
 905       nheaders = 0
 906       do while (associated(name_cursor))
 907          if (associated(name_cursor%fieldlist_head)) then
 908             if (is_time_dependent(name_cursor%fieldlist_head%fg_data) .and. &
 909                 .not. is_mask_field(name_cursor%fieldlist_head%fg_data)) then
 910                nheaders = nheaders + 1
 911             end if
 912          end if
 913          name_cursor => name_cursor%next
 914       end do
 915
 916       allocate(header_list(nheaders))
 917
 918       name_cursor => head
 919       nheaders = 0
 920       do while (associated(name_cursor))
 921          if (associated(name_cursor%fieldlist_head)) then
 922             if (is_time_dependent(name_cursor%fieldlist_head%fg_data) .and. &
 923                 .not. is_mask_field(name_cursor%fieldlist_head%fg_data)) then
 924                nheaders = nheaders + 1
 925                call dup(name_cursor%fieldlist_head%fg_data, header_list(nheaders))
 926             end if
 927          end if
 928          name_cursor => name_cursor%next
 929       end do
 930
 931    end subroutine storage_get_td_headers
 932
 933
 934    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 935    ! Name: storage_print_fields
 936    !
 937    ! Purpose:
 938    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 939    subroutine storage_print_fields()
 940
 941       use list_module
 942       use stringutil
 943
 944       implicit none
 945
 946       ! Local variables
 947       integer :: i, j, k, lmax, n_fields, n_levels, max_levels, itemp
 948       logical, allocatable, dimension(:,:) :: field_has_level
 949       integer, allocatable, dimension(:) :: all_levels
 950       integer, pointer, dimension(:) :: ilevels
 951       character (len=128), allocatable, dimension(:) :: fieldname_list
 952       character (len=9) :: ctemp
 953       type (fg_input), pointer, dimension(:) :: header_list
 954
 955       type (list) :: all_levs
 956
 957       !CWH Initialize local pointer variables
 958       nullify(ilevels)
 959       nullify(header_list)     !MGD initialization for header_list should not be necessary
 960
 961       call list_init(all_levs)
 962       call storage_get_td_headers(header_list)
 963       n_fields = size(header_list)
 964
 965       allocate(fieldname_list(n_fields))
 966
 967       max_levels = 0
 968
 969       do i=1,n_fields
 970          fieldname_list(i) = header_list(i)%header%field
 971          call storage_get_levels(header_list(i), ilevels)
 972          do j=1,size(ilevels)
 973             if (.not. list_search(all_levs, ikey=ilevels(j), ivalue=itemp)) then
 974                call list_insert(all_levs, ikey=ilevels(j), ivalue=ilevels(j))
 975             end if
 976          end do
 977          n_levels = size(ilevels)
 978          if (n_levels > max_levels) max_levels = n_levels
 979          if (associated(ilevels)) deallocate(ilevels)
 980       end do
 981
 982       max_levels = list_length(all_levs)
 983
 984       allocate(all_levels(max_levels))
 985       allocate(field_has_level(n_fields,max_levels))
 986
 987       field_has_level(:,:) = .false.
 988
 989       lmax = 0
 990       do i=1,n_fields
 991          call storage_get_levels(header_list(i), ilevels)
 992          n_levels = size(ilevels)
 993          do j=1,n_levels
 994             do k=1,lmax
 995                if (all_levels(k) == ilevels(j)) exit
 996             end do
 997             if (k > lmax) then
 998                all_levels(k) = ilevels(j)
 999                lmax = lmax + 1
1000             end if
1001             field_has_level(i,k) = .true.
1002          end do
1003          if (associated(ilevels)) deallocate(ilevels)
1004       end do
1005
1006       call mprintf(.true.,DEBUG,'        .',newline=.false.)
1007       do i=1,n_fields
1008          write(ctemp,'(a9)') fieldname_list(i)(1:9)
1009          call right_justify(ctemp,9)
1010          call mprintf(.true.,DEBUG,ctemp,newline=.false.)
1011       end do
1012       call mprintf(.true.,DEBUG,' ',newline=.true.)
1013       do j=1,max_levels
1014          write(ctemp,'(i9)') all_levels(j)
1015          call mprintf(.true.,DEBUG,'%s ',s1=ctemp,newline=.false.)
1016          do i=1,n_fields
1017             if (field_has_level(i,j)) then
1018                call mprintf(.true.,DEBUG,'        X',newline=.false.)
1019             else
1020                call mprintf(.true.,DEBUG,'        -',newline=.false.)
1021             end if
1022          end do
1023          call mprintf(.true.,DEBUG,' ',newline=.true.)
1024       end do
1025
1026       deallocate(all_levels)
1027       deallocate(field_has_level)
1028       deallocate(fieldname_list)
1029       deallocate(header_list)
1030
1031       call list_destroy(all_levs)
1032
1033    end subroutine storage_print_fields
1034
1035
1036    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1037    ! Name: find_missing_values
1038    !
1039    ! Purpose:
1040    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1041    subroutine find_missing_values()
1042
1043       implicit none
1044
1045       ! Local variables
1046       integer :: i, j
1047       logical :: found_missing
1048       type (head_node), pointer :: name_cursor
1049       type (data_node), pointer :: data_cursor
1050
1051       found_missing = .false.
1052
1053       name_cursor => head
1054       do while (associated(name_cursor))
1055
1056          if (associated(name_cursor%fieldlist_head)) then
1057             data_cursor => name_cursor%fieldlist_head
1058             do while ( associated(data_cursor) )
1059                if (.not. associated(data_cursor%fg_data%valid_mask)) then
1060                   call mprintf(.true.,INFORM, &
1061                                'Field %s does not have a valid mask and will not be checked for missing values', &
1062                                s1=data_cursor%fg_data%header%field)
1063                else
1064                   ILOOP: do i=1,data_cursor%fg_data%header%dim1(2)-data_cursor%fg_data%header%dim1(1)+1
1065                   JLOOP: do j=1,data_cursor%fg_data%header%dim2(2)-data_cursor%fg_data%header%dim2(1)+1
1066                      if (.not. bitarray_test(data_cursor%fg_data%valid_mask,i,j)) then
1067                         found_missing = .true.
1068                         call mprintf(.true.,WARN,'Field %s has missing values at level %i at (i,j)=(%i,%i)', &
1069                                      s1=data_cursor%fg_data%header%field, &
1070                                      i1=data_cursor%fg_data%header%vertical_level, &
1071                                      i2=i+data_cursor%fg_data%header%dim1(1)-1, &
1072                                      i3=j+data_cursor%fg_data%header%dim2(1)-1)
1073                         exit ILOOP
1074                      end if
1075                   end do JLOOP
1076                   end do ILOOP
1077                end if
1078                data_cursor => data_cursor%next
1079             end do
1080          end if
1081
1082          name_cursor => name_cursor%next
1083       end do
1084
1085       call mprintf(found_missing,ERROR,'Missing values encountered in interpolated fields. Stopping.')
1086
1087    end subroutine find_missing_values
1088
1089
1090    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1091    ! Name: storage_print_headers
1092    !
1093    ! Purpose:
1094    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1095    subroutine storage_print_headers()
1096
1097       implicit none
1098
1099       ! Local variables
1100       type (head_node), pointer :: name_cursor
1101       type (data_node), pointer :: data_cursor
1102
1103       call mprintf(.true.,DEBUG,'>>>> STORED FIELDS <<<<')
1104       call mprintf(.true.,DEBUG,'=======================')
1105
1106       ! We'll first see if there is already a list for this fieldname
1107       name_cursor => head
1108       do while (associated(name_cursor))
1109          call print_header(name_cursor%fg_data)
1110
1111          if (associated(name_cursor%fieldlist_head)) then
1112             data_cursor => name_cursor%fieldlist_head
1113             do while ( associated(data_cursor) )
1114                call mprintf(.true.,DEBUG,'  - %i', i1=get_level(data_cursor%fg_data))
1115                call mprintf(.true.,DEBUG,' ')
1116                data_cursor => data_cursor%next
1117             end do
1118          end if
1119
1120          name_cursor => name_cursor%next
1121       end do
1122
1123    end subroutine storage_print_headers
1124
1125 end module storage_module