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