module TimeArraySeriesModule

  use ArrayReadersModule, only: ReadArray
  use BlockParserModule,  only: BlockParserType
  use ConstantsModule,    only: LINELENGTH, UNDEFINED, STEPWISE, LINEAR,        &
                                LENTIMESERIESNAME, DZERO, DONE
  use GenericUtilitiesModule,   only: is_same
  use InputOutputModule,  only: GetUnit, openfile
  use KindModule,         only: DP, I4B
  use ListModule,         only: ListType, ListNodeType
  use SimVariablesModule, only: errmsg
  use SimModule,          only: count_errors, store_error, store_error_unit
  use TimeArrayModule,    only: TimeArrayType, ConstructTimeArray, &
                                AddTimeArrayToList, CastAsTimeArrayType, &
                                GetTimeArrayFromList
  use BaseDisModule,      only: DisBaseType
  use, intrinsic :: iso_fortran_env, only: IOSTAT_END

  implicit none
  private
  public  :: TimeArraySeriesType, ConstructTimeArraySeries, &
             CastAsTimeArraySeriesType, GetTimeArraySeriesFromList

  type TimeArraySeriesType
    ! -- Public members
    character(len=LENTIMESERIESNAME), public :: Name = ''
    ! -- Private members
    integer(I4B), private :: inunit = 0
    integer(I4B), private :: iout = 0
    integer(I4B), private :: iMethod = UNDEFINED
    real(DP), private :: sfac = DONE
    character(len=LINELENGTH), private :: dataFile = ''
    logical, private :: autoDeallocate = .true.
    type(ListType), pointer, private :: list => null()
    class(DisBaseType), pointer, private :: dis => null()
    type(BlockParserType), private :: parser
  contains
    ! -- Public procedures
    procedure, public :: tas_init
    procedure, public :: GetAverageValues
    procedure, public :: GetInunit
    procedure, public :: da => tas_da
    ! -- Private procedures
    procedure, private :: get_integrated_values
    procedure, private :: get_latest_preceding_node
    procedure, private :: get_values_at_time
    procedure, private :: get_surrounding_records
    procedure, private :: read_next_array
    procedure, private :: DeallocateBackward
  end type TimeArraySeriesType

contains

  ! -- Constructor for TimeArraySeriesType

  subroutine ConstructTimeArraySeries(newTas, filename)
! ******************************************************************************
! ConstructTimeArraySeries -- Allocate a new TimeArraySeriesType object.
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    type(TimeArraySeriesType), pointer, intent(out) :: newTas
    character(len=*), intent(in) :: filename
    ! -- local
    logical :: lex
! ------------------------------------------------------------------------------
    ! formats
    10 format('Error: Time-array-series file "',a,'" does not exist.')
    !
    ! -- Allocate a new object of type TimeArraySeriesType
    allocate(newTas)
    allocate(newTas%list)
    !
    ! -- Ensure that input file exists
    inquire(file=filename,exist=lex)
    if (.not. lex) then
      write(errmsg,10)trim(filename)
      call store_error(errmsg, terminate=.TRUE.)
    endif
    newTas%datafile = filename
    !
    return
  end subroutine ConstructTimeArraySeries

  ! -- Public procedures

  subroutine tas_init(this, fname, dis, iout, tasname, autoDeallocate)
! ******************************************************************************
! tas_init -- initialize the time array series
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(TimeArraySeriesType), intent(inout) :: this
    character(len=*), intent(in) :: fname
    class(DisBaseType), pointer, intent(inout) :: dis
    integer(I4B), intent(in) :: iout
    character(len=*), intent(inout) :: tasname
    logical, optional,          intent(in)    :: autoDeallocate
    ! -- local
    integer(I4B) :: istatus
    integer(I4B) :: ierr
    integer(I4B) :: inunit
    character(len=40) :: keyword, keyvalue
    logical :: found, continueread, endOfBlock
! ------------------------------------------------------------------------------
    !
    ! -- initialize some variables
    if (present(autoDeallocate)) this%autoDeallocate = autoDeallocate
    this%dataFile = fname
    allocate(this%list)
    !
    ! -- assign members
    this%dis => dis
    this%iout = iout
    !
    ! -- open time-array series input file
    inunit = GetUnit()
    this%inunit = inunit
    call openfile(inunit, 0, fname, 'TAS6')
    !
    ! -- initialize block parser
    call this%parser%Initialize(this%inunit, this%iout)
    !
    ! -- read ATTRIBUTES block
    continueread = .false.
    ierr = 0
    !
    ! -- get BEGIN line of ATTRIBUTES block
    call this%parser%GetBlock('ATTRIBUTES', found, ierr, &
      supportOpenClose=.true.)
    if (.not. found) then
      errmsg = 'Error: Attributes block not found in file: ' // &
              trim(fname)
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    endif
    !
    ! -- parse ATTRIBUTES entries
    do
      ! -- read line from input
      call this%parser%GetNextLine(endOfBlock)
      if (endOfBlock) exit
      !
      ! -- get the keyword
      call this%parser%GetStringCaps(keyword)
      !
      ! -- get the word following the keyword (the key value)
      call this%parser%GetStringCaps(keyvalue)
      select case (keyword)
      case ('NAME')
        this%Name = keyvalue
        tasname = keyvalue
      case ('METHOD')
        select case (keyvalue)
        case ('STEPWISE')
          this%iMethod = STEPWISE
        case ('LINEAR')
          this%iMethod = LINEAR
        case default
          errmsg = 'Unknown interpolation method: "' // trim(keyvalue) // '"'
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end select
      case ('AUTODEALLOCATE')
        this%autoDeallocate = (keyvalue == 'TRUE')
      case ('SFAC')
        read(keyvalue,*,iostat=istatus)this%sfac
        if (istatus /= 0) then
          errmsg = 'Error reading numeric SFAC value from "' // trim(keyvalue) &
                  // '"'
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        endif
      case default
        errmsg = 'Unknown option found in ATTRIBUTES block: "' // &
                trim(keyword) // '"'
        call store_error(errmsg)
        call this%parser%StoreErrorUnit()
      end select
    enddo
    !
    ! -- ensure that NAME and METHOD have been specified
    if (this%Name == '') then
      errmsg = 'Name not specified for time array series in file: ' // &
               trim(this%dataFile)
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    endif
    if (this%iMethod == UNDEFINED) then
      errmsg = 'Interpolation method not specified for time' // &
               ' array series in file: ' // trim(this%dataFile)
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    endif
    !
    ! -- handle any errors encountered so far
    if (count_errors()>0) then
      errmsg = 'Error(s) encountered initializing time array series from file: ' // &
               trim(this%dataFile)
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    endif
    !
    ! -- try to read first time array into linked list
    if (.not. this%read_next_array()) then
      errmsg = 'Error encountered reading time-array data from file: ' // &
               trim(this%dataFile)
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    endif
    !
    return
  end subroutine tas_init

  subroutine GetAverageValues(this, nvals, values, time0, time1)
! ******************************************************************************
! GetAverageValues -- populate an array time-weighted average value for a
!   specified time span.
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(TimeArraySeriesType), intent(inout) :: this
    integer(I4B),                    intent(in)    :: nvals
    real(DP), dimension(nvals), intent(inout) :: values
    real(DP),           intent(in)    :: time0
    real(DP),           intent(in)    :: time1
    ! -- local
    integer(I4B) :: i
    real(DP) :: timediff
! ------------------------------------------------------------------------------
    !
    timediff = time1 - time0
    if (timediff > 0) then
      call this%get_integrated_values(nvals, values, time0, time1)
      do i=1,nvals
        values(i) = values(i) / timediff
      enddo
    else
      ! -- time0 and time1 are the same, so skip the integration step.
      call this%get_values_at_time(nvals, values, time0)
    endif
    !
    return
  end subroutine GetAverageValues

  function GetInunit(this)
! ******************************************************************************
! GetInunit -- return unit number
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- return
    integer(I4B) :: GetInunit
    ! -- dummy
    class(TimeArraySeriesType) :: this
! ------------------------------------------------------------------------------
    !
    GetInunit = this%inunit
    !
    return
  end function GetInunit

  ! -- Private procedures

  subroutine get_surrounding_records(this, time, taEarlier, taLater)
! ******************************************************************************
! get_surrounding_records -- get_surrounding_records
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(TimeArraySeriesType), intent(inout) :: this
    real(DP),      intent(in)    :: time
    type(TimeArrayType), pointer, intent(inout) :: taEarlier
    type(TimeArrayType), pointer, intent(inout) :: taLater
    ! -- local
    real(DP) :: time0, time1
    type(ListNodeType), pointer :: currNode => null()
    type(ListNodeType), pointer :: node0 => null()
    type(ListNodeType), pointer :: node1 => null()
    type(TimeArrayType), pointer :: ta => null(), ta0 => null(), ta1 => null()
    class(*), pointer :: obj
! ------------------------------------------------------------------------------
    !
    taEarlier => null()
    taLater => null()
    !
    if (associated(this%list%firstNode)) then
      currNode => this%list%firstNode
    endif
    !
    ! -- If the next node is earlier than time of interest, advance along
    !    linked list until the next node is later than time of interest.
    do
      if (associated(currNode)) then
        if (associated(currNode%nextNode)) then
          obj => currNode%nextNode%GetItem()
          ta => CastAsTimeArrayType(obj)
          if (ta%taTime <= time) then
            currNode => currNode%nextNode
          else
            exit
          endif
        else
          ! -- read another array
          if (.not. this%read_next_array()) exit
        endif
      else
        exit
      endif
    enddo
    !
    if (associated(currNode)) then
      !
      ! -- find earlier record
      node0 => currNode
      obj => node0%GetItem()
      ta0 => CastAsTimeArrayType(obj)
      time0 = ta0%taTime
      do while (time0 > time)
        if (associated(node0%prevNode)) then
          node0 => node0%prevNode
          obj => node0%GetItem()
          ta0 => CastAsTimeArrayType(obj)
          time0 = ta0%taTime
        else
          exit
        endif
      enddo
      !
      ! -- find later record
      node1 => currNode
      obj => node1%GetItem()
      ta1 => CastAsTimeArrayType(obj)
      time1 = ta1%taTime
      do while (time1 < time)
        if (associated(node1%nextNode)) then
          node1 => node1%nextNode
          obj => node1%GetItem()
          ta1 => CastAsTimeArrayType(obj)
          time1 = ta1%taTime
        else
          ! -- get next array
          if (.not. this%read_next_array()) then
            ! -- end of file reached, so exit loop
            exit
          endif
        endif
      enddo
      !
    endif
    !
    if (time0 <= time) taEarlier => ta0
    if (time1 >= time) taLater => ta1
    !
    return
  end subroutine get_surrounding_records

  logical function read_next_array(this)
