precipitation__method_of_fragments.F90

Go to the documentation of this file.

!> @file
!! Contains the module @ref precipitation__method_of_fragments.
 
!>
!!  Module @ref precipitation__method_of_fragments
!!  provides support for creating synthetic daily precipitation
!!  given grids of monthly sum precipitation and a "fragments" file.
!!  The fragments file is generated from observations at discrete locations,
!!  the values of which range from 0 to 1, and the sum of which is 1.
!!  The fragment value is simply the daily observed precipitation value divided
!!  by the monthly sum of all observed precipitation values for that station.
!!
!!  In addition, this routine accepts another set of rainfall adjustment grids,
!!  needed in order to ensure the resulting precipitation totals fall in line with
!!  other published values, in the development case, the Rainfall Atlas of Hawaii.
 
module precipitation__method_of_fragments
 
  use iso_c_binding
  use constants_and_conversions, only  : asint, asfloat, true, false, random_start
  use data_catalog
  use data_catalog_entry
  use dictionary
  use exceptions
  use file_operations
  use kiss_random_number_generator
  use logfiles, only            : logs, log_all, log_debug
  use parameters
  use fstring
  use fstring_list
  use simulation_datetime
  use grid
  implicit none
 
  private
 
  public :: precipitation_method_of_fragments_initialize
  public :: read_daily_fragments
  public :: precipitation_method_of_fragments_calculate
 
  !> Module variable that holds the rainfall gage (zone) number
  integer (c_int), allocatable, public :: rain_gage_id(:)
 
  !> Module variable that holds the current day's rainfall fragment value
  real (c_float), allocatable, public  :: fragment_value(:)
 
  !> Module variable indicating which "simulation number" is active
  !! Only has meaning if the rainfall fragments are being applied via a predetermined
  !! sequence file
  integer (c_int), public :: simulation_number = 1
 
  !> Module variable that holds the rainfall adjustment factor
  real (c_float), allocatable, public  :: rainfall_adjust_factor(:)
 
  !> Module variable that holds a sequence of random numbers associated with the selection
  !! of the fragment set to use
  real (c_double), allocatable :: random_values(:,:)
 
  !> Module level variable used to create subsets of the FRAGMENT_SEQUENCES file
  logical (c_bool), allocatable :: sequence_selection(:)
 
  !> Module variable detemining whether fragment sequences are chosen at random or
  !! selected from an external file
  logical (c_bool) :: random_fragment_sequences = .true._c_bool
 
  !> Data structure that holds a single line of data from the input rainfall fragments file.
  type, public :: fragments_t
    integer (c_int) :: imonth
    integer (c_int) :: iraingagezone
    integer (c_int) :: ifragmentset
    real (c_float)  :: ffragmentvalue(31)
  end type fragments_t
 
  !> Pointer to a rainfall fragments data structure.
  type, public :: ptr_fragments_t
    type (fragments_t), pointer  :: pfragment => null()
  end type ptr_fragments_t
 
  !> Data structure to hold the current active rainfall fragments for
  !! a particular rain gage zone.
  type, public :: fragments_set_t
    integer (c_int) :: iraingagezone
    integer (c_int) :: inumberoffragments(12)
    integer (c_int) :: istartrecord(12)
  end type fragments_set_t
 
  !> Array of all fragments read in from the rainfall fragments file.
  type (fragments_t), allocatable, target, public       :: fragments(:)
 
  !> Subset of rainfall fragments file pointing to the currently active fragments.
  type (ptr_fragments_t), allocatable                   :: current_fragments(:,:)
 
  !> Array of fragments sets; fragments sets include indices to the start record
  !! associated with the fragment for each month; FRAGMENTS_SETS will have a
  !! number of elements equal to the number of rainfall gages in the model domain
  type (fragments_set_t), allocatable, public           :: fragments_sets(:)
 
  !> Data structure to hold static (pre-calculated) fragment selection numbers
  type, public :: fragments_sequence_t
    integer (c_int) :: sim_number
    integer (c_int) :: sim_month
    integer (c_int) :: sim_rainfall_zone
    integer (c_int) :: sim_year
    real (c_float)  :: sim_random_number
    integer (c_int) :: sim_selected_set
  end type fragments_sequence_t
 
  !> Array of fragment sequence sets
  type (fragments_sequence_t), allocatable, public  :: fragments_sequence(:)
 
  integer (c_int) :: lu_fragments_echo
 
