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, is_subgrid_var, 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       logical, intent(in) :: is_subgrid_var
 358       character (len=128), intent(out) :: field_name, mem_order, units, description, derived_from
 359       character (len=128), dimension(3), intent(out) :: dim_names
 360
 361 #include "wrf_io_flags.h"
 362 #include "wrf_status_codes.h"
 363
 364       ! Local variables
 365       type (data_node), pointer :: data_cursor
 366
 367       istatus = 1
 368
 369       derived_from = ''
 370
 371       if (.not. associated(next_output_field)) return
 372
 373       min_level = 1
 374       max_level = 0
 375       ndims = 2
 376
 377       do while (max_level == 0 .and. associated(next_output_field))
 378
 379          data_cursor => next_output_field%fieldlist_head
 380          if (associated(data_cursor)) then
 381             if (.not. is_mask_field(data_cursor%fg_data)) then
 382                do while ( associated(data_cursor) )
 383                   istatus = 0
 384                   max_level = max_level + 1
 385                   data_cursor => data_cursor%next
 386                end do
 387             end if
 388          end if
 389
 390          if (max_level == 0) next_output_field => next_output_field%next
 391       end do
 392
 393       if (max_level > 0 .and. associated(next_output_field)) then
 394
 395          if (max_level > 1) ndims = 3
 396          if (ndims == 2) then
 397             mem_order = 'XY '
 398             dim_names(3) = ' '
 399          else
 400             mem_order = 'XYZ'
 401             if (is_time_dependent(next_output_field%fg_data)) then
 402                dim_names(3) = ' '
 403                dim_names(3)(1:32) = next_output_field%fg_data%header%vertical_coord
 404             else
 405                write(dim_names(3),'(a11,i4.4)') 'z-dimension', max_level
 406             end if
 407          end if
 408          field_name = get_fieldname(next_output_field%fg_data)
 409          istagger = get_staggering(next_output_field%fg_data)
 410          if (istagger == M .or. istagger == HH .or. istagger == VV) then
 411             dim_names(1) = 'west_east'
 412             dim_names(2) = 'south_north'
 413          else if (istagger == U) then
 414             dim_names(1) = 'west_east_stag'
 415             dim_names(2) = 'south_north'
 416          else if (istagger == V) then
 417             dim_names(1) = 'west_east'
 418             dim_names(2) = 'south_north_stag'
 419          else if (istagger == CORNER) then
 420             dim_names(1) = 'west_east_stag'
 421             dim_names(2) = 'south_north_stag'
 422          else
 423             dim_names(1) = 'i-dimension'
 424             dim_names(2) = 'j-dimension'
 425          end if
 426          units = get_units(next_output_field%fg_data)
 427          description = get_description(next_output_field%fg_data)
 428          call get_subgrid_dim_name(nest_num, field_name, dim_names(1:2), &
 429                                    sr_x, sr_y, istatus)
 430
 431          next_output_field => next_output_field%next
 432       end if
 433
 434    end subroutine get_next_output_fieldname
 435
 436
 437    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 438    ! Name: get_subgrid_dim_name
 439    !
 440    ! Purpose:
 441    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 442    subroutine get_subgrid_dim_name(nest_num, field_name, dimnames, &
 443                                    sub_x, sub_y, istatus)
 444
 445       use gridinfo_module
 446
 447       implicit none
 448
 449       ! Arguments
 450       integer, intent(in) :: nest_num
 451       integer, intent(out) :: sub_x, sub_y, istatus
 452       character(len=128), intent(in) :: field_name
 453       character(len=128), dimension(2), intent(inout) :: dimnames
 454
 455       ! Local variables
 456       integer :: idx, nlen
 457
 458       sub_x = next_output_field%fg_data%header%sr_x
 459       sub_y = next_output_field%fg_data%header%sr_y
 460
 461       if (sub_x > 1) then
 462         dimnames(1) = trim(dimnames(1))//"_subgrid"
 463       end if
 464       if (sub_y > 1) then
 465         dimnames(2) = trim(dimnames(2))//"_subgrid"
 466       end if
 467
 468       istatus = 0
 469
 470    end subroutine get_subgrid_dim_name
 471
 472
 473    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 474    ! Name: get_next_output_field
 475    !
 476    ! Purpose:
 477    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 478    subroutine get_next_output_field(field_name, r_array, &
 479                                     start_i, end_i, start_j, end_j, min_level, max_level, istatus)
 480
 481       implicit none
 482
 483       ! Arguments
 484       integer, intent(out) :: start_i, end_i, start_j, end_j, min_level, max_level, istatus
 485       real, pointer, dimension(:,:,:) :: r_array
 486       character (len=128), intent(out) :: field_name
 487
 488 #include "wrf_io_flags.h"
 489 #include "wrf_status_codes.h"
 490
 491       ! Local variables
 492       integer :: k
 493       type (data_node), pointer :: data_cursor
 494       type (fg_input) :: temp_field
 495
 496       istatus = 1
 497
 498       if (.not. associated(next_output_field)) return
 499
 500       min_level = 1
 501       max_level = 0
 502
 503       do while (max_level == 0 .and. associated(next_output_field))
 504
 505          data_cursor => next_output_field%fieldlist_head
 506          if (associated(data_cursor)) then
 507             if (.not. is_mask_field(data_cursor%fg_data)) then
 508                do while ( associated(data_cursor) )
 509                   istatus = 0
 510                   max_level = max_level + 1
 511                   data_cursor => data_cursor%next
 512                end do
 513             end if
 514          end if
 515
 516          if (max_level == 0) next_output_field => next_output_field%next
 517       end do
 518
 519       if (max_level > 0 .and. associated(next_output_field)) then
 520
 521          start_i = 1
 522          end_i = next_output_field%fieldlist_head%field_shape(1)
 523          start_j = 1
 524          end_j = next_output_field%fieldlist_head%field_shape(2)
 525
 526          allocate(r_array(next_output_field%fieldlist_head%field_shape(1), &
 527                           next_output_field%fieldlist_head%field_shape(2), &
 528                           max_level) )
 529
 530          k = 1
 531          data_cursor => next_output_field%fieldlist_head
 532          do while ( associated(data_cursor) )
 533             call dup(data_cursor%fg_data, temp_field)
 534             call storage_get_field(temp_field, istatus)
 535             r_array(:,:,k) = temp_field%r_arr
 536             k = k + 1
 537             data_cursor => data_cursor%next
 538          end do
 539
 540          field_name = get_fieldname(next_output_field%fg_data)
 541
 542          next_output_field => next_output_field%next
 543       end if
 544
 545    end subroutine get_next_output_field
 546
 547
 548    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 549    ! Name: storage_delete_field
 550    !
 551    ! Purpose: Deletes the stored fg_input type whose header matches delete_me
 552    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 553    subroutine storage_delete_field(delete_me)
 554
 555       implicit none
 556
 557       ! Arguments
 558       type (fg_input), intent(in) :: delete_me
 559
 560       ! Local variables
 561       integer :: funit
 562       logical :: is_used
 563       character (len=64) :: fname
 564       type (head_node), pointer :: name_cursor
 565       type (data_node), pointer :: data_cursor
 566
 567       ! We'll first see if there is a list for this fieldname
 568       name_cursor => head
 569       do while (associated(name_cursor))
 570          if (primary_cmp(name_cursor%fg_data, delete_me) == EQUAL) exit
 571          name_cursor => name_cursor%next
 572       end do
 573
 574       if (.not. associated(name_cursor)) return
 575
 576       ! At this point, name_cursor points to a valid head node for fieldname
 577       data_cursor => name_cursor%fieldlist_head
 578       do while ( associated(data_cursor) )
 579          if (secondary_cmp(delete_me, data_cursor%fg_data) == EQUAL) then
 580
 581             if (data_cursor%filenumber > 0) then
 582                do funit=10,100
 583                   inquire(unit=funit, opened=is_used)
 584                   if (.not. is_used) exit
 585                end do
 586                write(fname,'(i9.9,a2,i3.3)') data_cursor%filenumber,'.p',my_proc_id
 587                open(funit,file=trim(fname),form='unformatted',status='old')
 588                close(funit,status='delete')
 589             else
 590                call remove_index(data_cursor%heap_index)
 591                memsize = memsize - size(data_cursor%fg_data%r_arr)
 592                deallocate(data_cursor%fg_data%r_arr)
 593             end if
 594             if (associated(data_cursor%fg_data%valid_mask)) call bitarray_destroy(data_cursor%fg_data%valid_mask)
 595             nullify(data_cursor%fg_data%valid_mask)
 596             if (associated(data_cursor%fg_data%modified_mask)) call bitarray_destroy(data_cursor%fg_data%modified_mask)
 597             nullify(data_cursor%fg_data%modified_mask)
 598
 599             ! Only item in the list
 600             if (.not. associated(data_cursor%next) .and. &
 601                 .not. associated(data_cursor%prev)) then
 602                nullify(name_cursor%fieldlist_head)
 603                nullify(name_cursor%fieldlist_tail)
 604                deallocate(data_cursor)
 605 ! DO WE REMOVE THIS HEADER NODE AT THIS POINT?
 606                return
 607
 608             ! Head of the list
 609             else if (.not. associated(data_cursor%prev)) then
 610                name_cursor%fieldlist_head => data_cursor%next
 611                nullify(data_cursor%next%prev)
 612                deallocate(data_cursor)
 613                return
 614
 615             ! Tail of the list
 616             else if (.not. associated(data_cursor%next)) then
 617                name_cursor%fieldlist_tail => data_cursor%prev
 618                nullify(data_cursor%prev%next)
 619                deallocate(data_cursor)
 620                return
 621
 622             ! Middle of the list
 623             else
 624                data_cursor%prev%next => data_cursor%next
 625                data_cursor%next%prev => data_cursor%prev
 626                deallocate(data_cursor)
 627                return
 628
 629             end if
 630
 631          end if
 632          data_cursor => data_cursor%next
 633       end do
 634
 635    end subroutine storage_delete_field
 636
 637
 638    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 639    ! Name: storage_delete_all_td
 640    !
 641    ! Purpose: Deletes the stored time-dependent data
 642    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 643    subroutine storage_delete_all_td()
 644
 645       implicit none
 646
 647       ! Local variables
 648       integer :: funit
 649       logical :: is_used
 650       character (len=64) :: fname
 651       type (head_node), pointer :: name_cursor
 652       type (data_node), pointer :: data_cursor, next_cursor
 653
 654       ! We'll first see if there is a list for this fieldname
 655       name_cursor => head
 656       do while (associated(name_cursor))
 657
 658          data_cursor => name_cursor%fieldlist_head
 659          do while ( associated(data_cursor) )
 660             if ( is_time_dependent(data_cursor%fg_data) ) then
 661
 662                if (data_cursor%filenumber > 0) then
 663                   do funit=10,100
 664                      inquire(unit=funit, opened=is_used)
 665                      if (.not. is_used) exit
 666                   end do
 667                   write(fname,'(i9.9,a2,i3.3)') data_cursor%filenumber,'.p',my_proc_id
 668                   open(funit,file=trim(fname),form='unformatted',status='old')
 669                   close(funit,status='delete')
 670                else
 671                   call remove_index(data_cursor%heap_index)
 672                   memsize = memsize - size(data_cursor%fg_data%r_arr)
 673                   deallocate(data_cursor%fg_data%r_arr)
 674                end if
 675                if (associated(data_cursor%fg_data%valid_mask)) call bitarray_destroy(data_cursor%fg_data%valid_mask)
 676                nullify(data_cursor%fg_data%valid_mask)
 677                if (associated(data_cursor%fg_data%modified_mask)) call bitarray_destroy(data_cursor%fg_data%modified_mask)
 678                nullify(data_cursor%fg_data%modified_mask)
 679
 680                ! We should handle individual cases, that way we can deal with a list
 681                !   that has both time independent and time dependent nodes in it.
 682
 683                ! Only item in the list
 684                if (.not. associated(data_cursor%next) .and. &
 685                    .not. associated(data_cursor%prev)) then
 686                   next_cursor => null()
 687                   nullify(name_cursor%fieldlist_head)
 688                   nullify(name_cursor%fieldlist_tail)
 689                   deallocate(data_cursor)
 690 ! DO WE REMOVE THIS HEADER NODE AT THIS POINT?
 691
 692                ! Head of the list
 693                else if (.not. associated(data_cursor%prev)) then
 694                   name_cursor%fieldlist_head => data_cursor%next
 695                   next_cursor => data_cursor%next
 696                   nullify(data_cursor%next%prev)
 697                   deallocate(data_cursor)
 698
 699                ! Tail of the list
 700                else if (.not. associated(data_cursor%next)) then
 701 ! THIS CASE SHOULD PROBABLY NOT OCCUR
 702                   name_cursor%fieldlist_tail => data_cursor%prev
 703                   next_cursor => null()
 704                   nullify(data_cursor%prev%next)
 705                   deallocate(data_cursor)
 706
 707                ! Middle of the list
 708                else
 709 ! THIS CASE SHOULD PROBABLY NOT OCCUR
 710                   next_cursor => data_cursor%next
 711                   data_cursor%prev%next => data_cursor%next
 712                   data_cursor%next%prev => data_cursor%prev
 713                   deallocate(data_cursor)
 714
 715                end if
 716
 717             end if
 718             data_cursor => next_cursor
 719          end do
 720
 721          name_cursor => name_cursor%next
 722       end do
 723
 724    end subroutine storage_delete_all_td
 725
 726
 727    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 728    ! Name: storage_get_levels
 729    !
 730    ! Purpose: Returns a list of all levels for the field indicated in the_header.
 731    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 732    subroutine storage_get_levels(the_header, list)
 733
 734       implicit none
 735
 736       ! Arguments
 737       integer, pointer, dimension(:) :: list
 738       type (fg_input), intent(in) :: the_header
 739
 740       ! Local variables
 741       integer :: n
 742       type (head_node), pointer :: name_cursor
 743       type (data_node), pointer :: data_cursor
 744
 745       if (associated(list)) deallocate(list)
 746       nullify(list)
 747
 748       ! We'll first see if there is a list for this header
 749       name_cursor => head
 750       do while (associated(name_cursor))
 751          if (primary_cmp(name_cursor%fg_data, the_header) == EQUAL) exit
 752          name_cursor => name_cursor%next
 753       end do
 754
 755       if (.not. associated(name_cursor)) return
 756
 757       n = 0
 758       ! At this point, name_cursor points to a valid head node for fieldname
 759       data_cursor => name_cursor%fieldlist_head
 760       do while ( associated(data_cursor) )
 761          n = n + 1
 762          if (.not. associated(data_cursor%next)) exit
 763          data_cursor => data_cursor%next
 764       end do
 765
 766       if (n > 0) allocate(list(n))
 767
 768       n = 1
 769       do while ( associated(data_cursor) )
 770          list(n) = get_level(data_cursor%fg_data)
 771          n = n + 1
 772          data_cursor => data_cursor%prev
 773       end do
 774
 775    end subroutine storage_get_levels
 776
 777
 778    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 779    ! Name: storage_delete_all
 780    !
 781    ! Purpose: Deletes all data, both time-independent and time-dependent.
 782    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 783    subroutine storage_delete_all()
 784
 785       implicit none
 786
 787       ! Local variables
 788       integer :: funit
 789       logical :: is_used
 790       character (len=64) :: fname
 791       type (head_node), pointer :: name_cursor
 792       type (data_node), pointer :: data_cursor
 793
 794       ! We'll first see if there is already a list for this fieldname
 795       name_cursor => head
 796       do while (associated(name_cursor))
 797
 798          if (associated(name_cursor%fieldlist_head)) then
 799             data_cursor => name_cursor%fieldlist_head
 800             do while ( associated(data_cursor) )
 801                name_cursor%fieldlist_head => data_cursor%next
 802
 803                if (data_cursor%filenumber > 0) then
 804                   do funit=10,100
 805                      inquire(unit=funit, opened=is_used)
 806                      if (.not. is_used) exit
 807                   end do
 808                   write(fname,'(i9.9,a2,i3.3)') data_cursor%filenumber,'.p',my_proc_id
 809                   open(funit,file=trim(fname),form='unformatted',status='old')
 810                   close(funit,status='delete')
 811                else
 812                   call remove_index(data_cursor%heap_index)
 813                   memsize = memsize - size(data_cursor%fg_data%r_arr)
 814                   deallocate(data_cursor%fg_data%r_arr)
 815                end if
 816                if (associated(data_cursor%fg_data%valid_mask)) call bitarray_destroy(data_cursor%fg_data%valid_mask)
 817                nullify(data_cursor%fg_data%valid_mask)
 818                if (associated(data_cursor%fg_data%modified_mask)) call bitarray_destroy(data_cursor%fg_data%modified_mask)
 819                nullify(data_cursor%fg_data%modified_mask)
 820
 821                deallocate(data_cursor)
 822                data_cursor => name_cursor%fieldlist_head
 823             end do
 824          end if
 825
 826          head => name_cursor%next
 827          deallocate(name_cursor)
 828          name_cursor => head
 829       end do
 830
 831       nullify(head)
 832       nullify(tail)
 833
 834       call heap_destroy()
 835
 836    end subroutine storage_delete_all
 837
 838
 839    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 840    ! Name: storage_get_all_headers
 841    !
 842    ! Purpose:
 843    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 844    subroutine storage_get_all_headers(header_list)
 845
 846       implicit none
 847
 848       ! Arguments
 849       type (fg_input), pointer, dimension(:) :: header_list
 850
 851       ! Local variables
 852       integer :: nheaders
 853       type (head_node), pointer :: name_cursor
 854
 855       nullify(header_list)
 856
 857       ! First find out how many time-dependent headers there are
 858       name_cursor => head
 859       nheaders = 0
 860       do while (associated(name_cursor))
 861          if (associated(name_cursor%fieldlist_head)) then
 862             if (.not. is_mask_field(name_cursor%fieldlist_head%fg_data)) then
 863                nheaders = nheaders + 1
 864             end if
 865          end if
 866          name_cursor => name_cursor%next
 867       end do
 868
 869       allocate(header_list(nheaders))
 870
 871       name_cursor => head
 872       nheaders = 0
 873       do while (associated(name_cursor))
 874          if (associated(name_cursor%fieldlist_head)) then
 875             if (.not. is_mask_field(name_cursor%fieldlist_head%fg_data)) then
 876                nheaders = nheaders + 1
 877                call dup(name_cursor%fieldlist_head%fg_data, header_list(nheaders))
 878             end if
 879          end if
 880          name_cursor => name_cursor%next
 881       end do
 882
 883    end subroutine storage_get_all_headers
 884
 885
 886    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 887    ! Name: storage_get_all_td_headers
 888    !
 889    ! Purpose:
 890    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 891    subroutine storage_get_td_headers(header_list)
 892
 893       implicit none
 894
 895       ! Arguments
 896       type (fg_input), pointer, dimension(:) :: header_list
 897
 898       ! Local variables
 899       integer :: nheaders
 900       type (head_node), pointer :: name_cursor
 901
 902       nullify(header_list)
 903
 904       ! First find out how many time-dependent headers there are
 905       name_cursor => head
 906       nheaders = 0
 907       do while (associated(name_cursor))
 908          if (associated(name_cursor%fieldlist_head)) then
 909             if (is_time_dependent(name_cursor%fieldlist_head%fg_data) .and. &
 910                 .not. is_mask_field(name_cursor%fieldlist_head%fg_data)) then
 911                nheaders = nheaders + 1
 912             end if
 913          end if
 914          name_cursor => name_cursor%next
 915       end do
 916
 917       allocate(header_list(nheaders))
 918
 919       name_cursor => head
 920       nheaders = 0
 921       do while (associated(name_cursor))
 922          if (associated(name_cursor%fieldlist_head)) then
 923             if (is_time_dependent(name_cursor%fieldlist_head%fg_data) .and. &
 924                 .not. is_mask_field(name_cursor%fieldlist_head%fg_data)) then
 925                nheaders = nheaders + 1
 926                call dup(name_cursor%fieldlist_head%fg_data, header_list(nheaders))
 927             end if
 928          end if
 929          name_cursor => name_cursor%next
 930       end do
 931
 932    end subroutine storage_get_td_headers
 933
 934
 935    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 936    ! Name: storage_print_fields
 937    !
 938    ! Purpose:
 939    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 940    subroutine storage_print_fields()
 941
 942       use list_module
 943       use stringutil
 944
 945       implicit none
 946
 947       ! Local variables
 948       integer :: i, j, k, lmax, n_fields, n_levels, max_levels, itemp
 949       logical, allocatable, dimension(:,:) :: field_has_level
 950       integer, allocatable, dimension(:) :: all_levels
 951       integer, pointer, dimension(:) :: ilevels
 952       character (len=128), allocatable, dimension(:) :: fieldname_list
 953       character (len=9) :: ctemp
 954       type (fg_input), pointer, dimension(:) :: header_list
 955
 956       type (list) :: all_levs
 957
 958       !CWH Initialize local pointer variables
 959       nullify(ilevels)
 960       nullify(header_list)     !MGD initialization for header_list should not be necessary
 961
 962       call list_init(all_levs)
 963       call storage_get_td_headers(header_list)
 964       n_fields = size(header_list)
 965
 966       allocate(fieldname_list(n_fields))
 967
 968       max_levels = 0
 969
 970       do i=1,n_fields
 971          fieldname_list(i) = header_list(i)%header%field
 972          call storage_get_levels(header_list(i), ilevels)
 973          do j=1,size(ilevels)
 974             if (.not. list_search(all_levs, ikey=ilevels(j), ivalue=itemp)) then
 975                call list_insert(all_levs, ikey=ilevels(j), ivalue=ilevels(j))
 976             end if
 977          end do
 978          n_levels = size(ilevels)
 979          if (n_levels > max_levels) max_levels = n_levels
 980          if (associated(ilevels)) deallocate(ilevels)
 981       end do
 982
 983       max_levels = list_length(all_levs)
 984
 985       allocate(all_levels(max_levels))
 986       allocate(field_has_level(n_fields,max_levels))
 987
 988       field_has_level(:,:) = .false.
 989
 990       lmax = 0
 991       do i=1,n_fields
 992          call storage_get_levels(header_list(i), ilevels)
 993          n_levels = size(ilevels)
 994          do j=1,n_levels
 995             do k=1,lmax
 996                if (all_levels(k) == ilevels(j)) exit
 997             end do
 998             if (k > lmax) then
 999                all_levels(k) = ilevels(j)