! ******************************************************************************
! read_next_array -- Read next time array from input file and append to list.
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(TimeArraySeriesType), intent(inout) :: this
    ! -- local
    integer(I4B) :: i, ierr, istart, istat, istop, lloc, nrow, ncol, nodesperlayer
    logical :: lopen, isFound
    type(TimeArrayType), pointer  :: ta => null()
! ------------------------------------------------------------------------------
    !
    istart = 1
    istat = 0
    istop = 1
    lloc = 1
    ! Get dimensions for supported discretization type
    if (this%dis%supports_layers()) then
      nodesperlayer = this%dis%get_ncpl()
      if(size(this%dis%mshape) == 3) then
        nrow = this%dis%mshape(2)
        ncol = this%dis%mshape(3)
      else
        nrow = 1
        ncol = this%dis%mshape(2)
      endif
    else
      errmsg = 'Time array series is not supported for selected discretization type.'
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    endif
    !
    read_next_array = .false.
    inquire(unit=this%inunit,opened=lopen)
    if (lopen) then
      call ConstructTimeArray(ta, this%dis)
      ! -- read a time and an array from the input file
      ! -- Get a TIME block and read the time
      call this%parser%GetBlock('TIME', isFound, ierr, &
        supportOpenClose=.false.)
      if (isFound) then
        ta%taTime = this%parser%GetDouble()
        ! -- Read the array
        call ReadArray(this%parser%iuactive, ta%taArray, this%Name, &
                        this%dis%ndim, ncol, nrow, 1, nodesperlayer, &
                        this%iout, 0, 0)
        !
        ! -- multiply values by sfac
        do i = 1, nodesperlayer
          ta%taArray(i) = ta%taArray(i) * this%sfac
        enddo
        !
        ! -- append the new time array to the list
        call AddTimeArrayToList(this%list, ta)
        read_next_array = .true.
        !
        ! -- make sure block is closed
        call this%parser%terminateblock()
      endif
    endif
    return ! Normal return
    !
    return
  end function read_next_array

  subroutine get_values_at_time(this, nvals, values, time)
! ******************************************************************************
! get_values_at_time -- Return an array of values for a specified time, same
!   units as time-series values.
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(TimeArraySeriesType), intent(inout) :: this
    integer(I4B), intent(in) :: nvals
    real(DP), dimension(nvals), intent(inout) :: values
    real(DP), intent(in) :: time ! time of interest
    ! -- local
    integer(I4B) :: i, ierr
    real(DP) :: ratio, time0, time1, timediff, timediffi, val0, val1, &
                        valdiff
    type(TimeArrayType), pointer :: taEarlier => null()
    type(TimeArrayType), pointer :: taLater => null()
    ! formats
    10 format('Error getting array at time ',g10.3, &
              ' for time-array series "',a,'"')
! ------------------------------------------------------------------------------
    !
    ierr = 0
    call this%get_surrounding_records(time,taEarlier,taLater)
    if (associated(taEarlier)) then
      if (associated(taLater)) then
        ! -- values are available for both earlier and later times
        if (this%iMethod == STEPWISE) then
          ! -- Just populate values from elements of earlier time array
          values =  taEarlier%taArray
        elseif (this%iMethod == LINEAR) then
          ! -- perform linear interpolation
          time0 = taEarlier%taTime
          time1 = taLater%tatime
          timediff = time1 - time0
          timediffi = time - time0
          if (timediff>0) then
            ratio = timediffi/timediff
          else
            ! -- should not happen if TS does not contain duplicate times
            ratio = 0.5d0
          endif
          ! -- Iterate through all elements and perform interpolation.
          do i=1,nvals
            val0 = taEarlier%taArray(i)
            val1 = taLater%taArray(i)
            valdiff = val1 - val0
            values(i) = val0 + (ratio*valdiff)
          enddo
        else
          ierr = 1
        endif
      else
        if (is_same(taEarlier%taTime, time)) then
          values = taEarlier%taArray
        else
          ! -- Only earlier time is available, and it is not time of interest;
          !    however, if method is STEPWISE, use value for earlier time.
          if (this%iMethod == STEPWISE) then
            values =  taEarlier%taArray
          else
            ierr = 1
          endif
        endif
      endif
    else
      if (associated(taLater)) then
        if (is_same(taLater%taTime, time)) then
          values = taLater%taArray
        else
          ! -- only later time is available, and it is not time of interest
          ierr = 1
        endif
      else
        ! -- Neither earlier nor later time is available.
        !    This should never happen!
        ierr = 1
      endif
    endif
    !
    if (ierr > 0) then
      write(errmsg,10)time,trim(this%Name)
      call store_error(errmsg)
      call store_error_unit(this%inunit)
    endif
    !
    return
  end subroutine get_values_at_time

  subroutine get_integrated_values(this, nvals, values, time0, time1)
! ******************************************************************************
! get_integrated_values -- Populates an array with integrated values for a
!    specified time span.  Units: (ts-value-unit)*time
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(TimeArraySeriesType), intent(inout) :: this
    integer(I4B),                    intent(in)    :: nvals
    real(DP), dimension(nvals), intent(inout) :: values
    real(DP),           intent(in)    :: time0
    real(DP),           intent(in)    :: time1
    ! -- local
    integer(I4B) :: i
    real(DP) :: area, currTime, nextTime, ratio0, ratio1, t0, &
                        t01, t1, timediff, value, value0, value1, valuediff
    logical :: ldone
    type(ListNodeType), pointer :: precNode => null()
    type(ListNodeType), pointer :: currNode => null(), nextNode => null()
    type(TimeArrayType), pointer :: currRecord => null(), nextRecord => null()
    class(*), pointer :: currObj => null(), nextObj => null()
    ! -- formats
10  format('Error encountered while performing integration', &
        ' for time-array series "',a,'" for time interval: ', &
        g12.5,' to ',g12.5)
! ------------------------------------------------------------------------------
    !
    values = DZERO
    value = DZERO
    ldone = .false.
    t1 = -DONE
    call this%get_latest_preceding_node(time0, precNode)
    if (associated(precNode)) then
      currNode => precNode
      do while (.not. ldone)
        currObj => currNode%GetItem()
        currRecord => CastAsTimeArrayType(currObj)
        currTime = currRecord%taTime
        if (currTime < time1) then
          if (.not. associated(currNode%nextNode)) then
            ! -- try to read the next array
            if (.not. this%read_next_array()) then
              write(errmsg,10)trim(this%Name),time0,time1
              call store_error(errmsg)
              call store_error_unit(this%inunit)
            endif
          endif
          if (associated(currNode%nextNode)) then
            nextNode => currNode%nextNode
            nextObj => nextNode%GetItem()
            nextRecord => CastAsTimeArrayType(nextObj)
            nextTime = nextRecord%taTime
            ! -- determine lower and upper limits of time span of interest
            !    within current interval
            if (currTime >= time0) then
              t0 = currTime
            else
              t0 = time0
            endif
            if (nextTime <= time1) then
              t1 = nextTime
            else
              t1 = time1
            endif
            ! -- For each element, find area of rectangle
            !    or trapezoid delimited by t0 and t1.
            t01 = t1 - t0
            select case (this%iMethod)
            case (STEPWISE)
              do i=1,nvals
                ! -- compute area of a rectangle
                value0 = currRecord%taArray(i)
                area = value0 * t01
                ! -- add area to integrated value
                values(i) = values(i) + area
              enddo
            case (LINEAR)
              do i=1,nvals
                ! -- compute area of a trapezoid
                timediff = nextTime - currTime
                ratio0 = (t0 - currTime) / timediff
                ratio1 = (t1 - currTime) / timediff
                valuediff = nextRecord%taArray(i) - currRecord%taArray(i)
                value0 = currRecord%taArray(i) + ratio0 * valuediff
                value1 = currRecord%taArray(i) + ratio1 * valuediff
                area = 0.5d0 * t01 * (value0 + value1)
                ! -- add area to integrated value
                values(i) = values(i) + area
              enddo
            end select
          else
            write(errmsg,10)trim(this%Name),time0,time1
            call store_error(errmsg)
            call store_error('(Probable programming error)', terminate=.TRUE.)
          endif
        else
          ! Current node time = time1 so should be done
          ldone = .true.
        endif
        !
        ! -- Are we done yet?
        if (t1 >= time1) then
          ldone = .true.
        else
          if (.not. associated(currNode%nextNode)) then
            ! -- try to read the next array
            if (.not. this%read_next_array()) then
              write(errmsg,10)trim(this%Name),time0,time1
              call store_error(errmsg)
              call this%parser%StoreErrorUnit()
            endif
          endif
          if (associated(currNode%nextNode)) then
            currNode => currNode%nextNode
          else
            write(errmsg,10)trim(this%Name),time0,time1
            call store_error(errmsg)
            call store_error('(Probable programming error)', terminate=.TRUE.)
          endif
        endif
      enddo
    endif
    !
    if (this%autoDeallocate) then
      if (associated(precNode)) then
        if (associated(precNode%prevNode))then
          call this%DeallocateBackward(precNode%prevNode)
        endif
      endif
    endif
    !
    return
  end subroutine get_integrated_values

  subroutine DeallocateBackward(this, fromNode)
! ******************************************************************************
! DeallocateBackward -- Deallocate fromNode and all previous nodes in list;
!   reassign firstNode.
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(TimeArraySeriesType),  intent(inout) :: this
    type(ListNodeType), pointer, intent(inout) :: fromNode
    !
    ! -- local
    type(ListNodeType),  pointer :: current => null()
    type(ListNodeType),  pointer :: prev => null()
    type(TimeArrayType), pointer :: ta => null()
    class(*),            pointer :: obj => null()
! ------------------------------------------------------------------------------
    !
    if (associated(fromNode)) then
      ! -- reassign firstNode
      if (associated(fromNode%nextNode)) then
        this%list%firstNode => fromNode%nextNode
      else
        this%list%firstNode => null()
      endif
      ! -- deallocate fromNode and all previous nodes
      current => fromNode
      do while (associated(current))
        prev => current%prevNode
        obj => current%GetItem()
        ta => CastAsTimeArrayType(obj)
        ! -- Deallocate the contents of this time array,
        !    then remove it from the list
        call ta%da()
        call this%list%RemoveNode(current, .true.)
        current => prev
      enddo
      fromNode => null()
    endif
    !
    return
  end subroutine DeallocateBackward

  subroutine get_latest_preceding_node(this, time, tslNode)