contains
 
  !> Initialize method of fragments.
  !!
  !! This routine accesses the "RAINFALL_ZONE" gridded data object and
  !! calls the routine to read in the rainfall fragments file. Values of RAINFALL_ZONE are stored
  !! in a module variable @ref RAIN_GAGE_ID for future reference.
  !!
  !! @params[in]   lActive   2-D boolean array defining active and inactive cells
 
  subroutine precipitation_method_of_fragments_initialize( lActive )
 
    logical (c_bool), intent(in)     :: lactive(:,:)
 
    ! [ LOCALS ]
    integer (c_int)                 :: istat
    type (data_catalog_entry_t), pointer :: prainfall_zone
    type (fstring_list_t)                 :: slstring
    integer (c_int)                 :: imaxrainzones
    integer (c_int), allocatable    :: isimulationnumbers(:)
    character (len=256)                  :: error_str
 
    ! look up the simulation number associated with the desired fragment sequence set
    call cf_dict%get_values( skey="FRAGMENTS_SEQUENCE_SIMULATION_NUMBER", ivalues=isimulationnumbers )
    if ( isimulationnumbers(1) > 0 )  simulation_number = isimulationnumbers(1)
 
    ! locate the data structure associated with the gridded rainfall zone entries
    prainfall_zone => dat%find("RAINFALL_ZONE")
    if ( .not. associated(prainfall_zone) ) &
        call die("A RAINFALL_ZONE grid must be supplied in order to make use of this option.",    &
          __file__, __line__)
 
    allocate( rain_gage_id( count(lactive) ), stat=istat )
    call assert( istat == 0, "Problem allocating memory", __file__, __line__ )
 
    call prainfall_zone%getvalues()
 
    ! map the 2D array of RAINFALL_ZONE values to the vector of active cells
    rain_gage_id = pack( prainfall_zone%pGrdBase%iData, lactive )
 
    allocate( rainfall_adjust_factor( count(lactive) ), stat=istat )
    call assert( istat == 0, "Problem allocating memory", __file__, __line__ )
 
    allocate( fragment_value( count(lactive) ), stat=istat )
    call assert( istat == 0, "Problem allocating memory", __file__, __line__ )
 
 
    ! look up the name of the fragments file in the control file dictionary
    call cf_dict%get_values( skey="FRAGMENTS_DAILY_FILE", slstring=slstring )
 
    ! use the first entry in the string list slString as the filename to open for
    ! use with the daily fragments routine
 
 
    call read_daily_fragments( slstring%get(1) )
    call slstring%clear()
 
    ! look up the name of the fragments SEQUENCE file in the control file dictionary
    call cf_dict%get_values( skey="FRAGMENTS_SEQUENCE_FILE", slstring=slstring )
 
    if ( .not. ( slstring%get(1) .strequal. "<NA>" ) )  then
      call read_fragments_sequence( slstring%get(1) )
      random_fragment_sequences = .false._c_bool
      allocate ( sequence_selection( count(fragments_sequence%sim_month > 0) ), stat=istat )
      call assert( istat == 0, "Problem allocating memory", __file__, __line__ )
    endif
 
    !> Now the fragments file is in memory. Create an ancillary data structure
    !> to keep track of which records correspond to various rain zones
 
    imaxrainzones = maxval(fragments%iRainGageZone)
 
    allocate ( fragments_sets( imaxrainzones ), stat=istat )
    call assert( istat == 0, "Problem allocating memory", __file__, __line__ )
 
    if ( .not. allocated(current_fragments) ) then
      allocate (current_fragments( imaxrainzones, 1 ), stat=istat, errmsg=error_str )
      call assert( istat == 0, "Problem allocating memory, stat="//ascharacter(istat)  &
        //"; msg: "//trim(error_str), __file__, __line__ )
    endif
 
    if ( .not. allocated( random_values) ) then
      allocate (random_values( imaxrainzones, 1 ), stat=istat )
      call assert( istat == 0, "Problem allocating memory, stat="//ascharacter(istat)  &
        //"; msg: "//trim(error_str), __file__, __line__ )
    endif
 
    if ( random_fragment_sequences ) call kiss64_initialize( random_start )
 
    call process_fragment_sets()
 
    open( newunit=lu_fragments_echo, file="Fragments_as_implemented_by_SWB.csv")
    write( lu_fragments_echo, fmt="(a, 30('fragment, '),'fragment')")          &
        "Simulation_Number, Month, Rain_Zone, Year, Random_Number, Fragment_Set,"
 
  end subroutine precipitation_method_of_fragments_initialize
 
  !--------------------------------------------------------------------------------------------------
 
  subroutine read_daily_fragments( sFilename )
 
    character (len=*), intent(in)    :: sfilename
 
    ! [ LOCALS ]
    character (len=512)   :: srecord, ssubstring
    integer (c_int)  :: istat
    integer (c_int)  :: icount
    integer (c_int)  :: iindex
    integer (c_int)  :: last_zone
    integer (c_int)  :: last_fragment
    integer (c_int)  :: last_month
    integer (c_int)  :: inumlines
    real (c_float)   :: ftempvalue
    type (ascii_file_t), allocatable :: fragments_file
 
    allocate(fragments_file)
 
    call fragments_file%open( sfilename = sfilename,         &
                              scommentchars = "#%!",         &
                              sdelimiters = "WHITESPACE",    &
                              lhasheader = .false._c_bool )
 
    inumlines = fragments_file%numLines()
 
    allocate(  fragments( inumlines ), stat=istat )
    call assert( istat == 0, "Problem allocating memory for fragments table", __file__, __line__ )
 
    icount = 0
    last_zone = 1
    last_fragment = 0
    last_month = 1
 
    do
 
      ! read in next line of file
      srecord = fragments_file%readLine()
 
      if ( fragments_file%isEOF() ) exit
 
      icount = icount + 1
 
      ! read in month number
      call chomp(srecord, ssubstring, fragments_file%sDelimiters )
 
      if ( len_trim(ssubstring) == 0 )                                    &
        call die( "Missing month number in the daily fragments file",     &
          __file__, __line__, "Problem occured on line number "        &
          //ascharacter(fragments_file%currentLineNum() )                 &
          //" of file "//dquote(sfilename) )
 
      fragments(icount)%iMonth = asint(ssubstring)
 
      ! reset the counter tracking the previous fragment ID
      if ( fragments(icount)%iMonth == (last_month + 1) ) then
        last_fragment = 0
        last_zone = 0
      endif
 
      if ( fragments(icount)%iMonth < last_month )                             &
        call die( "Out-of-order month value in the daily fragments file",      &
          __file__, __line__, "Problem occured on line number "             &
          //ascharacter(fragments_file%currentLineNum() )                      &
          //" of file "//dquote(sfilename) )
 
      last_month = fragments(icount)%iMonth
 
      ! read in rain gage zone
      call chomp(srecord, ssubstring, fragments_file%sDelimiters )
 
      if ( len_trim(ssubstring) == 0 )                                         &
        call die( "Missing rain gage zone number in the daily fragments file", &
          __file__, __line__, "Problem occured on line number "             &
          //ascharacter(fragments_file%currentLineNum() )                      &
          //" of file "//dquote(sfilename) )
 
      fragments(icount)%iRainGageZone = asint(ssubstring)
 
      if ( fragments(icount)%iRainGageZone < last_zone )                            &
        call die( "Rain gage zone number out of order in the daily fragments file", &
          __file__, __line__, "Problem occured on line number "                  &
          //ascharacter(fragments_file%currentLineNum() )                           &
          //" of file "//dquote(sfilename) )
 
      ! reset the counter tracking the previous fragment ID
      if ( fragments(icount)%iRainGageZone == (last_zone + 1) )                &
        last_fragment = 0
 
      last_zone = fragments(icount)%iRainGageZone
 
      ! read in fragment set number for this zone
      call chomp(srecord, ssubstring, fragments_file%sDelimiters )
 
      if ( len_trim(ssubstring) == 0 )                                         &
        call die( "Missing fragment set number in the daily fragments file",   &
          __file__, __line__, "Problem occured on line number "             &
          //ascharacter(fragments_file%currentLineNum() )                      &
          //" of file "//dquote(sfilename) )
 
      fragments(icount)%iFragmentSet = asint(ssubstring)
 
      if ( fragments(icount)%iFragmentSet /= (last_fragment + 1) )                           &
        call die( "Missing or out-of-order fragment value in the daily fragments file",      &
          __file__, __line__, "Problem occured on line number "                           &
          //ascharacter(fragments_file%currentLineNum() )                                    &
          //" of file "//dquote(sfilename) )
 
      last_fragment = fragments(icount)%iFragmentSet
 
      ! read in fragments for each day of the current month
      do iindex = 1, 31
 
        ! read in fragment for given day of month
        call chomp(srecord, ssubstring, fragments_file%sDelimiters )
 
        if ( len_trim(ssubstring) == 0 )                                       &
          call die( "Missing fragment value in the daily fragments file",      &
            __file__, __line__, "Problem occured on line number "           &
            //ascharacter(fragments_file%currentLineNum() )                    &
            //" of file "//dquote(sfilename) )
 
        ftempvalue = asfloat( ssubstring )
 
        ! This substitution is needed to prevent "-9999" or "9999" values embedded in the fragments file
        ! from creeping into calculations. In the event of a "9999", the appropriate substitution is
        ! zero, since the previous fragments for the month to this point should already sum to 1.0
        if ( ( ftempvalue < 0.0_c_float ) .or. ( ftempvalue > 1.0_c_float ) ) then
          fragments(icount)%fFragmentValue(iindex) = 0.0_c_float
        else
          fragments(icount)%fFragmentValue(iindex) = ftempvalue
        endif
 
      enddo
 
      if ( fragments(icount)%iMonth == 2 )  call normalize_february_fragment_sequence( icount )
 
    enddo
 
    call logs%write("Maximum rain gage zone number: "//ascharacter(maxval(fragments%iRainGageZone)), &
      itab=31, ilinesafter=1, iloglevel=log_all)
 
  end subroutine read_daily_fragments
 
!--------------------------------------------------------------------------------------------------
 
  !> after fragments file has been read in, iterate over a set of rainfall fragments
  ! and keep track of the index values that correspond with changes
  ! in month and rain gage numbers
  subroutine process_fragment_sets()
 
    integer (c_int)   :: iCount
    integer (c_int)   :: iIndex
    integer (c_int)   :: iRainGageZone
    integer (c_int)   :: iPreviousRainGageZone
    integer (c_int)   :: iMonth
    integer (c_int)   :: iPreviousMonth
    character (len=10)     :: sBuf0
    character (len=10)     :: sBuf1
    character (len=12)     :: sBuf2
    character (len=10)     :: sBuf3
    character (len=52)     :: sBuf4
 
    ! this counter is used to accumulate the number of fragments associated with the
    ! current raingage zone/month combination
    icount = 1
 
    iraingagezone = fragments( lbound( fragments, 1) )%iRainGageZone
    ipreviousraingagezone = iraingagezone
    ipreviousmonth = fragments( lbound( fragments, 1) )%iMonth
 
    ! at this point, iRainGageZone should be 1, and iPreviousMonth should be 1,
    ! assuming that the fragments were sorted properly upon input
 
    ! populate the first record of FRAGMENT_SETS
    fragments_sets( iraingagezone )%iRainGageZone = iraingagezone
    fragments_sets( iraingagezone )%iStartRecord(ipreviousmonth) = lbound( fragments, 1)
 
    ! now iterate through *all* fragments, keeping track of the starting record for each new rainfall gage
    ! zone number
    do iindex = lbound( fragments, 1) + 1, ubound( fragments, 1 )
 
      iraingagezone = fragments(iindex)%iRainGageZone
      imonth = fragments(iindex)%iMonth
 
      if ( iraingagezone /= ipreviousraingagezone ) then
        ! the previous record was the last one associated with the previous
        ! rainfall gage zone; do not count the current record as part of the
        ! collection of records associated with previous zone
 
        fragments_sets( ipreviousraingagezone )%iNumberOfFragments(ipreviousmonth) = icount
 
 
        fragments_sets( iraingagezone )%iRainGageZone = iraingagezone
        fragments_sets( iraingagezone )%iStartRecord(imonth) = iindex
        ! ! need to handle the last fragment set as a special case
        ! FRAGMENTS_SETS( iRainGageZone )%iNumberOfFragments(iMonth) = iCount
        icount = 1
 
      else
        icount = icount + 1
      endif
 
      ipreviousmonth = imonth
      ipreviousraingagezone = iraingagezone
 
    enddo
 
    ! need to handle the last month of the last fragment set as a special case
    fragments_sets( iraingagezone )%iNumberOfFragments(imonth) = icount
 
    call logs%write("### Summary of fragment sets in memory ###", &
       iloglevel=log_all, ilinesbefore=1, ilinesafter=1, lecho=false )
    call logs%write("gage number | month      | start index  | num records ")
    call logs%write("----------- | ---------- | ------------ | ------------")
    do iindex=1, ubound( fragments_sets, 1)
      do imonth=1,12
        write (sbuf0, fmt="(i10)") iindex
        write (sbuf1, fmt="(i10)") imonth
        write (sbuf2, fmt="(i12)") fragments_sets(iindex)%iStartRecord(imonth)
        write (sbuf3, fmt="(i10)") fragments_sets(iindex)%iNumberOfFragments(imonth)
        write (sbuf4, fmt="(a10,'  | ', a10,' | ', a12,' | ',a10)") adjustl(sbuf0),     &
               adjustl(sbuf1), adjustl(sbuf2), adjustl(sbuf3)
        call logs%write( sbuf4 )
      enddo
    end do
 
  end subroutine process_fragment_sets
 
!--------------------------------------------------------------------------------------------------
 
  !> eliminate rainfall on the 29th day of February; bump up all other values to ensure sum = 1
  subroutine normalize_february_fragment_sequence( iCount )
 
    integer (c_int), intent(in)   :: iCount
 
    ! [ LOCALS ]
    real (c_float) :: sum_fragments
 
    ! we only want to correct the fragment if it was actually generated during a
    ! leap year
    if (fragments(icount)%fFragmentValue(29) > 0.0_c_float) then
      sum_fragments = sum( fragments(icount)%fFragmentValue(1:28) )
 
      ! bump up all February fragments so that their 28-day sum is 1
      fragments(icount)%fFragmentValue(1:28) = fragments(icount)%fFragmentValue(1:28)   &
                                                 / sum_fragments
      ! zero out remaining values
      fragments(icount)%fFragmentValue(29:31) = 0.0_c_float
    endif
 
  end subroutine normalize_february_fragment_sequence
 
!--------------------------------------------------------------------------------------------------
 
  ! in order to compare Hawaii Water Budget results to SWB results,
  ! it is necessary to force SWB to use the same sequence of fragments
  ! that was used in the HWB simulations
  subroutine read_fragments_sequence( sFilename )
 
    character (len=*), intent(in)    :: sFilename
 
    ! [ LOCALS ]
    character (len=512)   :: sRecord, sSubstring
    integer (c_int)  :: iStat
    integer (c_int)  :: iCount
    integer (c_int)  :: iIndex
    integer (c_int)  :: iNumLines
    type (ASCII_FILE_T), allocatable :: SEQUENCE_FILE
    character (len=10)     :: sBuf0
    character (len=10)     :: sBuf1
    character (len=12)     :: sBuf2
    character (len=10)     :: sBuf3
    character (len=10)     :: sBuf4
    character (len=256)    :: sBuf5
    character (len=256)    :: error_str
    type (FSTRING_LIST_T)   :: slHeader
    integer (c_int)   :: max_rain_gage_number
    integer (c_int)   :: max_simulation_number
 
    allocate(sequence_file)
 
    call sequence_file%open( sfilename = sfilename,         &
                             scommentchars = "#%!",         &
                             sdelimiters = "WHITESPACE",    &
                             lhasheader = .true._c_bool )
 
    slheader = sequence_file%readHeader()
 
    inumlines = sequence_file%numLines()
 
    allocate(  fragments_sequence( inumlines ), stat=istat )
    call assert( istat == 0, "Problem allocating memory for fragments sequence table",    &
      __file__, __line__ )
 
    icount = 0
 
    do
 
      ! read in next line of file
      srecord = sequence_file%readLine()
 
      if ( sequence_file%isEOF() ) exit
 
      icount = icount + 1
 
      ! read in simulation number
      call chomp(srecord, ssubstring, sequence_file%sDelimiters )
 
      if ( len_trim(ssubstring) == 0 )                                              &
          call die( "Missing simulation number in the fragments sequence file",     &
            __file__, __line__, "Problem occured on line number "                &
            //ascharacter(sequence_file%currentLineNum() )                          &
            //" of file "//dquote(sfilename) )
 
      fragments_sequence(icount)%sim_number = asint(ssubstring)
 
      ! read in month
      call chomp(srecord, ssubstring, sequence_file%sDelimiters )
 
      if ( len_trim(ssubstring) == 0 )                                              &
        call die( "Missing month number in the fragments sequence file",            &
          __file__, __line__, "Problem occured on line number "                  &
          //ascharacter(sequence_file%currentLineNum() )                            &
          //" of file "//dquote(sfilename) )
 
      fragments_sequence(icount)%sim_month = asint(ssubstring)
 
      ! read in rainfall zone
      call chomp(srecord, ssubstring, sequence_file%sDelimiters )
 
      if ( len_trim(ssubstring) == 0 )                                              &
        call die( "Missing rainfall zone number in the fragments sequence file",    &
          __file__, __line__, "Problem occured on line number "                  &
          //ascharacter(sequence_file%currentLineNum() )                            &
          //" of file "//dquote(sfilename) )
 
      fragments_sequence(icount)%sim_rainfall_zone = asint(ssubstring)
 
 
      ! read in sim_year
      call chomp(srecord, ssubstring, sequence_file%sDelimiters )
 
      if ( len_trim(ssubstring) == 0 )                                              &
        call die( "Missing year number in the fragments sequence file",             &
          __file__, __line__, "Problem occured on line number "                  &
          //ascharacter(sequence_file%currentLineNum() )                            &
          //" of file "//dquote(sfilename) )
 
      fragments_sequence(icount)%sim_year = asint(ssubstring)
 
      ! read in sim_random_number
      call chomp(srecord, ssubstring, sequence_file%sDelimiters )
 
      if ( len_trim(ssubstring) == 0 )                                                   &
        call die( "Missing simulation random number in the fragments sequence file",     &
          __file__, __line__, "Problem occured on line number "                       &
          //ascharacter(sequence_file%currentLineNum() )                                 &
          //" of file "//dquote(sfilename) )
 
      fragments_sequence(icount)%sim_random_number = asfloat(ssubstring)
 
      ! read in simulation selected set (fragment set selected by HWB)
      call chomp(srecord, ssubstring, sequence_file%sDelimiters )
 
      if ( len_trim(ssubstring) == 0 )                                                   &
        call die( "Missing selected fragment set number in the fragments sequence file", &
          __file__, __line__, "Problem occured on line number "                       &
          //ascharacter(sequence_file%currentLineNum() )                                 &
          //" of file "//dquote(sfilename) )
 
      fragments_sequence(icount)%sim_selected_set = asint(ssubstring)
 
    enddo
 
    max_rain_gage_number = maxval(fragments_sequence(:)%sim_rainfall_zone,1)
    max_simulation_number = maxval(fragments_sequence(:)%sim_number,1)
 
    ! idea here is that if we are reading in a sequence file, there is a good chance
    ! that the users is running multiple simulations; this will allow for a
    ! separate pointer to be established for each rain gage/simulation number combination
    ! as well as a means to keep the random value sequences separate by simulation
 
    if ( allocated(current_fragments))  deallocate(current_fragments, stat=istat, &
      errmsg=error_str)
    call assert( istat == 0, "Problem deallocating memory, stat="//ascharacter(istat)  &
      //"; msg: "//trim(error_str), __file__, __line__ )
 
    allocate(current_fragments(max_rain_gage_number, max_simulation_number), stat=istat, &
      errmsg=error_str)
    call assert( istat == 0, "Problem allocating memory, stat="//ascharacter(istat)  &
      //"; msg: "//trim(error_str), __file__, __line__ )
 
 
    if ( allocated(random_values))  deallocate(random_values, stat=istat, errmsg=error_str)
    call assert( istat == 0, "Problem deallocating memory, stat="//ascharacter(istat)  &
      //"; msg: "//trim(error_str), __file__, __line__ )
 
    allocate(random_values(max_rain_gage_number, max_simulation_number), stat=istat, &
      errmsg=error_str )
    call assert( istat == 0, "Problem allocating memory, stat="//ascharacter(istat)  &
      //"; msg: "//trim(error_str), __file__, __line__ )
 
    call logs%write("### Summary of fragment sequence sets in memory ###", &
       iloglevel=log_debug, ilinesbefore=1, ilinesafter=1, lecho=false )
    call logs%write("sim number | rainfall zone   | month  | year   | selected set ")
    call logs%write("----------- | ---------- | ------------ | ------------|------------")
    do iindex=1, ubound( fragments_sequence, 1)
      write (sbuf0, fmt="(i10)") fragments_sequence( iindex )%sim_number
      write (sbuf1, fmt="(i10)") fragments_sequence( iindex )%sim_rainfall_zone
      write (sbuf2, fmt="(i12)") fragments_sequence( iindex )%sim_month
      write (sbuf3, fmt="(i10)") fragments_sequence( iindex )%sim_year
      write (sbuf4, fmt="(i10)") fragments_sequence( iindex )%sim_selected_set
 
      write (sbuf5, fmt="(a,'  | ', a,'  |  ', a,'  |  ',a,'  |  ',a)")                        &
        adjustl(sbuf0), adjustl(sbuf1), adjustl(sbuf2), adjustl(sbuf3), adjustl(sbuf4)
      call logs%write( trim( sbuf5 ) )
    end do
 
  end subroutine read_fragments_sequence
 
!--------------------------------------------------------------------------------------------------
 
  !> Update rainfall fragments on daily basis.
  !!
  !! If called when lShuffle is TRUE:
  !! 1) update random values
  !! 2) random values are used to select the next active fragment set
  !!    for the current RainGageZone
  !!
  !! *Each* time the routine is called, the appropriate fragment is
  !! selected from the current active fragment set and is assigned
  !! to all cells that share a common RainGageZone
 
  subroutine update_fragments( lShuffle )
 
    logical (c_bool), intent(in) :: lShuffle
 
    ! [ LOCALS ]
    integer (c_int) :: rain_zone
    integer (c_int) :: iMaxRainZones
    integer (c_int) :: iMonth
    integer (c_int) :: iDay
    integer (c_int) :: iYearOfSimulation
 
    integer (c_int) :: iNumberOfFragments
    integer (c_int)  :: iStartRecord
    integer (c_int) :: iEndRecord
    integer (c_int) :: iTargetRecord
    integer (c_int) :: iUBOUND_FRAGMENTS
    integer (c_int) :: iUBOUND_CURRENT_FRAGMENTS
 
 
    imaxrainzones = maxval(fragments%iRainGageZone)
    imonth = sim_dt%curr%iMonth
    iday = sim_dt%curr%iDay
    iyearofsimulation=sim_dt%iYearOfSimulation
 
    ! equal the number of fragments are in memory
    iubound_fragments = ubound( fragments, 1)
 
    ! should equal the number of rainfall gages in model domain
    iubound_current_fragments = ubound( current_fragments, 1)
 
    ! if by chance a mismatch in shape-to-grid results in an active cell with *NO* valid
    ! rain gage ID, we need to set the entire array to zero to quash any spurious values getting in
    fragment_value = 0.0_c_float
 
    do rain_zone = 1, imaxrainzones
 
      if ( lshuffle ) then
        ! find next fragment *record*
 
        ! update the module variable RANDOM_VALUES
        call update_random_values()
 
        istartrecord = fragments_sets( rain_zone )%iStartRecord(imonth)
        inumberoffragments = fragments_sets(rain_zone)%iNumberOfFragments(imonth)
        iendrecord = istartrecord + inumberoffragments - 1
        itargetrecord = istartrecord                                           &
                           + int(random_values(rain_zone, simulation_number)   &
                           * real( inumberoffragments ))
 
        if ( ( rain_zone > iubound_current_fragments ) .or. ( itargetrecord > iubound_fragments )   &
            .or. ( rain_zone < 1 ) .or. ( itargetrecord < 1) ) then
          call logs%write("Error detected in method of fragments routine; dump of current"       &
                          //" variables follows:", ilinesbefore=1)
          call logs%write("rain_zone         : "//ascharacter(rain_zone), itab=3 )
          call logs%write("simulation_number : "//ascharacter(rain_zone), itab=3 )
          call logs%write("iStartRecord      : "//ascharacter(istartrecord), itab=3 )
          call logs%write("iNumberOfFragments: "//ascharacter(inumberoffragments), itab=3 )
          call logs%write("iEndRecord        : "//ascharacter(iendrecord), itab=3 )
          call logs%write("iTargetRecord     : "//ascharacter(itargetrecord), itab=3 )
          call logs%write("ubound(CURRENT_FRAGMENTS, 1): "//ascharacter(iubound_current_fragments), &
                          itab=3 )
          call logs%write("ubound(FRAGMENTS, 1): "//ascharacter(iubound_fragments), itab=3 )
          call logs%write("RANDOM_VALUES(rain_zone,SIMULATION_NUMBER): "       &
            //ascharacter(random_values(rain_zone,simulation_number)), itab=3 )
          call die( "Miscalculation in target record: calculated record index is out of bounds", &
            __file__, __line__ )
        endif
 
        ! reassign fragment pointer for this rain zone to the newly selected record
        current_fragments(rain_zone, simulation_number)%pFragment => fragments( itargetrecord )
 
         write(lu_fragments_echo,fmt="(4(i5,','),f10.6,',',i5,',',30(f8.3,','),f8.3)")   &
                     simulation_number,                                                  &
                     fragments( itargetrecord)%iMonth,                                   &
                     fragments( itargetrecord)%iRainGageZone,                            &
                     iyearofsimulation,                                                  &
                     random_values(rain_zone, simulation_number),                        &
                     fragments( itargetrecord)%iFragmentSet,                             &
                     fragments( itargetrecord)%fFragmentValue
 
        ! call LOGS%write( trim(sBuf), iLogLevel=LOG_DEBUG, lEcho=FALSE )
 
      endif
 
      if ( ( current_fragments( rain_zone, simulation_number )%pFragment%fFragmentValue( iday ) < 0.0 ) &
         .or. ( current_fragments( rain_zone, simulation_number )%pFragment%fFragmentValue( iday ) > 1.0 ) ) then
 
        call logs%write("Error detected in method of fragments routine; dump of current variables"  &
                        //" follows:", ilinesbefore=1, iloglevel=log_all )
        call logs%write("rain_zone:"//ascharacter(rain_zone), itab=3 )
        call logs%write("iDay: "//ascharacter(iday), itab=3 )
        call logs%write("iRainGageZone: "//ascharacter(fragments( itargetrecord)%iRainGageZone), itab=3 )
        call logs%write("iFragmentSet: "//ascharacter(fragments( itargetrecord)%iFragmentSet), itab=3 )
        call logs%write("fFragmentValue: "//ascharacter(fragments( itargetrecord)%fFragmentValue(iday) ), itab=3 )
 
      endif
 
      ! call LOGS%write("frag: "//asCharacter(rain_zone)//"  day: "//asCharacter(iDay) &
      !    //"  value: "//asCharacter( CURRENT_FRAGMENTS( rain_zone )%pFragment%fFragmentValue( iDay ) ), &
      !    lEcho=FALSE )
 
      ! now place current days' fragment value into the matching cells
      where ( rain_gage_id == rain_zone )
 
        fragment_value = current_fragments( rain_zone, simulation_number )%pFragment%fFragmentValue( iday )
 
      endwhere
 
    enddo
 
  end subroutine update_fragments
 
!--------------------------------------------------------------------------------------------------
 
  subroutine update_random_values()
 
    ! [ LOCALS ]
    integer (c_int) :: iIndex, iIndex2
    logical (c_bool) :: lSequenceSelection
 
    if ( random_fragment_sequences ) then
 
      do iindex2=1,size(random_values,1)
 
      !call random_number( RANDOM_VALUES )
        random_values(iindex2, simulation_number) = kiss64_uniform_rng()
 
      enddo
 
    else
 
      random_values(:,simulation_number) = -9999999.9
 
      do iindex=1, size(fragments_sequence%sim_month, 1)
 
        lsequenceselection =       ( fragments_sequence(iindex)%sim_month == sim_dt%curr%iMonth )            &
                             .and. ( fragments_sequence(iindex)%sim_year == sim_dt%iYearOfSimulation )       &
                             .and. ( fragments_sequence(iindex)%sim_number == simulation_number )
 
        if ( .not. lsequenceselection ) cycle
 
        do iindex2=1,size(random_values,1)
 
          if ( fragments_sequence( iindex )%sim_rainfall_zone == iindex2 ) then
 
            random_values( iindex2, simulation_number ) = fragments_sequence( iindex )%sim_random_number
            exit
 
          endif
 
        enddo
 
      enddo
 
    endif
 
    if (any( random_values(:, simulation_number) < 0.0 ) ) then
 
      call logs%write("Error detected in method of fragments routine - random values " &
        //"not found in sequence file for rainfall zone(s):", ilinesbefore=1)
      do iindex=1,size(random_values, 1)
        if ( random_values(iindex, simulation_number) < 0.0 )  &
          call logs%write("simulation number, rainfall zone: "                 &
            //trim(ascharacter(simulation_number))//", "//trim(ascharacter(iindex)), itab=3 )
      enddo
 
    endif
 
  end subroutine update_random_values
 
!--------------------------------------------------------------------------------------------------
 
  subroutine precipitation_method_of_fragments_calculate( lActive )
 
    use ieee_exceptions, only        : ieee_invalid,          &
                                       ieee_divide_by_zero,   &
                                       ieee_overflow,         &
                                       ieee_set_halting_mode, &
                                       ieee_get_halting_mode
    use ieee_arithmetic, only        : ieee_is_nan
 
    logical (c_bool), intent(in)     :: lactive(:,:)
 
    ! [ LOCALS ]
    logical (c_bool), save       :: lfirstcall = true
 
    type (data_catalog_entry_t), pointer :: prainfall_adjust_factor
    logical :: old_invalid, old_div0, old_overflow
 
    ! attempting to keep code from choking on NaN comparisons during 
    ! file initialization
    call ieee_get_halting_mode(ieee_invalid,        old_invalid)
    call ieee_get_halting_mode(ieee_divide_by_zero, old_div0)
    call ieee_get_halting_mode(ieee_overflow,       old_overflow)
    
    ! Mask traps before calling into NetCDF/HDF5:
    call ieee_set_halting_mode(ieee_invalid,        .false.)
    call ieee_set_halting_mode(ieee_divide_by_zero, .false.)
    call ieee_set_halting_mode(ieee_overflow,       .false.)
 
    !! if it is the first day of the month, update the rainfall adjustment factor grid
    !! and update the fragments
    if ( sim_dt%curr%iDay == 1 .or. lfirstcall ) then
 
      ! locate the data structure associated with the gridded rainfall adjustment factor
      prainfall_adjust_factor => dat%find("RAINFALL_ADJUST_FACTOR")
 
      if ( .not. associated(prainfall_adjust_factor) ) &
          call die("A RAINFALL_ADJUST_FACTOR grid must be supplied in order to make use"     &
                   //" of this option.", __file__, __line__)
 
      call prainfall_adjust_factor%getvalues( sim_dt%curr )
 
      ! map the 2D array of RAINFALL_ADJUST_FACTOR values to the vector of active cells
      rainfall_adjust_factor = pack( prainfall_adjust_factor%pGrdBase%rData, lactive )
 
      call assert(.not. any(ieee_is_nan(rainfall_adjust_factor)),  &
                 smessage="NaN detected in rainfall adjustment factor grid.") 
 
      call update_fragments( lshuffle = true)
      lfirstcall = false
 
    else
 
      call update_fragments( lshuffle = false )
 
    endif
 
    ! Restore your preferred behavior afterwards:
    call ieee_set_halting_mode(ieee_invalid,        old_invalid)
    call ieee_set_halting_mode(ieee_divide_by_zero, old_div0)
    call ieee_set_halting_mode(ieee_overflow,       old_overflow)
 
  end subroutine precipitation_method_of_fragments_calculate
 
!--------------------------------------------------------------------------------------------------
 
end module precipitation__method_of_fragments