1000                lmax = lmax + 1
1001             end if
1002             field_has_level(i,k) = .true.
1003          end do
1004          if (associated(ilevels)) deallocate(ilevels)
1005       end do
1006
1007       call mprintf(.true.,DEBUG,'        .',newline=.false.)
1008       do i=1,n_fields
1009          write(ctemp,'(a9)') fieldname_list(i)(1:9)
1010          call right_justify(ctemp,9)
1011          call mprintf(.true.,DEBUG,ctemp,newline=.false.)
1012       end do
1013       call mprintf(.true.,DEBUG,' ',newline=.true.)
1014       do j=1,max_levels
1015          write(ctemp,'(i9)') all_levels(j)
1016          call mprintf(.true.,DEBUG,'%s ',s1=ctemp,newline=.false.)
1017          do i=1,n_fields
1018             if (field_has_level(i,j)) then
1019                call mprintf(.true.,DEBUG,'        X',newline=.false.)
1020             else
1021                call mprintf(.true.,DEBUG,'        -',newline=.false.)
1022             end if
1023          end do
1024          call mprintf(.true.,DEBUG,' ',newline=.true.)
1025       end do
1026
1027       deallocate(all_levels)
1028       deallocate(field_has_level)
1029       deallocate(fieldname_list)
1030       deallocate(header_list)
1031
1032       call list_destroy(all_levs)
1033
1034    end subroutine storage_print_fields
1035
1036
1037    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1038    ! Name: find_missing_values
1039    !
1040    ! Purpose:
1041    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1042    subroutine find_missing_values()
1043
1044       implicit none
1045
1046       ! Local variables
1047       integer :: i, j
1048       logical :: found_missing
1049       type (head_node), pointer :: name_cursor
1050       type (data_node), pointer :: data_cursor
1051
1052       found_missing = .false.
1053
1054       name_cursor => head
1055       do while (associated(name_cursor))
1056
1057          if (associated(name_cursor%fieldlist_head)) then
1058             data_cursor => name_cursor%fieldlist_head
1059             do while ( associated(data_cursor) )
1060                if (.not. associated(data_cursor%fg_data%valid_mask)) then
1061                   call mprintf(.true.,INFORM, &
1062                                'Field %s does not have a valid mask and will not be checked for missing values', &
1063                                s1=data_cursor%fg_data%header%field)
1064                else
1065                   ILOOP: do i=1,data_cursor%fg_data%header%dim1(2)-data_cursor%fg_data%header%dim1(1)+1
1066                   JLOOP: do j=1,data_cursor%fg_data%header%dim2(2)-data_cursor%fg_data%header%dim2(1)+1
1067                      if (.not. bitarray_test(data_cursor%fg_data%valid_mask,i,j)) then
1068                         found_missing = .true.
1069                         call mprintf(.true.,WARN,'Field %s has missing values at level %i at (i,j)=(%i,%i)', &
1070                                      s1=data_cursor%fg_data%header%field, &
1071                                      i1=data_cursor%fg_data%header%vertical_level, &
1072                                      i2=i+data_cursor%fg_data%header%dim1(1)-1, &
1073                                      i3=j+data_cursor%fg_data%header%dim2(1)-1)
1074                         exit ILOOP
1075                      end if
1076                   end do JLOOP
1077                   end do ILOOP
1078                end if
1079                data_cursor => data_cursor%next
1080             end do
1081          end if
1082
1083          name_cursor => name_cursor%next
1084       end do
1085
1086       call mprintf(found_missing,ERROR,'Missing values encountered in interpolated fields. Stopping.')
1087
1088    end subroutine find_missing_values
1089
1090
1091    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1092    ! Name: storage_print_headers
1093    !
1094    ! Purpose:
1095    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1096    subroutine storage_print_headers()
1097
1098       implicit none
1099
1100       ! Local variables
1101       type (head_node), pointer :: name_cursor
1102       type (data_node), pointer :: data_cursor
1103
1104       call mprintf(.true.,DEBUG,'>>>> STORED FIELDS <<<<')
1105       call mprintf(.true.,DEBUG,'=======================')
1106
1107       ! We'll first see if there is already a list for this fieldname
1108       name_cursor => head
1109       do while (associated(name_cursor))
1110          call print_header(name_cursor%fg_data)
1111
1112          if (associated(name_cursor%fieldlist_head)) then
1113             data_cursor => name_cursor%fieldlist_head
1114             do while ( associated(data_cursor) )
1115                call mprintf(.true.,DEBUG,'  - %i', i1=get_level(data_cursor%fg_data))
1116                call mprintf(.true.,DEBUG,' ')
1117                data_cursor => data_cursor%next
1118             end do
1119          end if
1120
1121          name_cursor => name_cursor%next
1122       end do
1123
1124    end subroutine storage_print_headers
1125
1126 end module storage_module