! ******************************************************************************
! get_latest_preceding_node -- Return pointer to ListNodeType object for the
!    node representing the latest preceding time in the time series
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(TimeArraySeriesType),  intent(inout) :: this
    real(DP),            intent(in)    :: time
    type(ListNodeType), pointer, intent(inout) :: tslNode
    ! -- local
    real(DP) :: time0
    type(ListNodeType),  pointer :: currNode => null()
    type(ListNodeType),  pointer :: node0 => null()
    type(TimeArrayType), pointer :: ta => null()
    type(TimeArrayType), pointer :: ta0 => null()
    class(*),            pointer :: obj => null()
! ------------------------------------------------------------------------------
    !
    tslNode => null()
    if (associated(this%list%firstNode)) then
      currNode => this%list%firstNode
    else
      call store_error('probable programming error in get_latest_preceding_node', &
                       terminate=.TRUE.)
    endif
    !
    continue
    ! -- If the next node is earlier than time of interest, advance along
    !    linked list until the next node is later than time of interest.
    do
      if (associated(currNode)) then
        if (associated(currNode%nextNode)) then
          obj => currNode%nextNode%GetItem()
          ta => CastAsTimeArrayType(obj)
          if (ta%taTime < time  .or. is_same(ta%taTime, time)) then
            currNode => currNode%nextNode
          else
            exit
          endif
        else
          ! -- read another record
          if (.not. this%read_next_array()) exit
        endif
      else
        exit
      endif
    enddo
    !
    if (associated(currNode)) then
      !
      ! -- find earlier record
      node0 => currNode
      obj => node0%GetItem()
      ta0 => CastAsTimeArrayType(obj)
      time0 = ta0%taTime
      do while (time0 > time)
        if (associated(node0%prevNode)) then
          node0 => node0%prevNode
          obj => node0%GetItem()
          ta0 => CastAsTimeArrayType(obj)
          time0 = ta0%taTime
        else
          exit
        endif
      enddo
    endif
    !
    if (time0 <= time) tslNode => node0
    !
    return
  end subroutine get_latest_preceding_node

  subroutine tas_da(this)
! ******************************************************************************
! tas_da -- deallocate
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(TimeArraySeriesType), intent(inout) :: this
    ! -- local
    integer :: i, n
    type(TimeArrayType), pointer :: ta => null()
! ------------------------------------------------------------------------------
    !
    ! -- Deallocate contents of each time array in list
    n = this%list%Count()
    do i=1,n
      ta => GetTimeArrayFromList(this%list, i)
      call ta%da()
    enddo
    !
    ! -- Deallocate the list of time arrays
    call this%list%Clear(.true.)
    deallocate(this%list)
    !
    return
  end subroutine tas_da

  ! -- Procedures not type-bound

  function CastAsTimeArraySeriesType(obj) result (res)
! ******************************************************************************
! CastAsTimeArraySeriesType -- Cast an unlimited polymorphic object as
!   class(TimeArraySeriesType)
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(*),   pointer, intent(inout) :: obj
    type(TimeArraySeriesType), pointer :: res
! ------------------------------------------------------------------------------
    !
    res => null()
    if (.not. associated(obj)) return
    !
    select type (obj)
    type is (TimeArraySeriesType)
      res => obj
    end select
    !
    return
  end function CastAsTimeArraySeriesType

  function GetTimeArraySeriesFromList(list, indx) result (res)
! ******************************************************************************
! GetTimeArraySeriesFromList -- get time array from list
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    type(ListType),          intent(inout) :: list
    integer,                 intent(in)    :: indx
    type(TimeArraySeriesType), pointer :: res
    ! -- local
    class(*), pointer :: obj
! ------------------------------------------------------------------------------
    !
    obj => list%GetItem(indx)
    res => CastAsTimeArraySeriesType(obj)
    !
    return
  end function GetTimeArraySeriesFromList

end module TimeArraySeriesModule


text:
module readaerobs
! module: readaerobs
!
!abstract: read aerosol data from  files.
! AOD data
!
!program history log:
! 2015-10 Buchard Initial version

use kinds, only: r_double,r_single,i_kind,r_kind
use params, only: expid,nhr_anal,nbackgrounds
use m_odsmeta, only: ktAOD,ktLogAOD
use m_tick, only: tick
use m_ods
use m_die
use Chem_RegistryMod
use Chem_BundleMod

integer,parameter :: my_nobs=-1
logical,parameter :: showjob=.false.

contains

subroutine get_num_aerobs(obspath,datestring,num_obs_tot,id)
    implicit none
    character (len=120), intent(in) :: obspath
    character (len=10), intent(in) :: datestring
    character(len=8), intent(in) :: id
    integer(i_kind), intent(out) :: num_obs_tot
    character(len=120) obsfile
    character(len=20) :: ftype

    character(len=*), parameter:: myname='get_num_aerobs'

    character(10) :: obstype
    character(12) :: geosdate
    character(14) :: geosdate2
    integer(i_kind) iunit, nhms, nymd, nymd0, nhms0, nsel, nb, rc

    type(ods_vect):: ods

! For now in the ods file, only observations data (AOD) after buddy check are available.

   read(datestring,'(i8,i2)') nymd0, nhms0
   nhms0 = nhms0 * 10000
   call tick (nymd0,nhms0,-6*3600)

   num_obs_tot=0; nb=1
   do while (nhr_anal(nb) > 0)
      nymd=nymd0;nhms=nhms0
      call tick (nymd,nhms,nhr_anal(nb)*3600)
      write(geosdate ,'(i8.8,a,i2.2,a)') nymd, '_', nhms/10000, 'z'
      write(geosdate2,'(i8.8,a,i4.4,a)') nymd, '_', nhms/100, 'z'

      obsfile = trim(adjustl(obspath))//"ensmean/"//trim(expid)//".ensmean_aod.obs."//geosdate//".ods"
      print *,'[r] reading ensmean aerosol file: ', trim(obsfile)

      call ods_get (obsfile, nymd, nhms, ftype, ods, rc)
      if ( rc /= 0 ) then
           num_obs_tot = 0
      end if

      call ODS_Select ( ods, ods%data%nobs, nsel, rc, &
                        kt=ktAOD, qcexcl=0 )

      if (ANY(ods%data%kt /= ktAOD)) then
         call die(myname, 'expecting only log(AOD) in observation files')
      endif
      if (my_nobs>0) then
        num_obs_tot = my_nobs
      else
        num_obs_tot = num_obs_tot + nsel
     endif

     call ODS_Clean(ods,rc)

     nb=nb+1
   end do
   write(6,'(a,2(1x,i7))') 'AOD observation count (total,qc-cleared): ', ods%data%nobs, num_obs_tot

end subroutine get_num_aerobs

subroutine get_aerobs_data(obspath, datestring, nobs_max,h_x_ensmean, h_xnobc, x_obs, x_err, &
           x_lon, x_lat, x_press, x_time, x_code, x_errorig, x_type,id, id2)

  implicit none
  character(len=120), intent(in) :: obspath
  character(len=10), intent(in)  :: datestring
  character(len=10), intent(in)  :: id2   ! member number format mem001 for ex
  character(len=10), intent(in)  :: id    ! m

  integer, intent(in) :: nobs_max

  real(r_single), dimension(nobs_max), intent(out) :: h_x_ensmean
  real(r_single), dimension(nobs_max), intent(out) :: h_xnobc
  real(r_single), dimension(nobs_max), intent(out) :: x_obs
  real(r_single), dimension(nobs_max), intent(out) :: x_err
  real(r_single), dimension(nobs_max), intent(out) :: x_lon
  real(r_single), dimension(nobs_max), intent(out) :: x_lat
  real(r_single), dimension(nobs_max), intent(out) :: x_press
  real(r_single), dimension(nobs_max), intent(out) :: x_time
  real(r_single), dimension(nobs_max), intent(out) :: x_errorig

  integer(i_kind), dimension(nobs_max), intent(out) :: x_code
  character(len=20), dimension(nobs_max), intent(out) ::  x_type

  character(len=*), parameter:: myname='get_aerobs_data'

  type(ods_vect):: mods ! mean
  type(ods_vect):: eods ! member e
  type(ods_vect),allocatable:: odsm(:)
  type(ods_vect),allocatable:: odse(:)
  type(ods_vect):: odsx         ! auxiliar

  logical :: twofiles

  real(r_single), allocatable, dimension(:) :: press, errorig
  real,parameter :: eps=0.01
  integer, allocatable,dimension(:) :: code
  character, allocatable, dimension(:) :: obstype
  integer :: nhms, nymd, nymd0, nhms0, rc, nobs_aer, nsel, nb

  character(len=120) :: fname
  character(len=120) :: obsfile
  character(len=12)  :: geosdate
  character(len=14)  :: geosdate2
  character(len=3) :: obtype
  character(len=20) :: ftype

  real,allocatable,dimension(:):: job

  twofiles = id /= id2

  allocate(odsm(nbackgrounds))
  allocate(odse(nbackgrounds))

  read(datestring,'(i8,i2)') nymd0, nhms0
  nhms0 = nhms0 * 10000
  call tick (nymd0,nhms0,-6*3600)

  nobs_aer=0; nb=1
  do while (nhr_anal(nb) > 0)
     nymd=nymd0;nhms=nhms0
     call tick (nymd,nhms,nhr_anal(nb)*3600)
     write(geosdate ,'(i8.8,a,i2.2,a)') nymd, '_', nhms/10000, 'z'
     write(geosdate2,'(i8.8,a,i4.4,a)') nymd, '_', nhms/100, 'z'

!    Read ods file containing AOD observations after buddy check and omf for each member and ensmean
!    -----------------------------------------------------------------------------------------------
     obsfile = trim(adjustl(obspath))//"ensmean/"//trim(expid)//".ensmean_aod.obs."//geosdate//".ods"
     if(trim(id2)=="mem001") print *,'[r] reading omf aerosol file:',id, obsfile

     call ods_get (obsfile, nymd, nhms, ftype, odsx, rc)
     if ( rc /= 0 ) then
            call die(myname, 'could not read ODS aerosol files')
     end if
     call ODS_Select ( odsx, odsx%data%nobs, nsel, rc, &
                       odss=odsm(nb), &
                       kt=ktAOD, qcexcl=0 )
!    odsm(nb)%data%time = nhms/10000!+odsm(nb)%data%time/60
     call ods_clean (odsx,rc)
     if (ANY(odsm(nb)%data%kt /= ktAOD)) then
        call die(myname, 'expecting only AOD in observation files')
     endif

     if(twofiles) then

       obsfile = trim(adjustl(obspath))//trim(id2)//"/"//"aer_omf.ext_Nc."//geosdate//".ods"
       print *,'[r] reading omf aerosol file twofile:', id2, obsfile

       call ods_get (obsfile, nymd, nhms, ftype, odsx, rc)
       if ( rc /= 0 ) then
          call die(myname, 'could not read ODS aerosol files')
       end if
       call ODS_Select ( odsx, odsx%data%nobs, nsel, rc, &
                         odss=odse(nb), &
                         kt=ktAOD, qcexcl=0 )
!      odse(nb)%data%time = nhms/10000!+odse(nb)%data%time/60
       call ods_clean (odsx,rc)
       if (ANY(odse(nb)%data%kt /= ktAOD)) then
          call die(myname, 'expecting only log(AOD) in observation files')
       endif
     endif

     ! Consistency check
     ! ------------------
     if(twofiles) then
       if (odsm(nb)%data%nobs /= odse(nb)%data%nobs) then
          print*, '[w], warning: the 2 ods files do not have same number of obs:', odsm(nb)%data%nobs, odse(nb)%data%nobs
          call die(myname, 'Aborting ...')
       endif
       if (ANY(odsm(nb)%data%qcexcl/=odse(nb)%data%qcexcl) ) then
          call die(myname, ': ensmean and members have inconsistent QC excl marks, aborting ...')
       endif
     endif

     nobs_aer = nobs_aer + odsm(nb)%data%nobs
     if(my_nobs>0) nobs_aer=my_nobs
     nb = nb + 1
  end do ! nb

  if ( nobs_aer /= nobs_max ) then
      call die(myname, 'inconsistent aerosol obs count, aborting ...')
  endif

  call ODS_Merge ( odsm, nbackgrounds, mods, rc )
  do nb=nbackgrounds,1,-1
     call ODS_Clean ( odsm(nb), rc )
  enddo
  deallocate(odsm)
  if(twofiles) then
     call ODS_Merge ( odse, nbackgrounds, eods, rc )
     do nb=nbackgrounds,1,-1
        call ODS_Clean ( odse(nb), rc )
     enddo
     deallocate(odse)
  endif

  allocate(press(nobs_aer), errorig(nobs_aer), obstype(nobs_aer), code(nobs_aer))
  code   (1:nobs_aer) = 1_i_kind
  press  (1:nobs_aer) = 1.e15_r_single ! total column quantities require oblev to be undef
  errorig(1:nobs_aer) = 0.20_r_single

! These are log of AOD ...
  h_x_ensmean = log(max(eps,eps+mods%data%obs-mods%data%omf))
  x_obs       = log(max(eps,eps+mods%data%obs))
  if(twofiles) then
    h_xnobc = log(max(eps,eps+eods%data%obs-eods%data%omf))
    if (showjob) then
       allocate(job(nobs_aer))
       job = ( log(max(eps,eps+eods%data%obs))-log(max(eps,eps+eods%data%obs-eods%data%omf)) )/errorig(1:nobs_aer)
       job = job*job
       print *, 'Member Jo(omb)/n = ', sum(job)/nobs_aer
       deallocate(job)
    endif
  else
    h_xnobc = h_x_ensmean
  endif

  ! output needed
  ! --------
  x_type = 'aod'
  x_code = code
  x_lat  = mods%data%lat(1:nobs_aer)
  x_lon  = mods%data%lon(1:nobs_aer)
  x_time = mods%data%time(1:nobs_aer)/60 ! hours relative to analysis time
  x_press= press
!  if (( mods%data%kx == 301).AND.(mods%data%kx == 302 )) then      ! modo modl   ! is there ods%data%index ????
!       x_errorig(i:) = 0.18
!  else if (( mods%data%kx == 311).AND.(mods%data%kx == 312 )) then ! mydo mydl
!       x_errorig(:) = 0.20
!  else if ( mods%data%kx == 313 ) then ! misr
!       x_errorig(:) = 0.15
!  else if ( mods%data%kx == 323 ) then ! aeronet
!       x_errorig(:) = 0.15
!  endif
  x_errorig = errorig*errorig
  x_err     = errorig*errorig

!  301  MODO  MODIS/TERRA-Ocean
!  302  MODL  MODIS/TERRA-Land
!  311  MYDO  MODIS/AQUA-Ocean
!  312  MYDL  MODIS/AQUA-Land
!  313  MISR  MISR
!  323  ANET  AERONET

  if (showjob) then
     allocate(job(nobs_aer))
     job = ( log(max(eps,eps+mods%data%obs))-log(max(eps,eps+mods%data%obs-mods%data%omf)) )/errorig(1:nobs_aer)
     job = job*job
     print *, 'Mean Jo(omb)/n = ', sum(job)/nobs_aer
     deallocate(job)
  endif

  call ODS_Clean(mods,rc)
  if(twofiles) then
     call ODS_Clean(eods,rc)
  endif
  deallocate(obstype, code, press, errorig)

end subroutine get_aerobs_data

end module readaerobs

max_stars_repo_name: GEOS-ESM/GEOSadas

max_issues_repo_name: GEOS-ESM/GEOSadas

max_forks_repo_name: GEOS-ESM/GEOSadas


text:
great_rsd(e1,f1).
great_rsd(t1,m1).
great_rsd(u1,t1).

ext: f

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


!
! Copyright 2005-2017 ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
!
! Description: how to get values using keys.
!
program read_message
use grib_api
implicit none
  integer                      ::  ifile,ofile
  integer                      ::  iret,igrib
  integer , dimension(50000)   :: buffer
  integer(kind=kindOfSize_t)   :: len1
  integer                      :: step,level

  call grib_open_file(ifile,'../../data/index.grib','r')
  call grib_open_file(ofile,'out.readmsg.grib','w')

! a grib message is read from file into buffer
  len1=size(buffer)*4
  call  grib_read_from_file(ifile,buffer,len1,iret)

  do while (iret/=GRIB_END_OF_FILE)

!   a new grib message is created from buffer
    call grib_new_from_message(igrib,buffer)

!   get as a integer
    call grib_get(igrib,'step', step)
    write(*,*) 'step=',step

    call grib_get(igrib,'level',level)
    write(*,*) 'level=',level

    call grib_release(igrib)

    call grib_write_bytes(ofile,buffer,len1)

!   a grib message is read from file into buffer
    len1=size(buffer)*4
    call  grib_read_from_file(ifile,buffer,len1,iret)

  enddo

  call grib_close_file(ifile)
  call grib_close_file(ofile)

end program read_message


text:
c
c =====================================
      subroutine chol_inv(a,x,b,n,ierr)
c =====================================

c   solve a.x = b using Choleski decomposition
c
c      real*8 a(n,n),x(n),b(n),p(100)     ! n <= 100
      real*4 a(n,n),x(n),b(n),p(100)     ! n <= 100
c
c  Choleski Decomposition  L. LT = A
c            L is returned in the lower triangle of A, except diagonal
c            elements   stored in p.
c
      ierr=0
      do i=1,n
         do j=1,n
            sum=a(i,j)
            do k=i-1,1,-1
               sum = sum - a(i,k)*a(j,k)
c!!!!	       print *,'SUM aik ajk',sum,a(i,k),a(j,k),i,j,k
            enddo
            if(i.eq.j) then
              if(sum.le.0.0) then
                ierr=-1
                print *,'  Choleski Decomposition Failed'
                go to 1
              else
                p(i)=sqrt(sum)
              endif
            else
              a(j,i) = sum/p(i)
            endif
         enddo
      enddo

   1  if(ierr.eq.0) then     !solve the linear system using forward
                             !and back substitution
c    Solve A x = L (LT x) = b
c                             First solve L y = b
c                             Then  solve LT x = y
c
      do i=1,n         ! solve L y = b, storing y in x
         sum=b(i)
         do k=i-1,1,-1
            sum = sum - a(i,k)*x(k)
         enddo
         x(i) = sum/p(i)
      enddo
c
      do i=n,1,-1      ! solve LT x = y
         sum=x(i)
         do k=i+1,n
            sum = sum - a(k,i)*x(k)
         enddo
         x(i) = sum/p(i)
      enddo
c
      endif

      end
c

max_stars_repo_name: mmbell/Airborne-Radar-QC

max_issues_repo_name: DalavanCloud/Airborne-Radar-QC

max_forks_repo_name: DalavanCloud/Airborne-Radar-QC


text:
      SUBROUTINE DISEDB (SAMP, START, UNIV, SUBSET)

C  DISEDB selects SAMP predefined bundles at random from the
C  edge bundle table which are distinct in all their attributes.
C  If there are not SAMP completely distinct bundles, DISEDB tries
C  to maximize the number of different attributes.
C
C  Input parameters:
C    SAMP        : number of predefined bundles to be picked
C    START       : starting index
C    UNIV        : last element in universe to pick from
C  Output parameters:
C    SUBSET      : array containing the selected distinct bundles

      COMMON /DIALOG/ DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT,
     2                DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS
      INTEGER         DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT, ERRIND
      REAL            DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS

C off/on switch for edge flag and error handling mode
      INTEGER    POFF,     PON
      PARAMETER (POFF = 0, PON = 1)

C  total number of attributes
      INTEGER    TOTATR
      PARAMETER (TOTATR=4)

      INTEGER    SAMP, START, UNIV, SUBSET(SAMP), I,J,JJ,K, OT
      INTEGER    ATRCNT, MAXATR,MAXLOC, JEF, JET,  JECI,
     1           KEF, KET, KECI, RNDINT, OPTATR

      REAL       JEW,KEW

      LOGICAL    ATRNEW(TOTATR)

      IF (UNIV+1-START .LT. SAMP) THEN
         CALL UNMSG ('Abort in DISEDB because size of universe ' //
     1               'is less than requested sample.')
      ENDIF

C  OPTATR is the best possible number of attributes that can be distinct
      OPTATR = TOTATR

C  this loop picks out samp elements
      DO 500 I = 1, SAMP
C  number of new attribute values so far
         MAXATR = -1
C  look at entire universe each time, starting at random location
         J = RNDINT(START, UNIV)
         DO 400 JJ = START, UNIV
            IF (J.GE.UNIV) THEN
C  cycle around to beginning
               J = START
            ELSE
               J = J+1
            ENDIF
            DO 50 OT = 1,TOTATR
               ATRNEW(OT) = .TRUE.
50          CONTINUE
C  get attributes of next candidate
            CALL PQPEDR (SPECWT, J, ERRIND, JEF, JET, JEW, JECI)
            CALL CHKINQ ('pqpedr', ERRIND)

C  check against all those picked so far;
            DO 300 K = 1, I-1
C  if already picked, get another j
               IF (J .EQ. SUBSET(K)) GOTO 400
               CALL PQPEDR (SPECWT, SUBSET(K),ERRIND,KEF,KET,KEW,KECI)
               CALL CHKINQ ('pqpedr', ERRIND)
C  check all attributes
               IF (JEF .EQ. KEF) THEN
                  IF (JEF .EQ. POFF) THEN
C  edge flag off implies all other attributes indistinguishable
                     ATRCNT = 0
                     GOTO 80
                  ELSE
                     ATRNEW(1) = .FALSE.
                  ENDIF
               ENDIF
               IF (JET .EQ. KET) ATRNEW(2) = .FALSE.
               IF (JEW .EQ. KEW) ATRNEW(3) = .FALSE.
               IF (JECI .EQ. KECI) ATRNEW(4) = .FALSE.
300         CONTINUE
C  count # of trues = # of distinct attributes
            ATRCNT = 0
            DO 75 OT = 1,TOTATR
               IF (ATRNEW(OT)) ATRCNT = ATRCNT + 1
75          CONTINUE
80          CONTINUE
C  take the best so far
            IF (ATRCNT .GT. MAXATR) THEN
               MAXATR = ATRCNT
               MAXLOC = J
            ENDIF
C  cannot get better than OPTATR
            IF (MAXATR .GE. OPTATR) GOTO 410
400      CONTINUE

C  put best one on the list
410      CONTINUE
         SUBSET(I) = MAXLOC
C  remember greatest number of distinct attributes
         OPTATR = MAXATR
500   CONTINUE

      END


text:
*DECK QK21
      SUBROUTINE QK21 (F, A, B, RESULT, ABSERR, RESABS, RESASC)
C***BEGIN PROLOGUE  QK21
C***PURPOSE  To compute I = Integral of F over (A,B), with error
C                           estimate
C                       J = Integral of ABS(F) over (A,B)
C***LIBRARY   SLATEC (QUADPACK)
C***CATEGORY  H2A1A2
C***TYPE      SINGLE PRECISION (QK21-S, DQK21-D)
C***KEYWORDS  21-POINT GAUSS-KRONROD RULES, QUADPACK, QUADRATURE
C***AUTHOR  Piessens, Robert
C             Applied Mathematics and Programming Division
C             K. U. Leuven
C           de Doncker, Elise
C             Applied Mathematics and Programming Division
C             K. U. Leuven
C***DESCRIPTION
C
C           Integration rules
C           Standard fortran subroutine
C           Real version
C
C           PARAMETERS
C            ON ENTRY
C              F      - Real
C                       Function subprogram defining the integrand
C                       FUNCTION F(X). The actual name for F needs to be
C                       Declared E X T E R N A L in the driver program.
C
C              A      - Real
C                       Lower limit of integration
C
C              B      - Real
C                       Upper limit of integration
C
C            ON RETURN
C              RESULT - Real
C                       Approximation to the integral I
C                       RESULT is computed by applying the 21-POINT
C                       KRONROD RULE (RESK) obtained by optimal addition
C                       of abscissae to the 10-POINT GAUSS RULE (RESG).
C
C              ABSERR - Real
C                       Estimate of the modulus of the absolute error,
C                       which should not exceed ABS(I-RESULT)
C
C              RESABS - Real
C                       Approximation to the integral J
C
C              RESASC - Real
C                       Approximation to the integral of ABS(F-I/(B-A))
C                       over (A,B)
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  R1MACH
C***REVISION HISTORY  (YYMMDD)
C   800101  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C***END PROLOGUE  QK21
C
      REAL A,ABSC,ABSERR,B,CENTR,DHLGTH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,
     1  FV1,FV2,HLGTH,RESABS,RESG,RESK,RESKH,RESULT,R1MACH,UFLOW,WG,WGK,
     2  XGK
      INTEGER J,JTW,JTWM1
      EXTERNAL F
C
      DIMENSION FV1(10),FV2(10),WG(5),WGK(11),XGK(11)
C
C           THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1).
C           BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
C           CORRESPONDING WEIGHTS ARE GIVEN.
C
C           XGK    - ABSCISSAE OF THE 21-POINT KRONROD RULE
C                    XGK(2), XGK(4), ...  ABSCISSAE OF THE 10-POINT
C                    GAUSS RULE
C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
C                    ADDED TO THE 10-POINT GAUSS RULE
C
C           WGK    - WEIGHTS OF THE 21-POINT KRONROD RULE
C
C           WG     - WEIGHTS OF THE 10-POINT GAUSS RULE
C
      SAVE XGK, WGK, WG
      DATA XGK(1),XGK(2),XGK(3),XGK(4),XGK(5),XGK(6),XGK(7),
     1  XGK(8),XGK(9),XGK(10),XGK(11)/
     2         0.9956571630258081E+00,     0.9739065285171717E+00,
     3     0.9301574913557082E+00,     0.8650633666889845E+00,
     4     0.7808177265864169E+00,     0.6794095682990244E+00,
     5     0.5627571346686047E+00,     0.4333953941292472E+00,
     6     0.2943928627014602E+00,     0.1488743389816312E+00,
     7     0.0000000000000000E+00/
C
      DATA WGK(1),WGK(2),WGK(3),WGK(4),WGK(5),WGK(6),WGK(7),
     1  WGK(8),WGK(9),WGK(10),WGK(11)/
     2     0.1169463886737187E-01,     0.3255816230796473E-01,
     3     0.5475589657435200E-01,     0.7503967481091995E-01,
     4     0.9312545458369761E-01,     0.1093871588022976E+00,
     5     0.1234919762620659E+00,     0.1347092173114733E+00,
     6     0.1427759385770601E+00,     0.1477391049013385E+00,
     7     0.1494455540029169E+00/
C
      DATA WG(1),WG(2),WG(3),WG(4),WG(5)/
     1     0.6667134430868814E-01,     0.1494513491505806E+00,
     2     0.2190863625159820E+00,     0.2692667193099964E+00,
     3     0.2955242247147529E+00/
C
C
C           LIST OF MAJOR VARIABLES
C           -----------------------
C
C           CENTR  - MID POINT OF THE INTERVAL
C           HLGTH  - HALF-LENGTH OF THE INTERVAL
C           ABSC   - ABSCISSA
C           FVAL*  - FUNCTION VALUE
C           RESG   - RESULT OF THE 10-POINT GAUSS FORMULA
C           RESK   - RESULT OF THE 21-POINT KRONROD FORMULA
C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B),
C                    I.E. TO I/(B-A)
C
C
C           MACHINE DEPENDENT CONSTANTS
C           ---------------------------
C
C           EPMACH IS THE LARGEST RELATIVE SPACING.
C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
C
C***FIRST EXECUTABLE STATEMENT  QK21
      EPMACH = R1MACH(4)
      UFLOW = R1MACH(1)
C
      CENTR = 0.5E+00*(A+B)
      HLGTH = 0.5E+00*(B-A)
      DHLGTH = ABS(HLGTH)
C
C           COMPUTE THE 21-POINT KRONROD APPROXIMATION TO
C           THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
C
      RESG = 0.0E+00
      FC = F(CENTR)
      RESK = WGK(11)*FC
      RESABS = ABS(RESK)
      DO 10 J=1,5
        JTW = 2*J
        ABSC = HLGTH*XGK(JTW)
        FVAL1 = F(CENTR-ABSC)
        FVAL2 = F(CENTR+ABSC)
        FV1(JTW) = FVAL1
        FV2(JTW) = FVAL2
        FSUM = FVAL1+FVAL2
        RESG = RESG+WG(J)*FSUM
        RESK = RESK+WGK(JTW)*FSUM
        RESABS = RESABS+WGK(JTW)*(ABS(FVAL1)+ABS(FVAL2))
   10 CONTINUE
      DO 15 J = 1,5
        JTWM1 = 2*J-1
        ABSC = HLGTH*XGK(JTWM1)
        FVAL1 = F(CENTR-ABSC)
        FVAL2 = F(CENTR+ABSC)
        FV1(JTWM1) = FVAL1
        FV2(JTWM1) = FVAL2
        FSUM = FVAL1+FVAL2
        RESK = RESK+WGK(JTWM1)*FSUM
        RESABS = RESABS+WGK(JTWM1)*(ABS(FVAL1)+ABS(FVAL2))
   15 CONTINUE
      RESKH = RESK*0.5E+00
      RESASC = WGK(11)*ABS(FC-RESKH)
      DO 20 J=1,10
        RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH))
   20 CONTINUE
      RESULT = RESK*HLGTH
      RESABS = RESABS*DHLGTH
      RESASC = RESASC*DHLGTH
      ABSERR = ABS((RESK-RESG)*HLGTH)
      IF(RESASC.NE.0.0E+00.AND.ABSERR.NE.0.0E+00)
     1  ABSERR = RESASC*MIN(0.1E+01,
     2  (0.2E+03*ABSERR/RESASC)**1.5E+00)
      IF(RESABS.GT.UFLOW/(0.5E+02*EPMACH)) ABSERR = MAX
     1  ((EPMACH*0.5E+02)*RESABS,ABSERR)
      RETURN
      END


text:
      subroutine openou(noutpt,nttyo,ufiln,uform,nrecl,ilu)
c
c     This subroutine opens an output type file. If a file of the
c     same name already exists, it is first destroyed. An unused
c     logical unit number is obtained.
c
c     This subroutine is called by:
c
c       EQ3NR/eq3nr.f
c       EQ6/eq6.f
c
c-----------------------------------------------------------------------
c
c     Input:
c
c       noutpt = the unit number of the output file
c       nttyo  = the unit number of the screen file
c       ufiln  = the name of the file to open
c       uform  = the file format, 'formatted' or 'unformatted'
c       nrecl  = the record length (number of characters per line)
c
c     Output:
c
c       ilu     = the logical unit number of opened file
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer noutpt,nttyo
c
      integer ilu,nrecl
c
c     Note- the character length for the variables ufiln and uform
c     can not be exactly specified.
c
      character*(*) ufiln,uform
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer j2,nerr
c
      integer ilnobl
c
      logical qex
c
      character*8 ustat
c
c-----------------------------------------------------------------------
c
      j2 = ilnobl(ufiln)
c
c     See if a file of the same name already exists. If so,
c     destroy it. This makes the logical unit number available.
c     If a file of the same name does not exist, get the next
c     available logical unit number.
c
      inquire(file=ufiln,exist=qex)
      if (qex) then
        ustat = 'old'
        call getlu(ilu,nerr)
        if (nerr .ne. 0) then
          if (noutpt .gt. 0) write (noutpt,1000) ustat,ufiln(1:j2)
          write (nttyo,1000) ustat,ufiln(1:j2)
 1000     format(/' * Error - (EQLIBU/openou) No logical unit number',
     $    /7x,'is available to open the ',a3,' file "',a,'".')
          stop
        endif
        open(ilu,file=ufiln,status=ustat,err=10)
        close(ilu,status='delete',err=15)
      else
        ustat = 'new'
        call getlu(ilu,nerr)
        if (nerr .ne. 0) then
          if (noutpt .gt. 0) write (noutpt,1000) ustat,ufiln(1:j2)
          write (nttyo,1000) ustat,ufiln(1:j2)
          stop
        endif
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Open the new file.
c
      if (nrecl .gt. 0) then
c
c       Use the specified record length.
c
        open(ilu,file=ufiln,form=uform,status='new',recl=nrecl,err=10)
      else
c
c       The record length is not specified. Open the file at the
c       default record length.
c
        open(ilu,file=ufiln,form=uform,status='new',err=10)
      endif
      go to 999
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
   10 if (noutpt .gt. 0) write (noutpt,1010) ustat,ufiln(1:j2)
      write (nttyo ,1010) ustat,ufiln(1:j2)
 1010 format(/" * Error - (EQLIBU/openou) Can't open the ",a3,' copy',
     $ /7x,'of the file "',a,'".')
      stop
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
   15 if (noutpt .gt. 0) write (noutpt,1020) ufiln(1:j2)
      write (nttyo,1020) ufiln(1:j2)
 1020 format(/" * Error - (EQLIBU/openou) Can't delete the old copy",
     $ /7x,'of the file "',a,'".')
      stop
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
  999 continue
      end


use IOClass
implicit none

real(real64),allocatable :: a(:,:)
type(IO_) :: f

a = f%import("test.csv",num_column=2)

a(1,1) = 2.0d0
print *, a(:,1)

call f%export("test2.csv",a)

end


text:
      subroutine remove_user_header_param(iuhead, nuhead, kuhead,
     *           cparam)
      !
      ! Removes a parameter and its associated value from the user header
      !
      ! iuhead  The user header integer array
      ! nuhead  The number of integers used in the user header (size)
      ! kuhead  The number of integers in the user header (capacity)
      ! cparam  The parameter to retrieve the value for
      !
      ! Mike Perryman
      ! USACE Hydrologic Engineering Center
      ! Nov 2021
      !
      implicit none
      !-------------------!
      ! define parameters !
      !-------------------!
      integer   iuhead(*), nuhead, kuhead
      character cparam*(*)
      !-----------------!
      ! local variables !
      !-----------------!
      integer   iuhead_copy(500), max_head_len, ifirst, ilast, i, plen
      character cuhead*2000
      logical   param_has_separator

      equivalence (iuhead_copy, cuhead)

      include 'zdssmz.h'

      !-------------------------------------------------------------------!
      ! copy the user header to a local variable so we can equivalence it !
      !-------------------------------------------------------------------!
      iuhead_copy = 0
      if (kuhead.gt.size(iuhead_copy)) then
        if (mlevel.ge.1) then
          write (munit,'(/a,i0,a,/,a,i0,a,/,a)')
     *    ' WARNING: Cannot modify user header >  ',
     *    size(iuhead_copy) * 4,
     *    ' characters long',
     *    ' User header is ',
     *    kuhead * 4,
     *    ' characters long',
     *    ' User header not modified'
        end if
        return
      end if
      max_head_len = min(kuhead, size(iuhead_copy))
      iuhead_copy(:max_head_len) = iuhead(:max_head_len)
      !--------------------------------------------------------!
      ! blank everything past the copy size or null terminator !
      !--------------------------------------------------------!
      ifirst = index(cuhead, char(0))
      if (ifirst.gt.0) then
        cuhead(ifirst:) = ' '
      else
        if (max_head_len.lt.size(iuhead_copy)) then
          ifirst = max_head_len*4
          cuhead(ifirst:) = ' '
        end if
      end if
      !-----------------------------------------!
      ! locate the parameter name in the header !
      !-----------------------------------------!
      plen = len_trim(cparam)
      ifirst = index(cuhead, cparam(1:plen)//':')
      !-----------------------------------------------!
      ! remove the parameter and its associated value !
      !-----------------------------------------------!
      if (ifirst.gt.0) then
        ilast = index(cuhead(ifirst:), ';')
        if (ilast.eq.0) then
          !------------------!
          ! at end of header !
          !------------------!
          cuhead(ifirst:ifirst) = ':'
          cuhead(ifirst+1:) = ' '
        else
          !-------------------------!
          ! in the middle of header !
          !-------------------------!
          ilast = ilast + ifirst - 1
          cuhead(ifirst:) = cuhead(ilast+1:)
        end if
      end if
      !-------------------------------!
      ! remove any leading delimiters !
      !-------------------------------!
      if (cuhead(1:1).eq.';') cuhead = cuhead(2:)
      iuhead(:kuhead) = 0
      iuhead(:max_head_len) = iuhead_copy(:max_head_len)
      nuhead = (len_trim(cuhead)-1)/4+1
      return
      end

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


text:
! FILE: sort.f
! PURPOSE: demonstrate the use of "qsort_inline.inc" and
! "qsort_inline_index.inc". These can be used as specific
! sort procedures under a common SORT generic name.
!---------------------------------------------------------------
! Sort a string array, with any string length.
subroutine sortp_string(array_size,index,string)
  integer, parameter :: QSORT_THRESHOLD = 32
  integer, intent(in) :: array_size
  integer, intent(out) :: index(array_size)
  character(len = *), intent(in) :: string(array_size)
  include "qsort_inline.inc"
contains
  include "qsort_inline_index.inc"
  logical &
  function less_than(a,b)
    integer, intent(in) :: a,b
    if ( string(index(a)) == string(index(b))  ) then
      less_than = ( index(a) < index(b) )
    else
      less_than = ( string(index(a)) < string(index(b)) )
    end if
  end function less_than
end subroutine sortp_string
!---------------------------------------------------------------
! Sort an index-array by its double array target
subroutine sortp_1r8(array_size,index,value)
  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: QSORT_THRESHOLD = 32
  integer, intent(in) :: array_size
  integer, intent(inout) :: index(array_size)
  real(dp), intent(in) :: value(array_size)
  include "qsort_inline.inc"
contains
  include "qsort_inline_index.inc"
  logical &
  function less_than(a,b)
    integer, intent(in) :: a,b
    real(dp), parameter :: small = 1.0e-15
    if ( abs(value(index(a))-value(index(b))) < small ) then
      less_than = index(a) < index(b)
    else
      less_than = value(index(a)) < value(index(b))
    end if
  end function less_than
end subroutine sortp_1r8
!---------------------------------------------------------------
! Sort a single-precision real array by index, with a fuzzy equality test
subroutine sortp_1r4(array_size,index,value)
  integer, parameter :: QSORT_THRESHOLD = 32
  integer, intent(in) :: array_size
  integer, intent(inout) :: index(array_size)
  real(4), intent(in) :: value(array_size)
  include "qsort_inline.inc"
contains
  include "qsort_inline_index.inc"
  logical &
  function less_than(a,b)
    integer, intent(in) :: a,b
    real(4), parameter :: small = 1.0e-6
    if ( abs(value(index(a))-value(index(b))) < small ) then
      less_than = index(a) < index(b)
    else
      less_than = value(index(a)) < value(index(b))
    end if
  end function less_than
end subroutine sortp_1r4
!---------------------------------------------------------------
! Sort an array of integers
subroutine sort_1i(array_size,i1)
  integer, parameter :: QSORT_THRESHOLD = 32
  integer, intent(in) :: array_size
  integer, intent(inout), dimension(array_size) :: i1
  include "qsort_inline.inc"
contains
  subroutine init()
  end subroutine init
  logical &
  function less_than(a,b)
    integer, intent(in) :: a,b
    if ( i1(a) == i1(b) ) then
      less_than = a < b
    else
      less_than = i1(a) < i1(b)
    end if
  end function less_than
  subroutine swap(a,b)
    integer, intent(in) :: a,b
    integer :: hold
    hold = i1(a); i1(a) = i1(b); i1(b) = hold
  end subroutine swap
! circular shift-right by one:
  subroutine rshift(left,right)
    integer, intent(in) :: left, right
    integer :: hold
    hold = i1(right); i1(left + 1:right) = i1(left:right-1); i1(left) = hold
  end subroutine rshift
end subroutine sort_1i
!---------------------------------------------------------------

max_stars_repo_name: SpakowitzLab/BasicWLC

max_issues_repo_name: riscalab/wlcsim

max_forks_repo_name: riscalab/wlcsim


text:
Located just downstream of downtown Winters, the Winters Swimming Hole found both Winters locals and adventurous Davisites enjoying its cool waters.  People can be found splashing into the water as they jump from a wellused rope swing.  People were also once seen sun bathing on the defunct concrete dam.

On October 10th, 2009 the dam was dismantled and the river was allowed to flow freely once again.

Directions
By Bicycles bicycle  Take the Russell Boulevard bike path to its completion, turn left on Stevenson Bridge Rd, right on Putah Creek road (7mi), and turn right over the pedestrian bridge next to Winters Rd.  Once over the bridge, turn right and walk your bike a few hundred yards down the pedestrian footpath until you see the rope swing hanging from a large Sycamore.  The trip turns out to be 14.25 miles total from campus to the rope swing.

Driving: Take Covell all the way to Winters.  Once in down town winters make a left past Buckhorn and you will see a Gazebo.  Park there/park your bike there).  Head past the Gazebo and there should be a trail.  Walk down the trail about 300 feet by Putah Creek and you will see the rope swing to your right.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
Address(Santa Rosa Street) is a residential street in West Davis.

 Intersecting Streets

  Hampton Drive
  Van Damme Drive
  Crystal Grove Lane


text:
	subroutine euler(angles,v)

*	Rotates vector given wrt (i,j,k) through Euler angles (i,w,Omega).

	implicit real*8 (a-h,o-z)
	real*8 v(3),angles(3),B(3,3),temp(3)

	cosi=cos(angles(1))
	sini=sin(angles(1))
	cosw=cos(angles(2))
	sinw=sin(angles(2))
	cosOm=cos(angles(3))
	sinOm=sin(angles(3))

        B(1,1) = cosw*cosOm - sinw*cosi*sinOm
        B(1,2) = cosw*sinOm + sinw*cosi*cosOm
        B(1,3) = sinw*sini
        B(2,1) = -sinw*cosOm - cosw*cosi*sinOm
        B(2,2) = -sinw*sinOm + cosw*cosi*cosOm
        B(2,3) = cosw*sini
        B(3,1) = sini*sinOm
        B(3,2) = -sini*cosOm
        B(3,3) = cosi

	do i=1,3
	   sum=0
	   sum2=0
	   do j=1,3
	      sum=sum + B(j,i)*v(j)
	   enddo
	   temp(i)=sum
	enddo
	do i=1,3
	   v(i)=temp(i)
	enddo

	end

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
 block data bd_printcommon
      implicit none
      integer depth
      common /printcommon/ depth
      data depth /1/
      end

program foobar
   print *,1
end program

max_stars_repo_name: raramakr/aomp

max_issues_repo_name: raramakr/aomp

max_forks_repo_name: raramakr/aomp


text:
Davis local obligatory Carls Jr. also happens to be paired with a Green Burrito, which is common practice in California and Arizona. According to http://news.yahoo.com Yahoo News, system wide there are 1,017 Carls Jr. fast food restaurants, 223 of which are dualbranded with Green Burrito. One of the great things about combining Carls Jr. with a Green Burrito is the ability to put salsa in your burgers or on your fries. They offer a 10% Senior Citizens senior Discounts and Freebies discount.  You can get a Hamburgers burger there, among http://www.carlsjr.com/menu/ other things.

In early 2014 a South Davis location was opened, in the husk of the former Wendys, which includes the long desired drive thru that the East Davis location lacks.

At least one student in 2005 believed the Davis franchise managed to be sufficiently worse than other Carls Jr. franchises and should be avoided at all costs. Hes probably graduated by now. And it is possible the management of the Davis franchise has changed in the intervening neardecade.  If said student was a food science & technology major, its also possible he could be the manager of Carls Jr. now. XMen such as Colossus love Carls Jr. https://www.youtube.com/watch?vPEoQoW8zZ_w link

Outside of the Western U.S., Carls Jr. is named Hardees, and has the same menu and logo, although menu items go by different names (the six dollar burger is a thickburger).  Check out the distribution of the two using http://www.fastfoodmaps.com/static.html Fast Food Maps or look at the map below.

I havent been to the Carls Jr. in Davis, but I like their selection of carbsmart burgers.  SS

  Carbsmart is fatstupid.

20050506 20:22:36 nbsp Carls Jr. is known as Hardees on the East Coast. They sell the same things like the $6 Dollar Burger but breakfast menu items are heavy on biscuits rather than burritos.  Users/CandiceBrown

20050506 23:45:20 nbsp Why isnt there one Carls Jr. that can make DECENT fries!?!? Users/PatrickSing

20050613 20:33:00 nbsp Whats up with not having drivethru? Fast food NEEDS drivethru Users/AshleyOrsaba

20051112 00:33:03 nbsp This Carls Jr. takes the cake on cleaniless and food quality. The cooks, Maricela specifically, do a fantastic job of making picture perfect burgers. There used to be a drive through for bikes at one point in time, which was since removed. Reason is speculated to be an employee breaking the window and the repair guy never replacing it. This branch of Carls Jr. is the highest ranking in the district, as well as the only one without a drive through. Users/TarynHornen

20051208 18:00:55 nbsp The 6 dollar burgers are obscene.  Childhood obesity, and early onset diabetes are a consequence of unhealthy food.  You can say to each his own, but unfortunatey we as a society face the burden of cost when an epidemic of diabetes takes place and the healthcare system becomes strained.  I predict a classaction lawsuit against fast food companies just like the tobacco industry faced in our lifetime Users/DudeNude
     Not to start a flame war, but if youre a parent and you buy your kid a $6 burger, then its not Carls Jr.s fault that your kid is fat and has diabetes.  I dont know of many 9 yearolds who go to Carls Jr. and have $6 out of their pocket to pony up for a burger. Users/CaseyBarc

20051211 03:42:43 nbsp I highly recommend the chocolate shakes and burgers from Carls Jr. in the Silo, but the fries are seriously lacking. Users/ThomasFifield

20061008 01:12:41 nbsp I went to the Carls in the Silo ... the food was nasty. nough said. SS

20061214 05:09:09 nbsp I ordered a hamburger there today. I saw the person making the burger lick some extra dressing off her finger and continue to make the burger. I ate it anyway. Users/DanielWorthington

20070512 18:20:24 nbsp I am an avid health freak, however I have searched the menu and they do have some healthier options such as salads. However, use fat free salad dressing or it will defeat the purpose. The Charbroiled BBQ Chicken Sandwich is actually pretty good tasting and doesnt have too much fat. Users/Jedron
  Getting full fat dressing does not defeat the purpose of a salad. The nutrition information you see is for using the whole dressing packet, which is an obscene amount of dressing (if you use half that much dressing, which is still more than a side of dressing at plutos, your salad wll certainly have less fat than the burgers on the menu). Also salads have more micronutrients per calorie than the burgers, regardless of dressing contents. However, if i want a salad im not going to Carls jr. I usually get a famous star, without mayo or cheese and extra lettuce and tomato, if i end up going here. Note they do not inventory the vegetable toppings so they do not cost extra. edit: due to an employees abuse of the not charging for extra veggies policy (I assume they put a salad on their burger) they now charge 20 cents per extra topping. Users/MattHh

20070512 21:31:51 nbsp Im not sure what is more pathetic, that Carls Jr here quotes a review from the Torrance Daily Breeze or that the Torrance Daily Breeze reviewed Carls Jr. Users/NickSchmalenberger
 haha thats pretty funny. They must have run out of places to review :). Users/MattHh

20070515 14:24:42 nbsp the CJs on Covell is pretty good actually. I like how its paired with the Green Burrito menu, so you can have more choices. The one at the Silo is TRASH. not only is the menu lacking, everything is overpriced. Users/DonaldJaye

20070719 17:55:33 nbsp Best hamburger if it is the 6 dollar burger be honest very good i like the western bacon with extra crispy bacon Users/dickjones

20070804 03:36:42 nbsp Love the big burgers actually i think best around besides the squeeze inn in Sacramento that place rocks. Carls is a great place to get a good fast food burger my opinion Users/Brians

20071227 01:31:34 nbsp Yumyums!!! Users/thelonepiper555

20080115 05:38:16 nbsp The night crew at the Carls Jr in the Silo can be pretty shoddy

During one trip there, I received the wrong soda, took it back and got a new one

Then I got the wrong meal, which I also took back after opening the burger and seeing that it was not mine

As i took it back, another girl went up to ask about her food. As it turns out, I was given her order. As the server walks off to get my order, the girl opens her bag, takes out her burger, and stares at the open wrapper for 5 full seconds. Then she seemed to decide, fk it and just left.

i barely contained my laughter Users/fredchen

20080123 21:32:58 nbsp Does the Silos Carls Jr. accept coupons? Users/sfjhffihojh

20080208 10:06:07 nbsp my local place for the classic hamburger! or the classic burrito! chicken stars are the best! Users/davisrox

20080228 20:22:22 nbsp Yes, the Silo one accepts coupons. I went there last night for a buy one get one free teriyaki burger, and it was great! I highly recommend them for their coupon deals.  Users/AnnaF

20080304 18:16:14 nbsp Buy 1 get 1 free
http://www.fatwallet.com/forums/topic.php?catid18&threadid799825 Users/CalvinFan

20080428 20:20:24 nbsp FAT Profit article. Very good read! http://www.fracknaps.com/read/?pSFFYTBAO Users/billybob

20080428 21:59:26 nbsp Wow, that picture up there makes Carls Jr. seem almost holy... Users/SunjeetBaadkar
 It is!

20080430 13:44:11 nbsp Whats up with the sign at the Silo that says: WARNING, this food causes cancer and birth defects. ?

Anyway, the people who have commented saying they like the carbsmart burgers or use the fatfree dressing are idiots. People who fool themselves into believing these things make fastfood healthy deserve to be fat. Please, rid the world of your stupidity by gorging yourself on food and getting all the healthrelated problems that come with diabetes. Users/Loccster

20090504 16:37:42 nbsp I love the burgers at Carls Jr, theyre so yummy. Users/thtly

20090504 20:26:10 nbsp The bourbon burger is gross, the bbq sauce makes it way too sweet. Users/NickSchmalenberger
  I agree. I tried it with a buy1get1free coupon, and it was all I could do to finish even one. I would have given the other one away, but I didnt know anyone that I disliked enough. Users/JoePomidor

20090505 09:45:10 nbsp The one on East Covell Boulevard is only good for one thing, and that is a place for students from Holmes Junior High School to hang out when they ditch class.  I heard they used to chase students out a few years back (makes me feel old thinking back like this).  But I figure theyd rather the business.   Users/MasonMurray

20090505 09:48:10 nbsp every time I eat fast food, I regret it shortly after. Users/RealComputers

20090818 23:14:44 nbsp they wouldnt let me use my fake coupons that everyone was getting from that Lakers game. you suck!

but dang... guacamole bacon six dollar burger...heavenly fattening.mmmmmmmmmMMM! Users/jasonwon

20090818 23:19:00 nbsp good portions! a real mans meal. none of that calorie counting girly crap. Users/bewbs

20100123 17:16:53 nbsp The Carls Jr in the silo has some of the best shakes! Users/elliotgoethe

20100712 14:38:56 nbsp I have visited Carls Jr. multiple times over the last few years and the staff have always been clean, efficient, and considerate.  Food does seem to take longer than at other fast food restaurants.  No worker has ever had difficulty understanding my order, nor have I had difficulty understanding them (and I can only speak English). Users/MattCundy

20100712 15:26:40 nbsp I live close to Carls Jr and bike here regularly. I have always had good experiences with the staff, I especially like the teriyaki burger and breakfast burger. The avocado burger is also very good. Users/NickSchmalenberger

20100713 01:36:20 nbsp I love the Carls Jr in North Davis. Boy do I love this place. What is it? Its just a good value for the burger, so damn tasty, attribution to the charbroiled goodness. The bike ride there is worth it, and definitely helps tone down the feeling of guilt when I engorge myself into one of their delicious Big Hamburgers, or sometimes I really like to reward myself after a week of monotonous lectures with a Teriyaki burger. Have your crisscut fries, because today, today resources are no object!

tl;dr I love this place. Users/davidennis

20100713 15:17:18 nbsp I consider it the best of the national chain hamburger joints in Davis. Users/JimStewart

20110430 18:16:42 nbsp Im very impressed by the six dollar burgers and the skinon fries. We just dont have good fast food burger places like that in Denmark (only McD and Burger King, and you have to go to a Cafe for a decent burger, hence the expression only buy burgers where you can also buy a beer with it). Unfortunately, Ive usually been at Carls jr at campus, where every other burger is cold. Does the north Davis place have that problem too? Users/MichaelNielsen

20110729 20:26:31 nbsp Sitting inside the Covell Blvd. location looking at a sign in the window in front of which five young Asian women were sitting. A head was blocking tofu, the last word of Nobody ever does a highfive after eating tofu.  Then a couple of them did a highfive, seemingly unconscious of the sign. Users/BruceHansen
 When I have tofu tonight, I plan to do a high five after. Users/NickSchmalenberger

20111110 23:16:06 nbsp chili cheese fries & $1 chicken sandwich!!  Users/ndah100

20120111 22:43:04 nbsp As of January 2012, it seems like the Silo Carls Jr. will open at 10AM and no longer be serving breakfast. Farewell, hash brown nuggets! I hardly knew thee. Users/AbbYu

20121028 15:18:02 nbsp when did the east davis one open? Users/davisrox
Days ago, Hella days...
  Been there since I was a kid so at least early 90s Daubert
   Been there since I was a kid too, so Id push that back to early 80s. Users/jefftolentino
     I hear the 80s, 80s is on the table, do I hear 70s? 70s folks, can anybody recall the 70s?  80s going once...
      I say 70s.  My son considered getting a job there while in high school and he graduated in 1980. That would narrow it to between 76 and early 80.

20140209 22:06:22 nbsp Service is fast, obviously, but the sad thing is no drive thru :( The milkshakes are made with real ice cream and are delicious. The spicy chicken sandwich is a bargain. Users/YoungResident

20140213 18:25:13 nbsp It was a good time biking here one New Years Eve in 2008 with Users/SteveDavison, after I drank some champagne and other stuff. I think we used a 2 for 1 burger coupon. Users/NickSchmalenberger


text:
subroutine units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2)
  use amr_commons
  use hydro_commons
  use cooling_module
  real(dp)::scale_nH,scale_T2,scale_t,scale_v,scale_d,scale_l
  !-----------------------------------------------------------------------
  ! Conversion factors from user units into cgs units
  ! For gravity runs, make sure that G=1 in user units.
  !-----------------------------------------------------------------------

  ! scale_d converts mass density from user units into g/cc
  scale_d = units_density
  if(cosmo) scale_d = omega_m * rhoc *(h0/100.)**2 / aexp**3

  ! scale_t converts time from user units into seconds
  scale_t = units_time
  if(cosmo) scale_t = aexp**2 / (h0*1d5/3.08d24)

  ! scale_l converts distance from user units into cm
  scale_l = units_length
  if(cosmo) scale_l = aexp * boxlen_ini * 3.08d24 / (h0/100)

  ! scale_v converts velocity in user units into cm/s
  scale_v = scale_l / scale_t

  ! scale_T2 converts (P/rho) in user unit into (T/mu) in Kelvin
  scale_T2 = mH/kB * scale_v**2

  ! scale_nH converts rho in user units into nH in H/cc
  scale_nH = X/mH * scale_d

end subroutine units

max_forks_repo_name: rknop/amuse


text:
!| The `numeric_real` module provides constants and procedures related to floating-point number.
!
! The constants include minimum and maximum values of the 4- and 8-byte floating-point number types
! provided by the Fortran standard,
! parameters for specifying the result kind of the intrinsic procedures like `transfera, and the machine epsilons.
! 16-byte floating-point number type is not supported yet.
!
! The procedures include functions to convert a floating-point number to a string in a specified format.
!
module numeric_real
    use, intrinsic :: iso_fortran_env
    implicit none
    private
    public :: to_string

    !&<
    real(real32), public, parameter :: result_type_real32 = 0.0_real32
        !! A constant to represent 4-byte floating-point number.
        !! This constant is used to specify the result kind of the intrinsic functions.
    real(real64), public, parameter :: result_type_real64 = 0.0_real64
        !! A constant to represent 8-byte floating-point number.
        !! This constant is used to specify the result kind of the intrinsic functions.
    !&>

    real(real32), public, parameter :: Real32_Positive_Min = real(Z'00800000', real32)
        !! The positive minimum value of 4-byte floating-point number `≈+1.17549435E-38`
    real(real32), public, parameter :: Real32_Positive_Max = real(Z'7F7FFFFF', real32)
        !! The positive maximum value of 4-byte floating-point number `≈+3.40282347E+38`
    real(real32), public, parameter :: Real32_Negative_Min = real(Z'80800000', real32)
        !! The negative minimum value of 4-byte floating-point number `≈-1.17549435E-38`
    real(real32), public, parameter :: Real32_Negative_Max = real(Z'FF7FFFFF', real32)
        !! The negative maximum value of 4-byte floating-point number `≈-3.40282347E+38`

    real(real64), public, parameter :: Real64_Positive_Min = real(Z'0010000000000000', real64)
        !! The positive maximum value of 8-byte floating-point number `≈+1.7976931348623157E+308`
    real(real64), public, parameter :: Real64_Positive_Max = real(Z'7FEFFFFFFFFFFFFF', real64)
        !! The positive minimum value of 8-byte floating-point number `≈+2.2250738585072014E-308`
    real(real64), public, parameter :: Real64_Negative_Min = real(Z'8010000000000000', real64)
        !! The negative maximum value of 8-byte floating-point number `≈-1.7976931348623157E+308`
    real(real64), public, parameter :: Real64_Negative_Max = real(Z'FFEFFFFFFFFFFFFF', real64)
        !! The negative minimum value of 8-byte floating-point number `≈-2.2250738585072014E-308`

    real(real32), public, parameter :: Real32_Machine_Epsilon = epsilon(result_type_real32)
        !! The machine epsilon of 4-byte floating-point number `≈1.19209290E-07`
    real(real64), public, parameter :: Real64_Machine_Epsilon = epsilon(result_type_real64)
        !! The machine epsilon of 8-byte floating-point number `≈2.2204460492503131E-016`

    interface to_string
        procedure :: real32_to_string
        procedure :: real64_to_string
    end interface

    ! set the default format to recover the internal binary representation when reading decimal representation in ASCII format.
    character(*), private, parameter :: default_real32_format = '(e18.10E3)'
    character(*), private, parameter :: default_real64_format = '(e30.19E4)'

    integer(int32), private, parameter :: MaxLength_of_IOMessage = 255
    integer(int32), private, parameter :: MaxLength_of_Buffer = 255

contains
    !| Converts a 4-byte floating-point number to a string based on a format passed as an argument.
    ! When the format is not given, `'(e18.10E3)'` is used as the format
    ! to recover the internal binary representation when reading decimal representation in ASCII format.
    ! A zero-length string is returned when the convert is failed.
    function real32_to_string(r32, format) result(str)
        implicit none
        !&<
        real(real32), intent(in)            :: r32
            !! 4-byte floating-point number
        character(*), intent(in), optional  :: format
            !! format for string conversion
        !&>
        character(:), allocatable :: str

        character(MaxLength_of_Buffer) :: buffer
        integer :: IOStatus
        character(MaxLength_of_IOMessage) :: IOMessage

        if (present(format)) then
            write (buffer, format, iostat=IOStatus, iomsg=IOMessage) r32
        else
            write (buffer, default_real32_format, iostat=IOStatus, iomsg=IOMessage) r32
        end if
        if (IOStatus /= 0) then
            write (error_unit, '(A)') trim(IOMessage)
            buffer = "" ! returns a zero-length string
        end if

        str = trim(adjustl(buffer))
    end function real32_to_string

    !| Converts an 8-byte floating-point number to a string based on a format passed as an argument.
    ! When the format is not given, `'(e30.19E4)'` is used as the format
    ! to recover the internal binary representation when reading decimal representation in ASCII format.
    ! A zero-length string is returned when the convert is failed.
    function real64_to_string(r64, format) result(str)
        implicit none
        !&<
        real(real64), intent(in)            :: r64
            !! 8-byte floating-point number
        character(*), intent(in), optional  :: format
            !! format for string conversion
        !&>
        character(:), allocatable :: str

        character(MaxLength_of_Buffer) :: buffer
        integer :: IOStatus
        character(MaxLength_of_IOMessage) :: IOMessage

        if (present(format)) then
            write (buffer, format, iostat=IOStatus, iomsg=IOMessage) r64
        else
            write (buffer, default_real64_format, iostat=IOStatus, iomsg=IOMessage) r64
        end if
        if (IOStatus /= 0) then
            write (error_unit, '(A)') trim(IOMessage)
            buffer = "" ! returns a zero-length string
        end if

        str = trim(adjustl(buffer))
    end function real64_to_string
end module numeric_real

max_stars_repo_name: degawa/fortran-numeric

max_issues_repo_name: degawa/fortran-numeric

max_forks_repo_name: degawa/fortran-numeric


text:
Davis Math Tutoring strives to discover your childs math needs and set goals for his/her achievement. Parents and students have enough stress, which is why I offer tutoring services in the convenience of your home or mine. I focus on your childs math weaknesses and help develop the skills necessary for sucess.

I have experience tutoring all ages and abilities from elementary  high school math and beyond!

Tutoring Offered:
Addition
Subtraction
Division
Multiplication
Fractions
Arithmetic
Algebra
Geometry
Trigonometry
Calculus
AP Courses
SAT Math Prep

Julie Schiffman received her Bachelors degree in Applied Mathematics with honors from UC San Diego in 2000. She worked for the UCSD Math Department as an undergraduate Teaching Assistant and was also awarded an Errett Bishop Math Scholarship. She earned an MBA in 2007 while living abroad.


Module Md_Types_Numeriques

  !Integer, parameter :: sp = selected_real_kind(6,37)   ! nombre de digit et puissance max
  !Integer, parameter :: dp = selected_real_kind(15,307)

  Integer, parameter :: sp  = kind(1.0)
  Integer, parameter :: dp  = kind(1.0d0)
  Integer, parameter :: spc = kind((1.0,1.0))
  Integer, parameter :: dpc = kind((1.0d0,1.0d0))

End module Md_Types_Numeriques