text:
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! H0 X
! H0 X   libAtoms+QUIP: atomistic simulation library
! H0 X
! H0 X   Portions of this code were written by
! H0 X     Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! H0 X     Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! H0 X
! H0 X   Copyright 2006-2010.
! H0 X
! H0 X   These portions of the source code are released under the GNU General
! H0 X   Public License, version 2, http://www.gnu.org/copyleft/gpl.html
! H0 X
! H0 X   If you would like to license the source code under different terms,
! H0 X   please contact Gabor Csanyi, gabor@csanyi.net
! H0 X
! H0 X   Portions of this code were written by Noam Bernstein as part of
! H0 X   his employment for the U.S. Government, and are not subject
! H0 X   to copyright in the USA.
! H0 X
! H0 X
! H0 X   When using this software, please cite the following reference:
! H0 X
! H0 X   http://www.libatoms.org
! H0 X
! H0 X  Additional contributions by
! H0 X    Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! H0 X
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

!X
!X Histogram helper functions in 1D
!X
!X Computes periodic and non-periodic histograms within certain bounds.
!X Interpolation can be chosen from
!X   INTERP_LINEAR:  Linear interpolation between neigboring grid
!X                   points
!X   INTERP_LUCY:    Lucy function interpolation between neighboring
!X                   grid points
!X TODO:
!X   INTERP_NONE:    No interpolation, simply sort into bins
!X
!X Provides routine for MPI communication.
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

#include "error.inc"

#define ASSERT(x)

module histogram1d_module
  use Error_module
  use System_module
  use Units_module
  use MPI_context_module

  private

  public :: INTERP_LINEAR, INTERP_LUCY
  integer, parameter  :: INTERP_LINEAR  = 0
  integer, parameter  :: INTERP_LUCY    = 1

  public :: Histogram1D
  type Histogram1D

     !
     ! General stuff
     !

     integer            :: interp     ! INTERP_LINEAR, INTERP_LUCY

     integer            :: n = -1     ! number of bins
     real(DP)           :: min_b      ! minimum value
     real(DP)           :: max_b      ! maximun value
     real(DP)           :: db         ! difference
     real(DP)           :: dbin       ! bin size

     logical            :: periodic

     !
     ! For smoothing function
     !

     real(DP)           :: sigma

     real(DP)           :: fac1
!     real(DP)           :: sigma_sq

     !
     ! Values
     !

     real(DP), allocatable  :: x(:)       ! x-values

     real(DP), allocatable  :: h(:)       ! values
     real(DP), allocatable  :: h_sq(:)    ! second moment

     real(DP), allocatable  :: h1(:)      ! histogram with norm 1 (number of values which have been added to each bin)

     integer                :: ng
     real(DP), pointer      :: smoothing_func(:) => NULL()

  endtype Histogram1D

  !
  ! Interface definition
  !

  public :: initialise
  interface initialise
     module procedure histogram1d_initialise, histogram1d_initialise_from_histogram
  endinterface

  public :: finalise
  interface finalise
     module procedure histogram1d_finalise
  endinterface

  public :: clear
  interface clear
     module procedure histogram1d_clear
  endinterface

  public :: set_bounds
  interface set_bounds
     module procedure histogram1d_set_bounds
  endinterface

  public :: write
  interface write
     module procedure histogram1d_write, histogram1d_write_mult
     module procedure histogram1d_write_character_fn, histogram1d_write_mult_character_fn
  endinterface

  public :: add
  interface add
     module procedure histogram1d_add, histogram1d_add_vals, histogram1d_add_vals_norms
     module procedure histogram1d_add_vals_mask, histogram1d_add_vals_norm_mask, histogram1d_add_vals_norms_mask
     module procedure histogram1d_add_histogram, histogram1d_add_mult_histograms
  endinterface

  public :: add_range
  interface add_range
     module procedure histogram1d_add_range, histogram1d_add_range_vals, histogram1d_add_range_vals_norms
  endinterface

  public :: average
  interface average
     module procedure histogram1d_average
  endinterface

  public :: entropy
  interface entropy
     module procedure histogram1d_entropy
  endinterface

  public :: expectation_value
  interface expectation_value
     module procedure histogram1d_expectation_value
  endinterface

  public :: mul
  interface mul
     module procedure histogram1d_mul, histograms_mul, histograms2_mul, histogram1d_mul_vals
  endinterface

  public :: div
  interface div
     module procedure histogram1d_div, histograms_div, histograms2_div, histogram1d_div_vals
  endinterface

  public :: normalize
  interface normalize
     module procedure histogram1d_normalize
  endinterface

  public :: reduce
  interface reduce
     module procedure histogram1d_reduce
  endinterface

  public :: smooth
  interface smooth
     module procedure histogram1d_smooth
  endinterface

  public :: sum_in_place
  interface sum_in_place
     module procedure histogram1d_sum_in_place
  endinterface

! Memory estimation not yet implemented in libAtoms
#if 0
  interface log_memory_estimate
     module procedure log_memory_estimate_histogram, log_memory_estimate_histogram2, log_memory_estimate_histogram3
  endinterface
#endif

contains

  !% Initialize the histogram
  !% Default is linear interpolation, non-periodic
  elemental subroutine histogram1d_initialise(this, n, min_b, max_b, sigma, &
       periodic)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    integer, intent(in)               :: n
    real(DP), intent(in)              :: min_b
    real(DP), intent(in)              :: max_b
    real(DP), intent(in), optional    :: sigma
    logical, intent(in), optional     :: periodic

    ! ---

    call finalise(this)

    this%n  = n

    if (present(sigma) .and. sigma > 0.0_DP) then

       this%interp    = INTERP_LUCY

       this%sigma     = sigma

!       this%fac1      = 1.0_DP/(sqrt(2*PI)*this%sigma)
!       this%sigma_sq  = 2*sigma**2
       this%fac1      = 5/(4*this%sigma)

    else

       this%interp    = INTERP_LINEAR

    endif

    if (present(periodic)) then
       this%periodic  = periodic
    else
       this%periodic  = .false.
    endif

    this%smoothing_func => NULL()

    allocate(this%x(n))
    allocate(this%h(n))
    allocate(this%h_sq(n))
    allocate(this%h1(n))

    this%min_b = min_b - 1.0_DP
    this%max_b = max_b - 1.0_DP

    call histogram1d_set_bounds(this, min_b, max_b)

    call histogram1d_clear(this)

  endsubroutine histogram1d_initialise

  !% Initialize the histogram
  subroutine histogram1d_initialise_from_histogram(this, that)
    implicit none

    type(Histogram1D), intent(out)  :: this
    type(Histogram1D), intent(in )  :: that

    ! ---

    call finalise(this)

    this%n         = that%n

    this%interp    = that%interp
    this%sigma     = that%sigma

    if (this%interp == INTERP_LUCY) then
       this%fac1      = 5/(4*this%sigma)
    endif

    this%periodic  = that%periodic

    this%smoothing_func => NULL()

    allocate(this%x(this%n))
    allocate(this%h(this%n))
    allocate(this%h_sq(this%n))
    allocate(this%h1(this%n))

    this%min_b  = that%min_b - 1.0_DP
    this%max_b  = that%max_b - 1.0_DP

    call histogram1d_set_bounds(this, that%min_b, that%max_b)

    call histogram1d_clear(this)

  endsubroutine histogram1d_initialise_from_histogram

  !% Clear histogram
  elemental subroutine histogram1d_clear(this)
    implicit none

    type(Histogram1D), intent(inout)  :: this

    ! ---

    this%h1    =  0.0_DP
    this%h     = 0.0_DP
    this%h_sq  = 0.0_DP

  endsubroutine histogram1d_clear

  !% Set histogram bounds
  elemental subroutine histogram1d_set_bounds(this, min_b, max_b)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: min_b
    real(DP), intent(in)              :: max_b

    ! ---

    integer   :: i
    real(DP)  :: x

    ! ---

    if (this%min_b /= min_b .or. this%max_b /= max_b) then

       this%min_b  = min_b
       this%max_b  = max_b
       this%db     = max_b - min_b

       if (this%periodic) then
          this%dbin   = this%db/this%n
       else
          this%dbin   = this%db/(this%n-1)
       endif

       do i = 1, this%n
          this%x(i)  = this%min_b + this%dbin*(i-0.5_DP)
       enddo

       !
       ! Tabulate the smoothing function
       !

       if (this%interp == INTERP_LUCY) then

          this%ng  = floor(this%sigma/this%dbin)

          if (.not. associated(this%smoothing_func)) then
             allocate(this%smoothing_func(-this%ng:this%ng))
          else if (ubound(this%smoothing_func, 1) /= this%ng) then
             deallocate(this%smoothing_func)
             allocate(this%smoothing_func(-this%ng:this%ng))
          endif

          do i = -this%ng, this%ng
!          this%smoothing_func(i)  = this%fac1 * exp(-((this%dbin*i)**2)/this%sigma_sq)
             x = abs( ( this%dbin*i )/this%sigma )
             this%smoothing_func(i)  = this%fac1 * (1+3*x)*(1-x)**3
          enddo

       endif

    endif

  endsubroutine histogram1d_set_bounds

  !% Delete a histogram table
  elemental subroutine histogram1d_finalise(this)
    implicit none

    type(Histogram1D), intent(inout)  :: this

    ! ---

#if 0
    if (allocated(this%x)) then
       deallocate(this%x)
    endif
    if (allocated(this%h)) then
       deallocate(this%h)
    endif
    if (allocated(this%h_sq)) then
       deallocate(this%h_sq)
    endif
    if (allocated(this%h1)) then
       deallocate(this%h1)
    endif

    if (associated(this%smoothing_func)) then
       deallocate(this%smoothing_func)
    endif
#endif

  endsubroutine histogram1d_finalise

  !% Output the histogram table to a file (by file unit)
  subroutine histogram1d_write_mult(this, file, xvalues, header, error)
    implicit none

    type(Histogram1D), intent(in)       :: this(:)
    type(InOutput), intent(in)          :: file
    logical, intent(in), optional       :: xvalues
    character(*), intent(in), optional  :: header(lbound(this, 1):ubound(this, 1))
    integer, intent(out), optional      :: error

    ! ---

    integer        :: i, j
    character(80)  :: fmt
    character(20)  :: ext_header(lbound(this, 1):ubound(this, 1))

    integer        :: n
    real(DP)       :: min_b, max_b, dbin, y(lbound(this, 1):ubound(this, 1))

    ! ---

    INIT_ERROR(error)

    n      = this(lbound(this, 1))%n
    min_b  = this(lbound(this, 1))%min_b
    max_b  = this(lbound(this, 1))%max_b
    dbin   = this(lbound(this, 1))%dbin

    do i = lbound(this, 1)+1, ubound(this, 1)
       if (this(i)%n /= n) then
          RAISE_ERROR("Number of histogram bins do not match.", error)
       endif

       if (this(i)%min_b /= min_b) then
          RAISE_ERROR("*min_b*s do not match.", error)
       endif

       if (this(i)%max_b /= max_b) then
          RAISE_ERROR("*max_b*s do not match.", error)
       endif
    enddo

    if (present(header)) then
       do i = lbound(this, 1), ubound(this, 1)
          write (fmt, '(I2.2)')  i+2
          ext_header(i)  = trim(fmt) // ":" // trim(header(i))
          ext_header(i)  = adjustr(ext_header(i))
       enddo

       write (fmt, '(A,I4.4,A)')  "(A5,5X,A20,", ubound(this, 1)-lbound(this, 1)+1, "A20)"
       write (file%unit, fmt)  "#01:i", "02:x", ext_header
    endif

    if (present(xvalues) .and. .not. xvalues) then
       write (fmt, '(A,I4.4,A)')  "(A,", n, "ES20.10)"
       write (file%unit, trim(fmt))  "# ", this(lbound(this, 1))%x
       write (fmt, '(A,I4.4,A)')  "(I10,", ubound(this, 1)-lbound(this, 1)+1, "ES20.10)"
    else
       write (fmt, '(A,I4.4,A)')  "(I10,", ubound(this, 1)-lbound(this, 1)+2, "ES20.10)"
    endif

    do i = 1, n
       do j = lbound(this, 1), ubound(this, 1)
          y(j) = this(j)%h(i)
       enddo

       if (present(xvalues) .and. .not. xvalues) then
          write (file%unit, trim(fmt))  i, y
       else
          write (file%unit, trim(fmt))  i, this(lbound(this, 1))%x(i), y
       endif
    enddo

  endsubroutine histogram1d_write_mult

  !% Output the histogram table to a file (by file name)
  subroutine histogram1d_write_mult_character_fn(this, fn, xvalues, header, error)
    implicit none

    type(Histogram1D), intent(in)       :: this(:)
    character(*), intent(in)            :: fn
    logical, intent(in), optional       :: xvalues
    character(*), intent(in), optional  :: header(lbound(this, 1):ubound(this, 1))
    integer, intent(out), optional      :: error

    ! ---

    type(InOutput)  :: file

    ! ---

    INIT_ERROR(error)

    call initialise(file, fn, action=OUTPUT, error=error)
    PASS_ERROR(error)
    call histogram1d_write_mult(this, file, xvalues, header, error=error)
    call finalise(file)
    PASS_ERROR(error)

  endsubroutine histogram1d_write_mult_character_fn

  !% Output the histogram table to a file (by file unit)
  subroutine histogram1d_write(this, file, xvalues, header, error)
    implicit none

    type(Histogram1D), intent(in)       :: this
    type(InOutput), intent(in)          :: file
    logical, intent(in), optional       :: xvalues
    character(*), intent(in), optional  :: header
    integer, intent(out), optional      :: error

    ! ---

    INIT_ERROR(error)

    if (present(header)) then
       call histogram1d_write_mult((/ this /), file, xvalues, (/ header /), error=error)
       PASS_ERROR(error)
    else
       call histogram1d_write_mult((/ this /), file, xvalues, error=error)
       PASS_ERROR(error)
    endif

  endsubroutine histogram1d_write

  !% Output the histogram table to a file (by file name)
  subroutine histogram1d_write_character_fn(this, fn, xvalues, header, error)
    implicit none

    type(Histogram1D), intent(in)       :: this
    character(*), intent(in)            :: fn
    logical, intent(in), optional       :: xvalues
    character(*), intent(in), optional  :: header
    integer, intent(out), optional      :: error

    ! ---

    INIT_ERROR(error)

    if (present(header)) then
       call histogram1d_write_mult_character_fn((/ this /), fn, xvalues, (/ header /), error=error)
       PASS_ERROR(error)
    else
       call histogram1d_write_mult_character_fn((/ this /), fn, xvalues, error=error)
       PASS_ERROR(error)
    endif

  endsubroutine histogram1d_write_character_fn

  !% Add a value to the histogram with linear interpolation
  subroutine histogram1d_add_linear(this, val, norm)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: val
    real(DP), intent(in)              :: norm

    ! ---

    integer   :: i1, i2
    real(DP)  :: d1, d2

    ! ---

    i1  = int(floor((val-this%min_b)/this%dbin))+1
    i2  = i1+1

    d1  = ( (i2-1) - (val-this%min_b)/this%dbin )/this%dbin
    d2  = ( (val-this%min_b)/this%dbin - (i1-1) )/this%dbin

    if ((i1 >= 0 .and. i1 <= this%n) .or. this%periodic) then
       ASSERT(d1 >= 0.0_DP .and. d1*this%dbin <= 1.0_DP)

       if (this%periodic) then
          i1  = modulo(i1-1, this%n)+1
       else
          if (i1 == 0) then
             i1 = 1
          endif
       endif

       this%h1(i1)    = this%h1(i1)   + d1
       this%h(i1)     = this%h(i1)    + d1*norm
       this%h_sq(i1)  = this%h_sq(i1) + d1*norm**2
    endif

    if ((i2 >= 1 .and. i2 <= this%n+1) .or. this%periodic) then
       ASSERT(d2 >= 0.0_DP .and. d2*this%dbin <= 1.0_DP)

       if (this%periodic) then
          i2  = modulo(i2-1, this%n)+1
       else
          if (i2 == this%n+1) then
             i2 = this%n
          endif
       endif

       this%h1(i2)    = this%h1(i2)   + d2
       this%h(i2)     = this%h(i2)    + d2*norm
       this%h_sq(i2)  = this%h_sq(i2) + d2*norm**2
    endif

  endsubroutine histogram1d_add_linear

  !% Add multiple values to the histogram with linear interpolation
  subroutine histogram1d_add_linear_vals(this, vals, norm)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: vals(:)
    real(DP), intent(in)              :: norm

    ! ---

    integer   :: i1(lbound(vals, 1):ubound(vals, 1)), i2(lbound(vals, 1):ubound(vals, 1))
    real(DP)  :: d1(lbound(vals, 1):ubound(vals, 1)), d2(lbound(vals, 1):ubound(vals, 1))

    integer   :: i

    ! ---

    i1  = int(floor((vals-this%min_b)/this%dbin))+1
    i2  = i1+1

    d1  = ( (i2-1) - (vals-this%min_b)/this%dbin )/this%dbin
    d2  = ( (vals-this%min_b)/this%dbin - (i1-1) )/this%dbin

    ASSERT(all(d1 >= 0.0_DP .and. d1*this%dbin <= 1.0_DP))
    ASSERT(all(d2 >= 0.0_DP .and. d2*this%dbin <= 1.0_DP))

    if (this%periodic) then

       i1  = modulo(i1-1, this%n)+1
       i2  = modulo(i2-1, this%n)+1

       do i = lbound(vals, 1), ubound(vals, 1)
          this%h1(i1(i))    = this%h1(i1(i))   + d1(i)
          this%h1(i2(i))    = this%h1(i2(i))   + d2(i)

          this%h(i1(i))     = this%h(i1(i))    + d1(i)*norm
          this%h(i2(i))     = this%h(i2(i))    + d2(i)*norm

          this%h_sq(i1(i))  = this%h_sq(i1(i)) + d1(i)*norm**2
          this%h_sq(i2(i))  = this%h_sq(i2(i)) + d2(i)*norm**2
       enddo

    else

       do i = lbound(vals, 1), ubound(vals, 1)
          if (i1(i) >= 1 .and. i1(i) <= this%n) then
             this%h1(i1(i))    = this%h1(i1(i))   + d1(i)
             this%h(i1(i))     = this%h(i1(i))    + d1(i)*norm
             this%h_sq(i1(i))  = this%h_sq(i1(i)) + d1(i)*norm**2
          else if (i1(i) == 0) then
             this%h1(1)    = this%h1(1)   + d1(i)
             this%h(1)     = this%h(1)    + d1(i)*norm
             this%h_sq(1)  = this%h_sq(1) + d1(i)*norm**2
          endif

          if (i2(i) >= 1 .and. i2(i) <= this%n) then
             this%h1(i2(i))    = this%h1(i2(i))   + d2(i)
             this%h(i2(i))     = this%h(i2(i))    + d2(i)*norm
             this%h_sq(i2(i))  = this%h_sq(i2(i)) + d2(i)*norm**2
          else if (i2(i) == this%n+1) then
             this%h1(this%n)    = this%h1(this%n)   + d2(i)
             this%h(this%n)     = this%h(this%n)    + d2(i)*norm
             this%h_sq(this%n)  = this%h_sq(this%n) + d2(i)*norm**2
          endif
       enddo

    endif

  endsubroutine histogram1d_add_linear_vals

  !% Add multiple values to the histogram with linear interpolation
  subroutine histogram1d_add_linear_vals_norms(this, vals, norms)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: vals(:)
    real(DP), intent(in)              :: norms(:)

    ! ---

    integer   :: i1(lbound(vals, 1):ubound(vals, 1)), i2(lbound(vals, 1):ubound(vals, 1))
    real(DP)  :: d1(lbound(vals, 1):ubound(vals, 1)), d2(lbound(vals, 1):ubound(vals, 1))

    real(DP)  :: d1n(lbound(vals, 1):ubound(vals, 1)), d2n(lbound(vals, 1):ubound(vals, 1))
    real(DP)  :: d1nn(lbound(vals, 1):ubound(vals, 1)), d2nn(lbound(vals, 1):ubound(vals, 1))

    integer   :: i

    ! ---

    i1  = int(floor((vals-this%min_b)/this%dbin))+1
    i2  = i1+1

    d1  = ( (i2-1) - (vals-this%min_b)/this%dbin )/this%dbin
    d2  = ( (vals-this%min_b)/this%dbin - (i1-1) )/this%dbin

    ASSERT(all(d1 >= 0.0_DP .and. d1*this%dbin <= 1.0_DP))
    ASSERT(all(d2 >= 0.0_DP .and. d2*this%dbin <= 1.0_DP))

    d1n   = d1*norms
    d2n   = d2*norms

    d1nn  = d1n*norms
    d2nn  = d2n*norms

    if (this%periodic) then

       i1  = modulo(i1-1, this%n)+1
       i2  = modulo(i2-1, this%n)+1

       do i = lbound(vals, 1), ubound(vals, 1)
          this%h1(i1(i))    = this%h1(i1(i))  + d1(i)
          this%h1(i2(i))    = this%h1(i2(i))  + d2(i)

          this%h(i1(i))     = this%h(i1(i))   + d1n(i)
          this%h(i2(i))     = this%h(i2(i))   + d2n(i)

          this%h_sq(i1(i))  = this%h_sq(i1(i)) + d1nn(i)
          this%h_sq(i2(i))  = this%h_sq(i2(i)) + d2nn(i)
       enddo

    else

       do i = lbound(vals, 1), ubound(vals, 1)
          if (i1(i) >= 1 .and. i1(i) <= this%n) then
             this%h1(i1(i))    = this%h1(i1(i))   + d1(i)
             this%h(i1(i))     = this%h(i1(i))    + d1n(i)
             this%h_sq(i1(i))  = this%h_sq(i1(i)) + d1nn(i)
          else if (i1(i) == 0) then
             this%h1(1)    = this%h1(1)   + d1(i)
             this%h(1)     = this%h(1)    + d1n(i)
             this%h_sq(1)  = this%h_sq(1) + d1nn(i)
          endif

          if (i2(i) >= 1 .and. i2(i) <= this%n) then
             this%h1(i2(i))    = this%h1(i2(i))   + d2(i)
             this%h(i2(i))     = this%h(i2(i))    + d2n(i)
             this%h_sq(i2(i))  = this%h_sq(i2(i)) + d2nn(i)
          else if (i2(i) == this%n+1) then
             this%h1(this%n)    = this%h1(this%n)   + d2(i)
             this%h(this%n)     = this%h(this%n)    + d2n(i)
             this%h_sq(this%n)  = this%h_sq(this%n) + d2nn(i)
          endif
       enddo

    endif

  endsubroutine histogram1d_add_linear_vals_norms

  !% Add multiple values to the histogram with linear interpolation
  !% and an additional mask
  subroutine histogram1d_add_linear_vals_norm_mask(this, vals, norm, mask)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: vals(:)
    real(DP), intent(in)              :: norm
    logical, intent(in)               :: mask(:)

    ! ---

    integer   :: i1(lbound(vals, 1):ubound(vals, 1)), i2(lbound(vals, 1):ubound(vals, 1))
    real(DP)  :: d1(lbound(vals, 1):ubound(vals, 1)), d2(lbound(vals, 1):ubound(vals, 1))

    real(DP)  :: d1n(lbound(vals, 1):ubound(vals, 1)), d2n(lbound(vals, 1):ubound(vals, 1))
    real(DP)  :: d1nn(lbound(vals, 1):ubound(vals, 1)), d2nn(lbound(vals, 1):ubound(vals, 1))

    integer   :: i

    ! ---

    i1  = int(floor((vals-this%min_b)/this%dbin))+1
    i2  = i1+1

    d1  = ( (i2-1) - (vals-this%min_b)/this%dbin )/this%dbin
    d2  = ( (vals-this%min_b)/this%dbin - (i1-1) )/this%dbin

    ASSERT(all(d1 >= 0.0_DP .and. d1*this%dbin <= 1.0_DP))
    ASSERT(all(d2 >= 0.0_DP .and. d2*this%dbin <= 1.0_DP))

    d1n   = d1*norm
    d2n   = d2*norm

    d1nn  = d1n*norm
    d2nn  = d2n*norm

    if (this%periodic) then

       i1  = modulo(i1-1, this%n)+1
       i2  = modulo(i2-1, this%n)+1

       do i = lbound(vals, 1), ubound(vals, 1)
          if (mask(i)) then
             this%h1(i1(i))    = this%h1(i1(i))  + d1(i)
             this%h1(i2(i))    = this%h1(i2(i))  + d2(i)

             this%h(i1(i))     = this%h(i1(i))   + d1n(i)
             this%h(i2(i))     = this%h(i2(i))   + d2n(i)

             this%h_sq(i1(i))  = this%h_sq(i1(i)) + d1nn(i)
             this%h_sq(i2(i))  = this%h_sq(i2(i)) + d2nn(i)
          endif
       enddo

    else

       do i = lbound(vals, 1), ubound(vals, 1)
          if (mask(i)) then
             if (i1(i) >= 1 .and. i1(i) <= this%n) then
                this%h1(i1(i))    = this%h1(i1(i))   + d1(i)
                this%h(i1(i))     = this%h(i1(i))    + d1n(i)
                this%h_sq(i1(i))  = this%h_sq(i1(i)) + d1nn(i)
             else if (i1(i) == 0) then
                this%h1(1)    = this%h1(1)   + d1(i)
                this%h(1)     = this%h(1)    + d1n(i)
                this%h_sq(1)  = this%h_sq(1) + d1nn(i)
             endif

             if (i2(i) >= 1 .and. i2(i) <= this%n) then
                this%h1(i2(i))    = this%h1(i2(i))   + d2(i)
                this%h(i2(i))     = this%h(i2(i))    + d2n(i)
                this%h_sq(i2(i))  = this%h_sq(i2(i)) + d2nn(i)
             else if (i2(i) == this%n+1) then
                this%h1(this%n)    = this%h1(this%n)   + d2(i)
                this%h(this%n)     = this%h(this%n)    + d2n(i)
                this%h_sq(this%n)  = this%h_sq(this%n) + d2nn(i)
             endif
          endif
       enddo

    endif

  endsubroutine histogram1d_add_linear_vals_norm_mask

  !% Add multiple values to the histogram with linear interpolation
  !% and an additional mask
  subroutine histogram1d_add_linear_vals_norms_mask(this, vals, norms, mask)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: vals(:)
    real(DP), intent(in)              :: norms(:)
    logical, intent(in)               :: mask(:)

    ! ---

    integer   :: i1(lbound(vals, 1):ubound(vals, 1)), i2(lbound(vals, 1):ubound(vals, 1))
    real(DP)  :: d1(lbound(vals, 1):ubound(vals, 1)), d2(lbound(vals, 1):ubound(vals, 1))

    real(DP)  :: d1n(lbound(vals, 1):ubound(vals, 1)), d2n(lbound(vals, 1):ubound(vals, 1))
    real(DP)  :: d1nn(lbound(vals, 1):ubound(vals, 1)), d2nn(lbound(vals, 1):ubound(vals, 1))

    integer   :: i

    ! ---

    i1  = int(floor((vals-this%min_b)/this%dbin))+1
    i2  = i1+1

    d1  = ( (i2-1) - (vals-this%min_b)/this%dbin )/this%dbin
    d2  = ( (vals-this%min_b)/this%dbin - (i1-1) )/this%dbin

    ASSERT(all(d1 >= 0.0_DP .and. d1*this%dbin <= 1.0_DP))
    ASSERT(all(d2 >= 0.0_DP .and. d2*this%dbin <= 1.0_DP))

    d1n   = d1*norms
    d2n   = d2*norms

    d1nn  = d1n*norms
    d2nn  = d2n*norms

    if (this%periodic) then

       i1  = modulo(i1-1, this%n)+1
       i2  = modulo(i2-1, this%n)+1

       do i = lbound(vals, 1), ubound(vals, 1)
          if (mask(i)) then
             this%h1(i1(i))    = this%h1(i1(i))  + d1(i)
             this%h1(i2(i))    = this%h1(i2(i))  + d2(i)

             this%h(i1(i))     = this%h(i1(i))   + d1n(i)
             this%h(i2(i))     = this%h(i2(i))   + d2n(i)

             this%h_sq(i1(i))  = this%h_sq(i1(i)) + d1nn(i)
             this%h_sq(i2(i))  = this%h_sq(i2(i)) + d2nn(i)
          endif
       enddo

    else

       do i = lbound(vals, 1), ubound(vals, 1)
          if (mask(i)) then
             if (i1(i) >= 1 .and. i1(i) <= this%n) then
                this%h1(i1(i))    = this%h1(i1(i))   + d1(i)
                this%h(i1(i))     = this%h(i1(i))    + d1n(i)
                this%h_sq(i1(i))  = this%h_sq(i1(i)) + d1nn(i)
             else if (i1(i) == 0) then
                this%h1(1)    = this%h1(1)   + d1(i)
                this%h(1)     = this%h(1)    + d1n(i)
                this%h_sq(1)  = this%h_sq(1) + d1nn(i)
             endif

             if (i2(i) >= 1 .and. i2(i) <= this%n) then
                this%h1(i2(i))    = this%h1(i2(i))   + d2(i)
                this%h(i2(i))     = this%h(i2(i))    + d2n(i)
                this%h_sq(i2(i))  = this%h_sq(i2(i)) + d2nn(i)
             else if (i2(i) == this%n+1) then
                this%h1(this%n)    = this%h1(this%n)   + d2(i)
                this%h(this%n)     = this%h(this%n)    + d2n(i)
                this%h_sq(this%n)  = this%h_sq(this%n) + d2nn(i)
             endif
          endif
       enddo

    endif

  endsubroutine histogram1d_add_linear_vals_norms_mask

  !% Add a value which is broadened by a smoothing function to the histogram
  subroutine histogram1d_add_smoothed(this, val, norm)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: val
    real(DP), intent(in)              :: norm

    ! ---

    real(DP)  :: expvals(this%n)

    integer   :: i1, i2, g1, g2
    real(DP)  :: x2

    ! ---

    x2  = (val-this%min_b)/this%dbin + 0.5_DP
    i1  = floor(x2)
    x2  = x2 - i1

    g1  = -this%ng
    g2  = this%ng-1

    expvals(1:2*this%ng)  = x2*this%smoothing_func(g1:g2) + (1.0_DP-x2)*this%smoothing_func(g1+1:g2+1)

    i1  = i1 - this%ng+1
    i2  = i1 + 2*this%ng - 1

    !
    ! Cut smoothing function if we're at the border
    !

    g1  = 1
    g2  = 2*this%ng

    if (i1 < 1) then
       g1 = g1 + (1 - i1)
       i1 = 1
    endif

    if (i2 > this%n) then
       g2 = g2 - (i2 - this%n)
       i2 = this%n
    endif

    this%h1(i1:i2)    = this%h1(i1:i2)   + expvals(g1:g2)
    this%h(i1:i2)     = this%h(i1:i2)    + expvals(g1:g2)*norm
    this%h_sq(i1:i2)  = this%h_sq(i1:i2) + expvals(g1:g2)*norm**2

  endsubroutine histogram1d_add_smoothed

  !**********************************************************************
  ! Add a value
  !**********************************************************************
  subroutine histogram1d_add(this, val, norm)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: val
    real(DP), intent(in), optional    :: norm

    ! ---

    real(DP)  :: n

    ! ---

    if (present(norm)) then
       n  = norm
    else
       n  = 1.0_DP
    endif

    if (this%interp == INTERP_LINEAR) then
       call histogram1d_add_linear(this, val, n)
    else
       call histogram1d_add_smoothed(this, val, n)
    endif

  endsubroutine histogram1d_add

  !% Add a list of values to the histogram
  subroutine histogram1d_add_vals(this, vals, norm)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: vals(:)
    real(DP), intent(in), optional    :: norm

    ! ---

    integer   :: i
    real(DP)  :: n

    ! ---

    if (present(norm)) then
       n  = norm
    else
       n  = 1.0_DP
    endif

    if (this%interp == INTERP_LINEAR) then

       call histogram1d_add_linear_vals(this, vals,  n)

    else

       do i = lbound(vals, 1), ubound(vals, 1)
          call histogram1d_add_smoothed(this, vals(i), n)
       enddo

    endif

  endsubroutine histogram1d_add_vals

  !% Add a list of values to the histogram
  subroutine histogram1d_add_vals_mask(this, vals, mask)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: vals(:)
    logical, intent(in)               :: mask(:)

    ! ---

    integer   :: i

    ! ---

    if (this%interp == INTERP_LINEAR) then

       do i = lbound(vals, 1), ubound(vals, 1)
          if (mask(i)) then
             call histogram1d_add_linear(this, vals(i), 1.0_DP)
          endif
       enddo

    else

       do i = lbound(vals, 1), ubound(vals, 1)
          if (mask(i)) then
             call histogram1d_add_smoothed(this, vals(i), 1.0_DP)
          endif
       enddo

    endif

  endsubroutine histogram1d_add_vals_mask

  !% Add a list of values to the histogram
  subroutine histogram1d_add_vals_norms(this, vals, norms)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: vals(:)
    real(DP), intent(in)              :: norms(:)

    ! ---

    integer   :: i

    ! ---

    if (this%interp == INTERP_LINEAR) then

       call histogram1d_add_linear_vals_norms(this, vals, norms)

    else

       do i = lbound(vals, 1), ubound(vals, 1)
          call histogram1d_add_smoothed(this, vals(i), norms(i))
       enddo

    endif

  endsubroutine histogram1d_add_vals_norms

  !% Add a list of values to the histogram
  subroutine histogram1d_add_vals_norm_mask(this, vals, norm, mask)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: vals(:)
    real(DP), intent(in)              :: norm
    logical, intent(in)               :: mask(:)

    ! ---

    integer   :: i

    ! ---

    if (this%interp == INTERP_LINEAR) then

       call histogram1d_add_linear_vals_norm_mask(this, vals, norm, mask)

    else

       do i = lbound(vals, 1), ubound(vals, 1)
          if (mask(i)) then
             call histogram1d_add_smoothed(this, vals(i), norm)
          endif
       enddo

    endif

  endsubroutine histogram1d_add_vals_norm_mask

  !% Add a list of values to the histogram
  subroutine histogram1d_add_vals_norms_mask(this, vals, norms, mask)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: vals(:)
    real(DP), intent(in)              :: norms(:)
    logical, intent(in)               :: mask(:)

    ! ---

    integer   :: i

    ! ---

    if (this%interp == INTERP_LINEAR) then

       call histogram1d_add_linear_vals_norms_mask(this, vals, norms, mask)

    else

       do i = lbound(vals, 1), ubound(vals, 1)
          if (mask(i)) then
             call histogram1d_add_smoothed(this, vals(i), norms(i))
          endif
       enddo

    endif

  endsubroutine histogram1d_add_vals_norms_mask

  !% Add a range of values to the histogram (linear interpolation only)
  subroutine histogram1d_add_range(this, vala, valb, norm, error)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: vala
    real(DP), intent(in)              :: valb
    real(DP), intent(in), optional    :: norm
    integer, intent(out), optional    :: error

    ! ---

    integer   :: i1a, i2a, i1b, i2b, i
    real(DP)  :: d1a, d1b, d2a, d2b
    real(DP)  :: va, vb, fac, n, n_sq

    ! ---

    INIT_ERROR(error)

    if (this%interp /= INTERP_LINEAR) then
       RAISE_ERROR("*histogram1d_add_range* can only be used with linear interpolation.", error)
    endif

    n  = 1.0_DP
    if (present(norm)) then
       n  = norm
    endif
    n_sq  = n**2

    if (vala > valb) then
       va  = valb
       vb  = vala
    else
       va  = vala
       vb  = valb
    endif

    i1a  = int(floor((va-this%min_b)/this%dbin))+1
    i2a  = i1a+1

    d1a  = ( (i2a-1) - (va-this%min_b)/this%dbin )
    d2a  = ( (va-this%min_b)/this%dbin - (i1a-1) )

    i1b  = int(floor((vb-this%min_b)/this%dbin))+1
    i2b  = i1b+1

    d1b  = ( (i2b-1) - (vb-this%min_b)/this%dbin )
    d2b  = ( (vb-this%min_b)/this%dbin - (i1b-1) )

    if (va == vb) then

       this%h1(i1a)    = this%h1(i1a)   + d1a
       this%h(i1a)     = this%h(i1a)    + d1a*n
       this%h_sq(i1a)  = this%h_sq(i1a) + d1a*n_sq

       this%h1(i2a)    = this%h1(i2a)   + d2a
       this%h(i2a)     = this%h(i2a)    + d2a*n
       this%h_sq(i2a)  = this%h_sq(i2a) + d2a*n_sq

    else

       fac  = 1.0_DP/(vb-va)

       if (i1a == i1b) then

          d1a  = fac * ( d2b * ( 1.0_DP - 0.5_DP*d2b ) - d2a * ( 1.0_DP - 0.5_DP*d2a ) )
          d2a  = fac * 0.5_DP * ( d2b**2 - d2a**2 )

          this%h1(i1a)    = this%h1(i1a)   + d1a
          this%h(i1a)     = this%h(i1a)    + d1a*n
          this%h_sq(i1a)  = this%h_sq(i1a) + d1a*n_sq

          this%h1(i2a)    = this%h1(i2a)   + d2a
          this%h(i2a)     = this%h(i2a)    + d2a*n
          this%h_sq(i2a)  = this%h_sq(i2a) + d2a*n_sq

       else

          forall(i = i2a:i1b-1)
             this%h1(i)      = this%h1(i)   + 0.5_DP*fac
             this%h(i)       = this%h(i)    + 0.5_DP*fac*n
             this%h_sq(i)    = this%h_sq(i) + 0.5_DP*fac*n_sq

             this%h1(i+1)    = this%h1(i+1)   + 0.5_DP*fac
             this%h(i+1)     = this%h(i+1)    + 0.5_DP*fac*n
             this%h_sq(i+1)  = this%h_sq(i+1) + 0.5_DP*fac*n_sq
          endforall

          d2a  = fac * d1a * ( 1 - 0.5_DP*d1a )
          d1a  = fac * 0.5_DP * d1a**2

          d1b  = fac * d2b * ( 1 - 0.5_DP*d2b )
          d2b  = fac * 0.5_DP * d2b**2

          this%h1(i1a)    = this%h1(i1a)   + d1a
          this%h(i1a)     = this%h(i1a)    + d1a*n
          this%h_sq(i1a)  = this%h_sq(i1a) + d1a*n_sq

          this%h1(i1b)    = this%h1(i1b)   + d1b
          this%h(i1b)     = this%h(i1b)    + d1b*n
          this%h_sq(i1b)  = this%h_sq(i1b) + d1b*n_sq

          this%h1(i2a)    = this%h1(i2a)   + d2a
          this%h(i2a)     = this%h(i2a)    + d2a*n
          this%h_sq(i2a)  = this%h_sq(i2a) + d2a*n_sq

          this%h1(i2b)    = this%h1(i2b)   + d2b
          this%h(i2b)     = this%h(i2b)    + d2b*n
          this%h_sq(i2b)  = this%h_sq(i2b) + d2b*n_sq

       endif

    endif

  endsubroutine histogram1d_add_range

  !% Add a range of values to the histogram (linear interpolation only)
  subroutine histogram1d_add_range_vals(this, valsa, valsb, norm, mask, error)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: valsa(:)
    real(DP), intent(in)              :: valsb(:)
    real(DP), intent(in), optional    :: norm
    logical, intent(in), optional     :: mask(:)
    integer, intent(out), optional    :: error

    ! ---

    integer   :: i, j
    integer   :: i1a(lbound(valsa, 1):ubound(valsa, 1)), i2a(lbound(valsa, 1):ubound(valsa, 1))
    integer   :: i1b(lbound(valsb, 1):ubound(valsb, 1)), i2b(lbound(valsb, 1):ubound(valsb, 1))
    real(DP)  :: d1a(lbound(valsa, 1):ubound(valsa, 1)), d2a(lbound(valsa, 1):ubound(valsa, 1))
    real(DP)  :: d1b(lbound(valsb, 1):ubound(valsb, 1)), d2b(lbound(valsb, 1):ubound(valsb, 1))
    real(DP)  :: va(lbound(valsa, 1):ubound(valsa, 1)), vb(lbound(valsb, 1):ubound(valsb, 1))
    real(DP)  :: fac(lbound(valsa, 1):ubound(valsa, 1))
    real(DP)  :: n, n_sq
    logical   :: m(lbound(valsa, 1):ubound(valsa, 1))

    ! ---

    INIT_ERROR(error)

    if (this%interp /= INTERP_LINEAR) then
       RAISE_ERROR("*histogram1d_add_range* can only be used with linear interpolation.", error)
    endif

    n  = 1.0_DP
    if (present(norm)) then
       n  = norm
    endif
    n_sq  = n**2

    where (valsa > valsb)
       va  = valsb
       vb  = valsa
    elsewhere
       va  = valsa
       vb  = valsb
    endwhere

    fac  = 1.0_DP
    where (va /= vb)
       fac  = 1.0_DP/(vb-va)
    endwhere

    m  = .true.
    if (present(mask)) then
       m  = mask
    endif

    i1a  = int(floor((va-this%min_b)/this%dbin))+1
    i2a  = i1a+1

    d1a  = ( (i2a-1) - (va-this%min_b)/this%dbin )
    d2a  = ( (va-this%min_b)/this%dbin - (i1a-1) )

    i1b  = int(floor((vb-this%min_b)/this%dbin))+1
    i2b  = i1b+1

    d1b  = ( (i2b-1) - (vb-this%min_b)/this%dbin )
    d2b  = ( (vb-this%min_b)/this%dbin - (i1b-1) )

    do j = lbound(valsa, 1), ubound(valsa, 1)

       if (m(j)) then

          if (va(j) == vb(j)) then

             this%h1(i1a(j))    = this%h1(i1a(j))   + d1a(j)
             this%h(i1a(j))     = this%h(i1a(j))    + d1a(j)*n
             this%h_sq(i1a(j))  = this%h_sq(i1a(j)) + d1a(j)*n_sq

             this%h1(i2a(j))    = this%h1(i2a(j))   + d2a(j)
             this%h(i2a(j))     = this%h(i2a(j))    + d2a(j)*n
             this%h_sq(i2a(j))  = this%h_sq(i2a(j)) + d2a(j)*n_sq

          else

             if (i1a(j) == i1b(j)) then

                d1a(j)  = fac(j) * ( d2b(j) * ( 1.0_DP - 0.5_DP*d2b(j) ) - d2a(j) * ( 1.0_DP - 0.5_DP*d2a(j) ) )
                d2a(j)  = fac(j) * 0.5_DP * ( d2b(j)**2 - d2a(j)**2 )

                this%h1(i1a(j))    = this%h1(i1a(j))   + d1a(j)
                this%h(i1a(j))     = this%h(i1a(j))    + d1a(j)*n
                this%h_sq(i1a(j))  = this%h_sq(i1a(j)) + d1a(j)*n_sq

                this%h1(i2a(j))    = this%h1(i2a(j))   + d2a(j)
                this%h(i2a(j))     = this%h(i2a(j))    + d2a(j)*n
                this%h_sq(i2a(j))  = this%h_sq(i2a(j)) + d2a(j)*n_sq

             else

                forall(i = i2a(j):i1b(j)-1)
                   this%h1(i)     = this%h1(i)     + 0.5_DP*fac(j)
                   this%h(i)      = this%h(i)      + 0.5_DP*fac(j)*n
                   this%h_sq(i)   = this%h_sq(i)   + 0.5_DP*fac(j)*n_sq

                   this%h1(i+1)   = this%h1(i+1)   + 0.5_DP*fac(j)
                   this%h(i+1)    = this%h(i+1)    + 0.5_DP*fac(j)*n
                   this%h_sq(i+1) = this%h_sq(i+1) + 0.5_DP*fac(j)*n_sq
                endforall

                d2a(j)  = fac(j) * d1a(j) * ( 1 - 0.5_DP*d1a(j) )
                d1a(j)  = fac(j) * 0.5_DP * d1a(j)**2

                d1b(j)  = fac(j) * d2b(j) * ( 1 - 0.5_DP*d2b(j) )
                d2b(j)  = fac(j) * 0.5_DP * d2b(j)**2

                this%h1(i1a(j))    = this%h1(i1a(j))   + d1a(j)
                this%h(i1a(j))     = this%h(i1a(j))    + d1a(j)*n
                this%h_sq(i1a(j))  = this%h_sq(i1a(j)) + d1a(j)*n_sq

                this%h1(i1b(j))    = this%h1(i1b(j))   + d1b(j)
                this%h(i1b(j))     = this%h(i1b(j))    + d1b(j)*n
                this%h_sq(i1b(j))  = this%h_sq(i1b(j)) + d1b(j)*n_sq

                this%h1(i2a(j))    = this%h1(i2a(j))   + d2a(j)
                this%h(i2a(j))     = this%h(i2a(j))    + d2a(j)*n
                this%h_sq(i2a(j))  = this%h_sq(i2a(j)) + d2a(j)*n_sq

                this%h1(i2b(j))    = this%h1(i2b(j))   + d2b(j)
                this%h(i2b(j))     = this%h(i2b(j))    + d2b(j)*n
                this%h_sq(i2b(j))  = this%h_sq(i2b(j)) + d2b(j)*n_sq

             endif

          endif

       endif

    enddo

  endsubroutine histogram1d_add_range_vals

  !% Add a range of values to the histogram (linear interpolation only)
  subroutine histogram1d_add_range_vals_norms(this, valsa, valsb, norms, mask, error)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: valsa(:)
    real(DP), intent(in)              :: valsb(:)
    real(DP), intent(in)              :: norms(:)
    logical, intent(in), optional     :: mask(:)
    integer, intent(out), optional    :: error

    ! ---

    integer   :: i, j
    integer   :: i1a(lbound(valsa, 1):ubound(valsa, 1)), i2a(lbound(valsa, 1):ubound(valsa, 1))
    integer   :: i1b(lbound(valsb, 1):ubound(valsb, 1)), i2b(lbound(valsb, 1):ubound(valsb, 1))
    real(DP)  :: d1a(lbound(valsa, 1):ubound(valsa, 1)), d2a(lbound(valsa, 1):ubound(valsa, 1))
    real(DP)  :: d1b(lbound(valsb, 1):ubound(valsb, 1)), d2b(lbound(valsb, 1):ubound(valsb, 1))
    real(DP)  :: va(lbound(valsa, 1):ubound(valsa, 1)), vb(lbound(valsb, 1):ubound(valsb, 1))
    real(DP)  :: fac(lbound(valsa, 1):ubound(valsa, 1))
    logical   :: m(lbound(valsa, 1):ubound(valsa, 1))

    ! ---

    INIT_ERROR(error)

    if (this%interp /= INTERP_LINEAR) then
       RAISE_ERROR("*histogram1d_add_range* can only be used with linear interpolation.", error)
    endif

    where (valsa > valsb)
       va  = valsb
       vb  = valsa
    elsewhere
       va  = valsa
       vb  = valsb
    endwhere

    fac  = 1.0_DP
    where (va /= vb)
       fac  = 1.0_DP/(vb-va)
    endwhere

    m  = .true.
    if (present(mask)) then
       m  = mask
    endif

    i1a  = int(floor((va-this%min_b)/this%dbin))+1
    i2a  = i1a+1

    d1a  = ( (i2a-1) - (va-this%min_b)/this%dbin )
    d2a  = ( (va-this%min_b)/this%dbin - (i1a-1) )

    i1b  = int(floor((vb-this%min_b)/this%dbin))+1
    i2b  = i1b+1

    d1b  = ( (i2b-1) - (vb-this%min_b)/this%dbin )
    d2b  = ( (vb-this%min_b)/this%dbin - (i1b-1) )

    do j = lbound(valsa, 1), ubound(valsa, 1)

       if (m(j)) then

          if (va(j) == vb(j)) then

             this%h1(i1a(j))    = this%h1(i1a(j))   + d1a(j)
             this%h(i1a(j))     = this%h(i1a(j))    + d1a(j)*norms(j)
             this%h_sq(i1a(j))  = this%h_sq(i1a(j)) + d1a(j)*norms(j)*norms(j)

             this%h1(i2a(j))    = this%h1(i2a(j))   + d2a(j)
             this%h(i2a(j))     = this%h(i2a(j))    + d2a(j)*norms(j)
             this%h_sq(i2a(j))  = this%h_sq(i2a(j)) + d2a(j)*norms(j)*norms(j)

          else

             if (i1a(j) == i1b(j)) then

                d1a(j)  = fac(j) * ( d2b(j) * ( 1.0_DP - 0.5_DP*d2b(j) ) - d2a(j) * ( 1.0_DP - 0.5_DP*d2a(j) ) )
                d2a(j)  = fac(j) * 0.5_DP * ( d2b(j)**2 - d2a(j)**2 )

                this%h1(i1a(j))    = this%h1(i1a(j))   + d1a(j)
                this%h(i1a(j))     = this%h(i1a(j))    + d1a(j)*norms(j)
                this%h_sq(i1a(j))  = this%h_sq(i1a(j)) + d1a(j)*norms(j)*norms(j)

                this%h1(i2a(j))    = this%h1(i2a(j))   + d2a(j)
                this%h(i2a(j))     = this%h(i2a(j))    + d2a(j)*norms(j)
                this%h_sq(i2a(j))  = this%h_sq(i2a(j)) + d2a(j)*norms(j)*norms(j)

             else

                forall(i = i2a(j):i1b(j)-1)
                   this%h1(i)      = this%h1(i)     + 0.5_DP*fac(j)
                   this%h(i)       = this%h(i)      + 0.5_DP*fac(j)*norms(j)
                   this%h_sq(i)    = this%h_sq(i)   + 0.5_DP*fac(j)*norms(j)*norms(j)

                   this%h1(i+1)    = this%h1(i+1)   + 0.5_DP*fac(j)
                   this%h(i+1)     = this%h(i+1)    + 0.5_DP*fac(j)*norms(j)
                   this%h_sq(i+1)  = this%h_sq(i+1) + 0.5_DP*fac(j)*norms(j)*norms(j)
                endforall

                d2a(j)  = fac(j) * d1a(j) * ( 1 - 0.5_DP*d1a(j) )
                d1a(j)  = fac(j) * 0.5_DP * d1a(j)**2

                d1b(j)  = fac(j) * d2b(j) * ( 1 - 0.5_DP*d2b(j) )
                d2b(j)  = fac(j) * 0.5_DP * d2b(j)**2

                this%h1(i1a(j))    = this%h1(i1a(j))   + d1a(j)
                this%h(i1a(j))     = this%h(i1a(j))    + d1a(j)*norms(j)
                this%h_sq(i1a(j))  = this%h_sq(i1a(j)) + d1a(j)*norms(j)*norms(j)

                this%h1(i1b(j))    = this%h1(i1b(j))   + d1b(j)
                this%h(i1b(j))     = this%h(i1b(j))    + d1b(j)*norms(j)
                this%h_sq(i1b(j))  = this%h_sq(i1b(j)) + d1b(j)*norms(j)*norms(j)

                this%h1(i2a(j))    = this%h1(i2a(j))   + d2a(j)
                this%h(i2a(j))     = this%h(i2a(j))    + d2a(j)*norms(j)
                this%h_sq(i2a(j))  = this%h_sq(i2a(j)) + d2a(j)*norms(j)*norms(j)

                this%h1(i2b(j))    = this%h1(i2b(j))   + d2b(j)
                this%h(i2b(j))     = this%h(i2b(j))    + d2b(j)*norms(j)
                this%h_sq(i2b(j))  = this%h_sq(i2b(j)) + d2b(j)*norms(j)*norms(j)

             endif

          endif

       endif

    enddo

  endsubroutine histogram1d_add_range_vals_norms

  !% Add two histograms
  subroutine histogram1d_add_histogram(this, that, fac)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    type(Histogram1D), intent(in)     :: that
    real(DP), intent(in), optional    :: fac

    ! ---

    this%h1  = this%h1 + that%h1
    if (present(fac)) then
       this%h   = this%h + fac*that%h
    else
       this%h   = this%h + that%h
    endif

  endsubroutine histogram1d_add_histogram

  !% Add a two vectors of multiple histograms histograms
  subroutine histogram1d_add_mult_histograms(this, that, fac)
    implicit none

    type(Histogram1D), intent(inout)  :: this(:)
    type(Histogram1D), intent(in)     :: that(lbound(this,1):ubound(this,1))
    real(DP), intent(in), optional    :: fac

    ! ---

    integer  :: i

    ! ---

    do i = lbound(this,1), ubound(this,1)
       this(i)%h1  = this(i)%h1 + that(i)%h1
       if (present(fac)) then
          this(i)%h   = this(i)%h + fac*that(i)%h
       else
          this(i)%h   = this(i)%h + that(i)%h
       endif
    enddo

  endsubroutine histogram1d_add_mult_histograms

  !% Compute average value in each bin (does only makes sense if
  !% norm != 1 was used in the add functions)
  subroutine histogram1d_average(this, mpi)
    implicit none

    type(Histogram1D),           intent(inout)  :: this
    type(MPI_context), optional, intent(in)     :: mpi

    ! ---

    if (present(mpi)) then
       call sum_in_place(mpi, this%h1)
       call sum_in_place(mpi, this%h)
       call sum_in_place(mpi, this%h_sq)
    endif

    where (this%h1 == 0.0_DP)
       this%h1 = 1.0_DP
    endwhere

    this%h     = this%h / this%h1
    this%h_sq  = this%h_sq / this%h1

  endsubroutine histogram1d_average

  !% Compute Shannon entropy of this histogram. Needs to be called
  !% after normalize.
  elemental function histogram1d_entropy(this) result(val)
    implicit none

    type(Histogram1D), intent(in)  :: this
    real(DP)                       :: val

    ! ---

    integer   :: i
    real(DP)  :: S

    ! ---

    S  = 0.0_DP
    do i = 1, this%n
       if (this%h(i) > 0.0_DP) then
          S  = S - this%h(i)*log(this%h(i)*this%dbin)
       endif
    enddo

    val  = S*this%dbin

  endfunction histogram1d_entropy

  !% Compute expectation value
  elemental function histogram1d_expectation_value(this) result(val)
    implicit none

    type(Histogram1D), intent(in)  :: this
    real(DP)                       :: val

    ! ---

    real(DP)  :: n

    ! ---

    n = sum(this%h)

    if (n == 0.0_DP) then
       val = 0.0_DP
    else
       val  = sum(this%x * this%h / n)
    endif

  endfunction histogram1d_expectation_value

  !% Multiply the histogram by a value - for normalization
  subroutine histogram1d_mul(this, val)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: val

    ! ---

    this%h  = this%h * val

  endsubroutine histogram1d_mul

  !% Multiply multiple histograms by a value - for normalization
  subroutine histograms_mul(this, val)
    implicit none

    type(Histogram1D), intent(inout)  :: this(:)
    real(DP), intent(in)              :: val

    ! ---

    integer  :: i

    ! ---

    do i = lbound(this, 1), ubound(this, 1)
       this(i)%h  = this(i)%h * val
    enddo

  endsubroutine histograms_mul

  !% Multiply multiple histograms by a value - for normalization
  subroutine histograms2_mul(this, val)
    implicit none

    type(Histogram1D), intent(inout)  :: this(:, :)
    real(DP), intent(in)              :: val

    ! ---

    integer  :: i, j

    ! ---

    do i = lbound(this, 1), ubound(this, 1)
       do j = lbound(this, 2), ubound(this, 2)
          this(i, j)%h  = this(i, j)%h * val
       enddo
    enddo

  endsubroutine histograms2_mul

  !% Multiply the histogram by multiple value - for normalization
  subroutine histogram1d_mul_vals(this, vals)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: vals(this%n)

    ! ---

    this%h  = this%h * vals

  endsubroutine histogram1d_mul_vals

  !% Divide histogram by a value
  subroutine histogram1d_div(this, val)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: val

    ! ---

    this%h  = this%h / val

  endsubroutine histogram1d_div

  !% Divide multiple histograms by a value
  subroutine histograms_div(this, val)
    implicit none

    type(Histogram1D), intent(inout)  :: this(:)
    real(DP), intent(in)              :: val

    ! ---

    integer  :: i

    ! ---

    do i = lbound(this, 1), ubound(this, 1)
       this(i)%h  = this(i)%h / val
    enddo

  endsubroutine histograms_div

  !% Divide multiple histograms by a value
  subroutine histograms2_div(this, val)
    implicit none

    type(Histogram1D), intent(inout)  :: this(:, :)
    real(DP), intent(in)              :: val

    ! ---

    integer  :: i, j

    ! ---

    do i = lbound(this, 1), ubound(this, 1)
       do j = lbound(this, 2), ubound(this, 2)
          this(i, j)%h  = this(i, j)%h / val
       enddo
    enddo

  endsubroutine histograms2_div

  !% Divide histogram by multiple values
  subroutine histogram1d_div_vals(this, vals)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: vals(this%n)

    ! ---

    this%h  = this%h / vals

  endsubroutine histogram1d_div_vals

  !% Normalize histogram
  elemental subroutine histogram1d_normalize(this)
    implicit none

    type(Histogram1D), intent(inout)  :: this

    ! ---

    real(DP)  :: n

    ! ---

    n = sum(this%h*this%dbin)

    if (n /= 0.0_DP) then
       this%h  = this%h / n
    endif

  endsubroutine histogram1d_normalize

  !% OpenMP reduction - *thpriv* is a threadprivate histogram,
  !% this needs to be called within an *omp parallel* construct.
  subroutine histogram1d_reduce(this, thpriv)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    type(Histogram1D), intent(in)     :: thpriv

    ! ---

    !$omp single
    call initialise(this, thpriv)
    !$omp end single

    !$omp critical
    this%h     = this%h    + thpriv%h
    this%h_sq  = this%h_sq + thpriv%h_sq
    this%h1    = this%h1   + thpriv%h1
    !$omp end critical

  endsubroutine histogram1d_reduce

  !% Smooth histogram by convolution with a Gaussian
  subroutine histogram1d_smooth(this, sigma)
    implicit none

    type(Histogram1D), intent(inout)  :: this
    real(DP), intent(in)              :: sigma

    ! ---

    integer   :: i
    real(DP)  :: fac, sigma_sq, h(this%n)

    ! ---

    h         = 0.0_DP

    fac       = this%dbin/(sqrt(2*PI)*sigma)
    sigma_sq  = sigma**2

    do i = 1, this%n
       h  = h + fac*this%h(i)*exp(-((this%x - this%x(i))**2/(2*sigma_sq)))
    enddo

    this%h  = h

  endsubroutine histogram1d_smooth

  !% Sum histogram from different processors onto root
  subroutine histogram1d_sum_in_place(mpi, this)
    implicit none

    type(MPI_context), intent(in)     :: mpi
    type(Histogram1D), intent(inout)  :: this

    ! ---

    call sum_in_place(mpi, this%h)
    call sum_in_place(mpi, this%h_sq)
    call sum_in_place(mpi, this%h1)

  endsubroutine histogram1d_sum_in_place

#if 0
  !**********************************************************************
  ! Memory estimate logging
  !**********************************************************************
  subroutine log_memory_estimate_histogram(this)
    implicit none

    type(Histogram1D), intent(in)  :: this

    ! ---

    call log_memory_estimate(this%x)
    call log_memory_estimate(this%h)
    call log_memory_estimate(this%h_sq)
    call log_memory_estimate(this%h1)

  endsubroutine log_memory_estimate_histogram

  !**********************************************************************
  ! Memory estimate logging
  !**********************************************************************
  subroutine log_memory_estimate_histogram2(this)
    implicit none

    type(Histogram1D), intent(in)  :: this(:)

    ! ---

    integer  :: i

    ! ---

    do i = lbound(this, 1), ubound(this, 1)
       call log_memory_estimate_histogram(this(i))
    enddo

  endsubroutine log_memory_estimate_histogram2

  !**********************************************************************
  ! Memory estimate logging
  !**********************************************************************
  subroutine log_memory_estimate_histogram3(this)
    implicit none

    type(Histogram1D), intent(in)  :: this(:, :)

    ! ---

    integer  :: i, j

    ! ---

    do j = lbound(this, 2), ubound(this, 2)
       do i = lbound(this, 1), ubound(this, 1)
          call log_memory_estimate_histogram(this(i, j))
       enddo
    enddo

  endsubroutine log_memory_estimate_histogram3
#endif

endmodule histogram1d_module

max_stars_repo_name: albapa/QUIP

max_issues_repo_name: albapa/QUIP

max_forks_repo_name: albapa/QUIP


text:
C     LAST UPDATE 24/02/94
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE FLIST(KUNIT)
      IMPLICIT NONE
C
C Purpose: Lists coordinates and intensities held in common and writes
C          to file FIX.OUT.
C
C Calls   RDCOM:
C Called by: FIX
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C FIX include file:
C
      INCLUDE 'FIXPAR.COM'
C
C Arguments:
C
      INTEGER KUNIT
C
C Local variables:
C
      INTEGER I,J,NTYPE
      REAL RTOD,D
      CHARACTER*7 TYPE(3)
C
      DATA TYPE /'POINT', 'POLYGON', 'SECTOR'/
      DATA RTOD /57.29578/
C
C-----------------------------------------------------------------------
      WRITE(KUNIT,1005)
C
C========List known parameters
C
      IF(GOTWAV)THEN
         WRITE(KUNIT,1000)WAVE
      ELSE
         WRITE(KUNIT,1010)
      ENDIF
      IF(GOTSDD)THEN
         WRITE(KUNIT,1020)SDD
      ELSE
         WRITE(KUNIT,1030)
      ENDIF
      IF(GOTCEN)THEN
         WRITE(KUNIT,1040)XC,YC
      ELSE
         WRITE(KUNIT,1050)
      ENDIF
      IF(GOTROT)THEN
         WRITE(KUNIT,1060)RTOD*ROTX,RTOD*ROTY,RTOD*ROTZ
      ELSE
         WRITE(KUNIT,1070)
      ENDIF
      IF(GOTTIL)THEN
         WRITE(KUNIT,1080)-RTOD*TILT
      ELSE
         WRITE(KUNIT,1090)
      ENDIF
C
C========List points
C
      IF(NPTS.EQ.0)RETURN
      WRITE(KUNIT,1100)
      DO 10 I=1,NPTS
         IF(NVALS(I).LT.0)THEN
            NTYPE = 3
         ELSEIF(NVALS(I).GT.1)THEN
            NTYPE = 2
         ELSE
            NTYPE = 1
         ENDIF
         D = SQRT(RZ(1,I)*RZ(1,I)+RZ(2,I)*RZ(2,I))
         WRITE(KUNIT,1110)TYPE(NTYPE),I,(FXY(J,I),J=1,2),
     &                    (SXY(J,I),J=1,2),(RZ(J,I),J=1,2),D,AVPIX(I)
 10   CONTINUE
C
 1005 FORMAT(/1X,'New listing...'/)
 1000 FORMAT(1X,'Wavelength = ',F6.3)
 1010 FORMAT(1X,'Don''t know wavelength')
 1020 FORMAT(1X,'Specimen-to-detector distance = ',F8.1)
 1030 FORMAT(1X,'Don''t know specimen-to-detector distance')
 1040 FORMAT(1X,'Centre coordinates = ',2F8.1)
 1050 FORMAT(1X,'Don''t know centre coordinates')
 1060 FORMAT(1X,'Detector rotation = ',F9.4/
     &       1X,'         twist    = ',F9.4/
     &       1X,'         tilt     = ',F9.4)
 1070 FORMAT(1X,'Don''t know detector rotation angles')
 1080 FORMAT(1X,'Specimen tilt angle = ',F6.1)
 1090 FORMAT(1X,'Don''t know specimen tilt angle')
 1100 FORMAT(/2X,'No.',6X,'Filespace',8X,'Standard',11X,'Reciprocal',
     &       6X,'Av. pixel',3X,'Background',3X,'Intensity',6X,'Type'/
     &       2X,'         X       Y',8X,'X       Y',10X,'R        Z')
 1110 FORMAT(1X,A7,1X,I3,1X,2F8.1,1X,2F8.1,4X,3F9.5,
     &       3X,G12.5)
      END

ext: f


text:
      subroutine open_topdrawer_file
      implicit none
      logical useitmax
      common/cuseitmax/useitmax
      open(unit=99,file='MADatNLO.top',status='unknown')
      useitmax=.false.
      return
      end

      subroutine close_topdrawer_file
      implicit none
      close(99)
      return
      end

      subroutine open_root_file
      implicit none
      logical useitmax
      common/cuseitmax/useitmax
      call rinit('MADatNLO.root')
      useitmax=.true.
      return
      end

      subroutine close_root_file
      implicit none
      call rwrit()
      call rclos()
      return
      end

      subroutine HwU_write_file
      implicit none
      double precision xnorm
      open (unit=99,file='MADatNLO.HwU',status='unknown')
      xnorm=1d0
      call HwU_output(99,xnorm)
      close (99)
      return
      end

max_stars_repo_name: JaySandesara/madminer

max_issues_repo_name: JaySandesara/madminer

max_forks_repo_name: JaySandesara/madminer


text:
C
C	$Id: gflas4.f,v 1.8 2008-07-27 00:17:13 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE GFLAS4(ID,FNAME)
C
C_BEGIN PROLOGUE GFLAS4
C
C_PURPOSE
C    A call to this subroutine does the necessary bookkeeping
C    to allow a dataset of plotting instructions created in a previous
C    job step with a GFLAS1/GFLAS2 sequence to be used in the
C    current program with a GFLAS3 call.  GFLAS4 can be used
C    only with the NCAR GKS package.
C
C
C_DESCRIPTION
C    A call to this subroutine does the necessary bookkeeping
C    to allow a dataset of plotting instructions created in a previous
C    job step with a GFLAS1/GFLAS2 sequence to be used in the
C    current program with a GFLAS3 call.  This allows for inserting
C    previously-created plotting instructions into the output
C    pipeline at any user-selected point.  GFLAS4 can be used only
C    with the NCAR GKS package.
C
C_ARGUMENTS
C
C ID      (INPUT,SINGLE,VARIABLE)  An integer variable in the range
C         [0,99] which is used as a user-defined identifier for the
C         dataset of graphics instructions contained in FNAME (see
C         below.)  ID can be used as an argument to subsequent
C         GFLAS3 calls.  ID cannot be equal to an identifier
C         used as an argument to GFLAS1 in the current program.
C FNAME   (INPUT,SINGLE,VARIABLE)  A character string containing
C         the dataset name of the dataset which contains plotting
C         instructions generated previously with a GFLAS1/GFLAS2
C         sequence.  FNAME can be at most seven characters in length.
C
C_I/O
C    (None)
C
C_ROUTINES CALLED
C    PLOTIT, SETER
C
C_COMMON BLOCKS
C    GFLASH
C
C_MAXIMUM GKS LEVEL
C    0A
C
C_LANGUAGE
C    FORTRAN
C
C_REVISION DATE
C    890117 (YYMMDD)
C
C_HISTORY
C    This subroutine was coded at NCAR by Fred Clare.  Its purpose
C    is to duplicate the functionality of the previous NCAR System
C    Plot Package entry FLASH4 in a GKS environment.
C
C_REFERENCES
C    "The System Plot Package", NCAR Technical Note number
C    NCAR-TN/162+IA, May, 1983.
C
C    "NCAR Graphics User's Guide", Version 2.00, NCAR Technical Note
C    number NCAR/TN-283+IA, August, 1987.
C
C_LONG DESCRIPTION
C    GFLAS4 finds the number of records in FNAME and does the
C    necessary bookkeeping to establish the dataset as valid for
C    subsequent use by GFLAS3.  GFLAS4 utilizes a GKS ESCAPE
C    function unique to the NCAR GKS package.  The function
C    identifier for this ESCAPE function is -1392.  If one tries
C    to use GFLAS4 with a non-NCAR GKS package, he will get an
C    error to the effect that no ESCAPE function exists with
C    function idendifier -1392.
C
C_EXAMPLES
C    Example 1 --  Establish a previously existing dataset of plotting
C                  instructions -- MF09 -- for use with GFLAS3.
C          SUBROUTINE INSDAT
C          CALL GFLAS4(9,'MF09')
C          CALL GFLAS3(9)
C          RETURN
C          END
C
C_END PROLOGUE GFLAS4
C
C-------------------------------------------------------------
C
C  Common for FLASH package.
C
      COMMON /GFLASH/MODEF,IOPWKS(100),IOACT(100),NUMOP,IWISSI
C
      CHARACTER*(*) FNAME
      CHARACTER*80 IDR(2),ODR(1)
C
      IF (MODEF .EQ. 1) THEN
        GO TO 106
      ELSE
        MODEF = 4
      ENDIF
C
C  Check if GKS is open.
C
      CALL GQOPS(IOPSTA)
      IF (IOPSTA .LT. 1) GO TO 105
C
C  Make sure that this GKS implementation has WISS (level 2).
C
      CALL GQEWK(1,IER,NUMB,NTYP)
      DO 10 I=1,NUMB
        CALL GQEWK(I,IER,NUMB,NTYP)
        CALL GQWKCA(NTYP,IER,NCAT)
        IF (NCAT .EQ. 3) GO TO 20
   10 CONTINUE
      GO TO 102
   20 CONTINUE
C
C  Make sure that WISS is open (IWISSI is initialized to an
C  arbitrary illegal value).
C
      IWISSI = -1335
      CALL GQOPWK(1,IER,IOL,IWKID)
      NUMOP = IOL
      IF (NUMOP .LT. 1) GO TO 103
      IF (NUMOP .GT. 100) GO TO 104
      DO 30 I=1,NUMOP
      CALL GQOPWK(I,IER,IOL,IOPWKS(I))
      IOACT(I) = 0
      CALL GQWKS(IOPWKS(I),IER,IOACT(I))
      CALL GQWKC(IOPWKS(I),IER,ICON,ITYPE)
      CALL GQWKCA(ITYPE,IER,ICAT)
      IF (ICAT .EQ. 3) THEN
        IWISSI = IOPWKS(I)
        CALL GQWKC(IWISSI,IER,ICONW,ITYPE)
      ENDIF
   30 CONTINUE
      IF (IWISSI .EQ. -1335) GO TO 103
C
C  Encode parameters into input data record.
C
      IDR(1) = ' '
      IDR(2) = ' '
      IFID = -1392
      ILEN = LEN(FNAME)
      WRITE(IDR(1)( 1:10),500) IFID
      WRITE(IDR(1)(11:20),500) ID
      WRITE(IDR(1)(21:23),501) ICONW
      IF (ILEN .LE. 57) THEN
        IDR(1)(24:24+ILEN-1) = FNAME(1:ILEN)
      ELSE IF (ILEN .LE. 137) THEN
        IDR(1)(24:24+ILEN-1) = FNAME(1:57)
        IDR(2)(1:ILEN-57) = FNAME(58:ILEN)
      ELSE
        GO TO 107
      ENDIF
      CALL PLOTIT(0,0,2)
      IF (ILEN .LE. 57) THEN
        CALL GESC(IFID,1,IDR,1,IDUM,ODR)
      ELSE
        CALL GESC(IFID,2,IDR,1,IDUM,ODR)
      ENDIF
      CALL PLOTIT(0,0,2)
      RETURN
C
  102 CONTINUE
      CALL GECLKS
      CALL SETER(' GFLAS4 -- THE NON NCAR GKS PACKAGE IS NOT LEVEL 2',
     -           10,2)
      STOP
  103 CONTINUE
      CALL GECLKS
      CALL SETER(' GFLAS4 -- WISS MUST BE OPEN BEFORE CALL TO GFLAS4',
     -           11,2)
      STOP
  104 CONTINUE
      CALL GECLKS
      CALL SETER(' GFLAS4 -- MAXIMUM NUMBER OF OPEN WORKSTATIONS ALLOWED
     - IS 100',12,2)
      STOP
  105 CONTINUE
      CALL GECLKS
      CALL SETER(' GFLAS4 -- GKS MUST BE OPEN BEFORE A CALL TO GFLAS4',
     -  13,2)
  106 CONTINUE
      CALL GECLKS
      CALL SETER(' GFLAS4 -- CANNOT BE CALLED WHILE GFLAS1 IS STILL OPEN
     -',14,2)
  107 CONTINUE
      CALL GECLKS
      CALL SETER(' GFLAS4 -- Maximum length of a file name is 137 charac
     -ters',15,2)
C
  500 FORMAT(I10)
  501 FORMAT(I3)
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


! { dg-do run }
! Fixed by the patch for PRs 60357 and 61275
!
! Contributed by Stefan Mauerberger  <stefan.mauerberger@gmail.com>
!
PROGRAM main
  IMPLICIT NONE
  TYPE :: test_typ
    REAL, ALLOCATABLE :: a
  END TYPE
  TYPE(test_typ) :: my_test_typ
  my_test_typ = test_typ (a = 1.0)
  if (abs (my_test_typ%a - 1.0) .gt. 1e-6) STOP 1
END PROGRAM main


text:
//////////
/ n1/i386/tables/prefac.f
/ Code table preface.
/ Useful abbreviations for code tables.
/ Must be the first file in the tabgen input list.
/ i386.
////////

#define	BYTE	(FS8|FU8)			/* [u]char		*/
#define	WORD	(FS16|FU16)			/* [u]short		*/
#define	SHORT	(BYTE|WORD)			/* [u]char or [u]short	*/
#define	LONG	(FS32|FU32)			/* [u]int or [u]long	*/
#define	PTX	(FPTR|FPTB)			/* pointer		*/
#define	DWORD	(LONG|PTX)			/* 32-bit incl. pointer	*/
#define	SINT	(FS32|FS16|FS8)			/* signed integers	*/
#define	UINT	(FU32|FU16|FU8)			/* unsigned integers	*/
#define	FLD	(FFLD32|FFLD16|FFLD8)		/* bitfields		*/
#define	FLOAT	(FF32|FF64)			/* float or double	*/
#define	NFLT	(BYTE|WORD|DWORD)		/* integer or pointer	*/
#define	NDPARG	(FLOAT|FS16|FS32)		/* 80x87 args		*/
#define	PEREL	(PEQ|PNE)			/* equality relations	*/
#define	PSREL	(PEQ|PNE|PGT|PGE|PLT|PLE)	/* signed relations	*/
#define	PREL	(PEQ|PNE|PGT|PGE|PLT|PLE|PUGT|PUGE|PULT|PULE)	/* relations */
#define	PVALUE	(PRVALUE|PLVALUE)		/* value contexts	*/

//////////
/ end of n1/i386/tables/prefac.f
//////////

max_stars_repo_name: gspu/Coherent

max_issues_repo_name: a-mail-group/coherent

max_forks_repo_name: a-mail-group/coherent


text:
! RUN: bbc %s -o "-" | FileCheck %s

! Test real intrinsic operation lowering to FIR.

! CHECK-LABEL:eq0_test
LOGICAL FUNCTION eq0_test(x0, x1)
REAL(4) :: x0
REAL(4) :: x1
! CHECK-DAG:[[reg1:%[0-9]+]] = fir.load %arg0
! CHECK-DAG:[[reg2:%[0-9]+]] = fir.load %arg1
! CHECK:[[reg3:%[0-9]+]] = arith.cmpf oeq, [[reg1]], [[reg2]] : f32
! CHECK:fir.convert [[reg3]] {{.*}} -> !fir.logical<4>
eq0_test = x0 .EQ. x1
END FUNCTION

! CHECK-LABEL:ne1_test
LOGICAL FUNCTION ne1_test(x0, x1)
REAL(4) :: x0
REAL(4) :: x1
! CHECK-DAG:[[reg1:%[0-9]+]] = fir.load %arg0
! CHECK-DAG:[[reg2:%[0-9]+]] = fir.load %arg1
! CHECK:[[reg3:%[0-9]+]] = arith.cmpf une, [[reg1]], [[reg2]] : f32
! CHECK:fir.convert [[reg3]] {{.*}} -> !fir.logical<4>
ne1_test = x0 .NE. x1
END FUNCTION

! CHECK-LABEL:lt2_test
LOGICAL FUNCTION lt2_test(x0, x1)
REAL(4) :: x0
REAL(4) :: x1
! CHECK-DAG:[[reg1:%[0-9]+]] = fir.load %arg0
! CHECK-DAG:[[reg2:%[0-9]+]] = fir.load %arg1
! CHECK:[[reg3:%[0-9]+]] = arith.cmpf olt, [[reg1]], [[reg2]] : f32
! CHECK:fir.convert [[reg3]] {{.*}} -> !fir.logical<4>
lt2_test = x0 .LT. x1
END FUNCTION

! CHECK-LABEL:le3_test
LOGICAL FUNCTION le3_test(x0, x1)
REAL(4) :: x0
REAL(4) :: x1
! CHECK-DAG:[[reg1:%[0-9]+]] = fir.load %arg0
! CHECK-DAG:[[reg2:%[0-9]+]] = fir.load %arg1
! CHECK:[[reg3:%[0-9]+]] = arith.cmpf ole, [[reg1]], [[reg2]] : f32
! CHECK:fir.convert [[reg3]] {{.*}} -> !fir.logical<4>
le3_test = x0 .LE. x1
END FUNCTION

! CHECK-LABEL:gt4_test
LOGICAL FUNCTION gt4_test(x0, x1)
REAL(4) :: x0
REAL(4) :: x1
! CHECK-DAG:[[reg1:%[0-9]+]] = fir.load %arg0
! CHECK-DAG:[[reg2:%[0-9]+]] = fir.load %arg1
! CHECK:[[reg3:%[0-9]+]] = arith.cmpf ogt, [[reg1]], [[reg2]] : f32
! CHECK:fir.convert [[reg3]] {{.*}} -> !fir.logical<4>
gt4_test = x0 .GT. x1
END FUNCTION

! CHECK-LABEL:ge5_test
LOGICAL FUNCTION ge5_test(x0, x1)
REAL(4) :: x0
REAL(4) :: x1
! CHECK-DAG:[[reg1:%[0-9]+]] = fir.load %arg0
! CHECK-DAG:[[reg2:%[0-9]+]] = fir.load %arg1
! CHECK:[[reg3:%[0-9]+]] = arith.cmpf oge, [[reg1]], [[reg2]] : f32
! CHECK:fir.convert [[reg3]] {{.*}} -> !fir.logical<4>
ge5_test = x0 .GE. x1
END FUNCTION

! CHECK-LABEL:add6_test
REAL(4) FUNCTION add6_test(x0, x1)
REAL(4) :: x0
REAL(4) :: x1
! CHECK-DAG:[[reg1:%[0-9]+]] = fir.load %arg0
! CHECK-DAG:[[reg2:%[0-9]+]] = fir.load %arg1
! CHECK:addf [[reg1]], [[reg2]] : f32
add6_test = x0 + x1
END FUNCTION

! CHECK-LABEL:sub7_test
REAL(4) FUNCTION sub7_test(x0, x1)
REAL(4) :: x0
REAL(4) :: x1
! CHECK-DAG:[[reg1:%[0-9]+]] = fir.load %arg0
! CHECK-DAG:[[reg2:%[0-9]+]] = fir.load %arg1
! CHECK:subf [[reg1]], [[reg2]] : f32
sub7_test = x0 - x1
END FUNCTION

! CHECK-LABEL:mult8_test
REAL(4) FUNCTION mult8_test(x0, x1)
REAL(4) :: x0
REAL(4) :: x1
! CHECK-DAG:[[reg1:%[0-9]+]] = fir.load %arg0
! CHECK-DAG:[[reg2:%[0-9]+]] = fir.load %arg1
! CHECK:mulf [[reg1]], [[reg2]] : f32
mult8_test = x0 * x1
END FUNCTION

! CHECK-LABEL:div9_test
REAL(4) FUNCTION div9_test(x0, x1)
REAL(4) :: x0
REAL(4) :: x1
! CHECK-DAG:[[reg1:%[0-9]+]] = fir.load %arg0
! CHECK-DAG:[[reg2:%[0-9]+]] = fir.load %arg1
! CHECK:divf [[reg1]], [[reg2]] : f32
div9_test = x0 / x1
END FUNCTION

max_stars_repo_name: lxbndr/llvm-project

max_issues_repo_name: lxbndr/llvm-project

max_forks_repo_name: lxbndr/llvm-project


text:
C
C      $Id: xy09f.f,v 1.6 2010-03-15 22:49:25 haley Exp $
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C                Copyright (C)  1995                                   C
C        University Corporation for Atmospheric Research               C
C                All Rights Reserved                                   C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  File:       xy09f.f
C
C  Author:     Mary Haley
C          National Center for Atmospheric Research
C          PO 3000, Boulder, Colorado
C
C  Date:       Thu Feb  9 10:00:38 MST 1995
C
C  Description:    This example is similar to the ncargex Autograph
C                  example "agex06".  It shows how to create different
C                  kinds of axes using a combination of TickMark and
C                  XyPlot resources.
C
      external NhlFAppClass
      external NhlFCairoWindowWorkstationClass
      external NhlFNcgmWorkstationClass
      external NhlFPSWorkstationClass
      external NhlFPDFWorkstationClass
      external NhlFCairoPSPDFWorkstationClass
      external NhlFCairoImageWorkstationClass
      external NhlFCoordArraysClass
      external NhlFXyPlotClass
C
C Define the number of points in each curve and the number of colors.
C
      parameter(NPTS=501,NCOLORS=6)
      parameter(PI100=.031415926535898)

      integer appid,xworkid,plotid,dataid
      integer rlist, i
      character*8 plot_name
C
C Create data arrays for XyPlot.
C
      real xdra(NPTS), ydra(NPTS)
C
C Modify the color map.  Color indices '1' and '2' are the background
C and foreground colors respectively.
C
      integer len(2)
      data len/3,NCOLORS/
      real cmap(3,NCOLORS)
      data cmap/0.0,0.0,0.0,
     +          1.0,1.0,1.0,
     +          1.0,0.5,0.0,
     +          0.0,1.0,0.5,
     +          0.5,0.0,1.0,
     +          0.6,0.2,0.2/

      CHARACTER*7  wks_type
C
C Define the workstation type
C
      wks_type = "x11"
C
C
C Initialize the HLU library and set up resource template.
C
      call NhlFInitialize
      call NhlFRLCreate(rlist,'setrl')
C
C Create Application object.  The Application object name is used to
C determine the name of the resource file, which is "xy09.res" in
C this case.
C
      call NhlFRLClear(rlist)
      call NhlFRLSetString(rlist,'appDefaultParent','True',ierr)
      call NhlFRLSetString(rlist,'appUsrDir','./',ierr)
      call NhlFCreate(appid,'xy09',NhlFAppClass,0,rlist,ierr)

      if (wks_type.eq."ncgm".or.wks_type.eq."NCGM") then
C
C Create an NCGM workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkMetaName','./xy09f.ncgm',ierr)
         call NhlFRLSetMDFloatArray(rlist,'wkColorMap',cmap,2,len,ierr)
         call NhlFCreate(xworkid,'xy09Work',
     +        NhlFNcgmWorkstationClass,0,rlist,ierr)
      else if (wks_type.eq."x11".or.wks_type.eq."X11") then
C
C Create an XWorkstation object.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkPause','True',ierr)
         call NhlFRLSetMDFloatArray(rlist,'wkColorMap',cmap,2,len,ierr)
         call NhlFCreate(xworkid,'xy09Work',
     +                NhlFCairoWindowWorkstationClass,
     +                0,rlist,ierr)
      else if (wks_type.eq."oldps".or.wks_type.eq."OLDPS") then
C
C Create an older-style PS workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkPSFileName','./xy09f.ps',ierr)
         call NhlFRLSetMDFloatArray(rlist,'wkColorMap',cmap,2,len,ierr)
         call NhlFCreate(xworkid,'xy09Work',
     +        NhlFPSWorkstationClass,0,rlist,ierr)
      else if (wks_type.eq."oldpdf".or.wks_type.eq."OLDPDF") then
C
C Create an older-style PDF workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkPDFFileName','./xy09f.pdf',ierr)
         call NhlFRLSetMDFloatArray(rlist,'wkColorMap',cmap,2,len,ierr)
         call NhlFCreate(xworkid,'xy09Work',
     +        NhlFPDFWorkstationClass,0,rlist,ierr)
      else if (wks_type.eq."pdf".or.wks_type.eq."PDF".or.
     +         wks_type.eq."ps".or.wks_type.eq."PS") then
C
C Create a cairo PS/PDF workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkFileName','./xy09f',ierr)
         call NhlFRLSetString(rlist,'wkFormat',wks_type,ierr)
         call NhlFRLSetMDFloatArray(rlist,'wkColorMap',cmap,2,len,ierr)
         call NhlFCreate(xworkid,'xy09Work',
     +        NhlFcairoPSPDFWorkstationClass,0,rlist,ierr)
      else if (wks_type.eq."png".or.wks_type.eq."PNG") then
C
C Create a cairo PNG workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkFileName','./xy09f',ierr)
         call NhlFRLSetString(rlist,'wkFormat',wks_type,ierr)
         call NhlFRLSetMDFloatArray(rlist,'wkColorMap',cmap,2,len,ierr)
         call NhlFCreate(xworkid,'xy09Work',
     +        NhlFcairoImageWorkstationClass,0,rlist,ierr)
      endif
C
C Initialize data for the XyPlot object.
C
      do 10 i = 1,NPTS
         theta = PI100*real(i-1)
         xdra(i)=500.+.9*real(i-1)*cos(theta)
         ydra(i)=500.+.9*real(i-1)*sin(theta)
 10   continue
C
C Define the data object.
C
      call NhlFRLClear(rlist)
      call NhlFRLSetFloatArray(rlist,'caXArray',xdra,NPTS,ierr)
      call NhlFRLSetFloatArray(rlist,'caYArray',ydra,NPTS,ierr)
      call NhlFCreate(dataid,'xyData',NhlFCoordArraysClass,
     +                0,rlist,ierr)
C
C Create and draw four XyPlot objects.
C
      do 30 i = 1,4
         call NhlFRLClear(rlist)
         call NhlFRLSetInteger(rlist,'xyCoordData',dataid,ierr)
         write(plot_name,40)i
         call NhlFCreate(plotid,plot_name,NhlFXyPlotClass,xworkid,
     +        rlist,ierr)
         call NhlFDraw(plotid,ierr)
 30   continue
 40   format('xyPlot',i1)

      call NhlFFrame(xworkid,ierr)
C
C NhlDestroy destroys the given id and all of its children
C so destroying "xworkid" will also destroy "plotid".
C
      call NhlFRLDestroy(rlist)
      call NhlFDestroy(xworkid,ierr)
      call NhlFDestroy(appid,ierr)
C
C Restores state.
C
      call NhlFClose
      stop
      end


SUBROUTINE REVINM(EVIVAR, FILNAM, FULLPATH)
    !
    !_begin_doc
    !  RCS ID string
    !  $Id: revinm.F,v 1.1.1.1 1997/03/28 20:58:33 marki Exp $
    !
    !  Documentation for subroutine REVINM
    !
    !  Purpose:  Given the enviroment variable EVIVAR and the filename FILNAM
    !  --------  returns the full path FULLPATH to the file.
    !
    !  Calling Sequence:
    !  ----------------
    !
    !  Input Parameters:  (EVIVAR - C - Name of enviroment variable)
    !  ----------------    FILNAM - C - Name of file
    !
    !  Output Parameters:  (FULLPATH - C - EVIVAR/FILNAM )
    !  -----------------
    !
    !  Called from:
    !  ------------
    !
    !  Other routines:
    !  ---------------
    !
    !  Notes:
    !  ------
    !
    !  Author:   Arne Freyberger      Created:  Wed Mar 15 13:39:27 EST 1995
    !  -------
    !
    !  Major revisions:
    !  ----------------
    !  4/4/95   Removed messaging so that code can be used with the
    !  apf      bulk of the RECUTL/RECSIS package.
    !_end_doc
    !
    IMPLICIT NONE
    SAVE
    !
    !_begin_inc
    !  include files :
    !  ---------------------
    ! BOS common block  uncomment the next line for BOS include file
    !#include "bcs.inc"
    !_end_inc
    !
    !_begin_var
    !  input/output variables:
    !  -----------------------
    !
    !  Local pre-defined variables:
    !  ---------------------------
    !  RCS information:
    CHARACTER*132  CFILE, CREVIS, CSTATE, CDATE, CAUTHO, CRCSID
    PARAMETER (CFILE = '$RCSfile: revinm.F,v $')
    PARAMETER (CREVIS = '$Revision: 1.1.1.1 $')
    PARAMETER (CSTATE = '$State: Exp $')
    PARAMETER (CDATE = '$Date: 1997/03/28 20:58:33 $')
    PARAMETER (CAUTHO = '$Author: marki $')
    DATA CRCSID/&
            '$Id: revinm.F,v 1.1.1.1 1997/03/28 20:58:33 marki Exp $'&
            /
    !  Module information:
    CHARACTER*(*)  CRNAME, CRAUTH
    CHARACTER*100  CRMESS
    PARAMETER (CRNAME = 'REVINM')
    PARAMETER (CRAUTH = 'Arne Freyberger')
    !
    !  Local User defined variables:
    !  -----------------------------
    INTEGER LENOCC
    EXTERNAL LENOCC
    CHARACTER*(*) EVIVAR, FILNAM, FULLPATH
    INTEGER ICEND, ICBEG, IADD, ICEND2
    !_end_var
    !
    !  executable code for routine REVINM:
    !  -------------------------------------
    !
    !
    CALL GETENV(EVIVAR, FULLPATH)
    ICEND = LENOCC(FULLPATH)
    IF (FULLPATH(1:1) .EQ. ' ') THEN
        FULLPATH = FILNAM
        ICEND2 = LENOCC(EVIVAR)
        WRITE(CRMESS, 10)EVIVAR(1:ICEND2)
        !        CALL RECLOG(CRNAME,'W',CRMESS)
        RETURN
    ENDIF
    ICEND = ICEND + 1
    !
    FULLPATH(ICEND:ICEND) = '/'
    !
    IADD = LENOCC(FILNAM)
    ICBEG = ICEND + 1
    ICEND = ICBEG + IADD
    FULLPATH(ICBEG:ICEND) = FILNAM(1:IADD)
    !
    10    FORMAT('Enviroment variable:  ', A, '  not defined')
    RETURN
END
!
!------------------------------------------------------------------------------


text:
great_rsd(aa1,x1).
great_rsd(dd1,x1).
great_rsd(cc1,u1).
great_rsd(p1,cc1).
great_rsd(ee1,x1).
great_rsd(v1,x1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
!===== T. Dennen, May 2014
!===== Initialise array of triangle integral values.
!===== For use with qqb->4gamma and qqb->2j2gamma
!===== coefficients by same author.
      subroutine Triint_init(i1,i2,i3,i4,i5,i6,
     & Triint,ord)
      implicit none
      include 'constants.f'
      include 'sprods_com.f'
      include 'scale.f'
      integer i1,i2,i3,i4,i5,i6
      integer ord
      double complex Triint(90), qlI3
      double precision t

      t(i1,i2,i3) = s(i1,i2) + s(i2,i3) + s(i3,i1)

      Triint(:) = czip

      Triint(1) = qlI3(0d0,0d0,s(i1,i2),0d0,0d0,0d0,musq,ord)
      Triint(2) = qlI3(0d0,s(i1,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(3) = qlI3(s(i5,i6),0d0,0d0,0d0,0d0,0d0,musq,ord)
      Triint(4) = qlI3(0d0,s(i1,i5),0d0,0d0,0d0,0d0,musq,ord)
      Triint(5) = qlI3(s(i4,i6),0d0,0d0,0d0,0d0,0d0,musq,ord)
      Triint(6) = qlI3(0d0,s(i1,i4),0d0,0d0,0d0,0d0,musq,ord)
      Triint(7) = qlI3(s(i4,i5),0d0,0d0,0d0,0d0,0d0,musq,ord)
      Triint(8) = qlI3(s(i3,i6),0d0,0d0,0d0,0d0,0d0,musq,ord)
      Triint(9) = qlI3(0d0,s(i1,i3),0d0,0d0,0d0,0d0,musq,ord)
      Triint(10) = qlI3(s(i3,i5),0d0,0d0,0d0,0d0,0d0,musq,ord)
      Triint(11) = qlI3(s(i3,i4),0d0,0d0,0d0,0d0,0d0,musq,ord)
      Triint(12) = qlI3(s(i2,i6),0d0,0d0,0d0,0d0,0d0,musq,ord)
      Triint(13) = qlI3(s(i2,i5),0d0,0d0,0d0,0d0,0d0,musq,ord)
      Triint(14) = qlI3(s(i2,i4),0d0,0d0,0d0,0d0,0d0,musq,ord)
      Triint(15) = qlI3(s(i2,i3),0d0,0d0,0d0,0d0,0d0,musq,ord)
      Triint(16) = qlI3(0d0,s(i2,i3),t(i4,i5,i6),0d0,0d0,0d0,musq,ord)
      Triint(17) = qlI3(0d0,t(i2,i3,i4),s(i5,i6),0d0,0d0,0d0,musq,ord)
      Triint(18) = qlI3(s(i1,i2),0d0,t(i4,i5,i6),0d0,0d0,0d0,musq,ord)
      Triint(19) = qlI3(s(i1,i2),t(i3,i4,i5),0d0,0d0,0d0,0d0,musq,ord)
      Triint(20) = qlI3(t(i1,i2,i3),0d0,s(i5,i6),0d0,0d0,0d0,musq,ord)
      Triint(21) = qlI3(t(i1,i2,i3),s(i4,i5),0d0,0d0,0d0,0d0,musq,ord)
      Triint(22) = qlI3(s(i1,i2),t(i3,i4,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(23) = qlI3(t(i1,i2,i3),s(i4,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(24) = qlI3(0d0,t(i2,i3,i5),s(i4,i6),0d0,0d0,0d0,musq,ord)
      Triint(25) = qlI3(s(i1,i2),t(i3,i5,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(26) = qlI3(0d0,t(i2,i3,i6),s(i4,i5),0d0,0d0,0d0,musq,ord)
      Triint(27) = qlI3(0d0,s(i2,i4),t(i3,i5,i6),0d0,0d0,0d0,musq,ord)
      Triint(28) = qlI3(t(i1,i2,i4),0d0,s(i5,i6),0d0,0d0,0d0,musq,ord)
      Triint(29) = qlI3(t(i1,i2,i4),s(i3,i5),0d0,0d0,0d0,0d0,musq,ord)
      Triint(30) = qlI3(t(i1,i2,i4),s(i3,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(31) = qlI3(0d0,t(i2,i4,i5),s(i3,i6),0d0,0d0,0d0,musq,ord)
      Triint(32) = qlI3(0d0,t(i2,i4,i6),s(i3,i5),0d0,0d0,0d0,musq,ord)
      Triint(33) = qlI3(0d0,s(i2,i5),t(i3,i4,i6),0d0,0d0,0d0,musq,ord)
      Triint(34) = qlI3(t(i1,i2,i5),0d0,s(i4,i6),0d0,0d0,0d0,musq,ord)
      Triint(35) = qlI3(t(i1,i2,i5),s(i3,i4),0d0,0d0,0d0,0d0,musq,ord)
      Triint(36) = qlI3(t(i1,i2,i5),s(i3,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(37) = qlI3(0d0,t(i2,i5,i6),s(i3,i4),0d0,0d0,0d0,musq,ord)
      Triint(38) = qlI3(0d0,s(i2,i6),t(i3,i4,i5),0d0,0d0,0d0,musq,ord)
      Triint(39) = qlI3(t(i1,i2,i6),0d0,s(i4,i5),0d0,0d0,0d0,musq,ord)
      Triint(40) = qlI3(t(i1,i2,i6),s(i3,i4),0d0,0d0,0d0,0d0,musq,ord)
      Triint(41) = qlI3(t(i1,i2,i6),s(i3,i5),0d0,0d0,0d0,0d0,musq,ord)
      Triint(42) = qlI3(s(i1,i3),0d0,t(i4,i5,i6),0d0,0d0,0d0,musq,ord)
      Triint(43) = qlI3(s(i1,i3),t(i2,i4,i5),0d0,0d0,0d0,0d0,musq,ord)
      Triint(44) = qlI3(s(i1,i3),t(i2,i4,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(45) = qlI3(s(i1,i3),t(i2,i5,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(46) = qlI3(t(i1,i3,i4),0d0,s(i5,i6),0d0,0d0,0d0,musq,ord)
      Triint(47) = qlI3(t(i1,i3,i4),s(i2,i5),0d0,0d0,0d0,0d0,musq,ord)
      Triint(48) = qlI3(t(i1,i3,i4),s(i2,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(49) = qlI3(t(i1,i3,i5),0d0,s(i4,i6),0d0,0d0,0d0,musq,ord)
      Triint(50) = qlI3(t(i1,i3,i5),s(i2,i4),0d0,0d0,0d0,0d0,musq,ord)
      Triint(51) = qlI3(t(i1,i3,i5),s(i2,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(52) = qlI3(t(i1,i3,i6),0d0,s(i4,i5),0d0,0d0,0d0,musq,ord)
      Triint(53) = qlI3(t(i1,i3,i6),s(i2,i4),0d0,0d0,0d0,0d0,musq,ord)
      Triint(54) = qlI3(t(i1,i3,i6),s(i2,i5),0d0,0d0,0d0,0d0,musq,ord)
      Triint(55) = qlI3(s(i1,i4),0d0,t(i3,i5,i6),0d0,0d0,0d0,musq,ord)
      Triint(56) = qlI3(s(i1,i4),t(i2,i3,i5),0d0,0d0,0d0,0d0,musq,ord)
      Triint(57) = qlI3(s(i1,i4),t(i2,i3,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(58) = qlI3(s(i1,i4),t(i2,i5,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(59) = qlI3(t(i1,i4,i5),0d0,s(i3,i6),0d0,0d0,0d0,musq,ord)
      Triint(60) = qlI3(t(i1,i4,i5),s(i2,i3),0d0,0d0,0d0,0d0,musq,ord)
      Triint(61) = qlI3(t(i1,i4,i5),s(i2,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(62) = qlI3(t(i1,i4,i6),0d0,s(i3,i5),0d0,0d0,0d0,musq,ord)
      Triint(63) = qlI3(t(i1,i4,i6),s(i2,i3),0d0,0d0,0d0,0d0,musq,ord)
      Triint(64) = qlI3(t(i1,i4,i6),s(i2,i5),0d0,0d0,0d0,0d0,musq,ord)
      Triint(65) = qlI3(s(i1,i5),0d0,t(i3,i4,i6),0d0,0d0,0d0,musq,ord)
      Triint(66) = qlI3(s(i1,i5),t(i2,i3,i4),0d0,0d0,0d0,0d0,musq,ord)
      Triint(67) = qlI3(s(i1,i5),t(i2,i3,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(68) = qlI3(s(i1,i5),t(i2,i4,i6),0d0,0d0,0d0,0d0,musq,ord)
      Triint(69) = qlI3(t(i1,i5,i6),0d0,s(i3,i4),0d0,0d0,0d0,musq,ord)
      Triint(70) = qlI3(t(i1,i5,i6),s(i2,i3),0d0,0d0,0d0,0d0,musq,ord)
      Triint(71) = qlI3(t(i1,i5,i6),s(i2,i4),0d0,0d0,0d0,0d0,musq,ord)
      Triint(72) = qlI3(s(i1,i6),0d0,t(i3,i4,i5),0d0,0d0,0d0,musq,ord)
      Triint(73) = qlI3(s(i1,i6),t(i2,i3,i4),0d0,0d0,0d0,0d0,musq,ord)
      Triint(74) = qlI3(s(i1,i6),t(i2,i3,i5),0d0,0d0,0d0,0d0,musq,ord)
      Triint(75) = qlI3(s(i1,i6),t(i2,i4,i5),0d0,0d0,0d0,0d0,musq,ord)
      Triint(76) = qlI3(s(i1,i2),s(i3,i4),s(i5,i6),0d0,0d0,0d0,musq,ord)
      Triint(77) = qlI3(s(i1,i2),s(i3,i5),s(i4,i6),0d0,0d0,0d0,musq,ord)
      Triint(78) = qlI3(s(i1,i2),s(i3,i6),s(i4,i5),0d0,0d0,0d0,musq,ord)
      Triint(79) = qlI3(s(i1,i3),s(i2,i4),s(i5,i6),0d0,0d0,0d0,musq,ord)
      Triint(80) = qlI3(s(i1,i3),s(i2,i5),s(i4,i6),0d0,0d0,0d0,musq,ord)
      Triint(81) = qlI3(s(i1,i3),s(i2,i6),s(i4,i5),0d0,0d0,0d0,musq,ord)
      Triint(82) = qlI3(s(i1,i4),s(i2,i3),s(i5,i6),0d0,0d0,0d0,musq,ord)
      Triint(83) = qlI3(s(i1,i4),s(i2,i5),s(i3,i6),0d0,0d0,0d0,musq,ord)
      Triint(84) = qlI3(s(i1,i4),s(i2,i6),s(i3,i5),0d0,0d0,0d0,musq,ord)
      Triint(85) = qlI3(s(i1,i5),s(i2,i3),s(i4,i6),0d0,0d0,0d0,musq,ord)
      Triint(86) = qlI3(s(i1,i5),s(i2,i4),s(i3,i6),0d0,0d0,0d0,musq,ord)
      Triint(87) = qlI3(s(i1,i5),s(i2,i6),s(i3,i4),0d0,0d0,0d0,musq,ord)
      Triint(88) = qlI3(s(i1,i6),s(i2,i3),s(i4,i5),0d0,0d0,0d0,musq,ord)
      Triint(89) = qlI3(s(i1,i6),s(i2,i4),s(i3,i5),0d0,0d0,0d0,musq,ord)
      Triint(90) = qlI3(s(i1,i6),s(i2,i5),s(i3,i4),0d0,0d0,0d0,musq,ord)

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
      PROGRAM TEST
      USE FMZM
      IMPLICIT NONE

!  Least squares fit for the coefficients in the asymptotic series for the Jth harmonic number.

!  H(J) = 1 + 1/2 + 1/3 + ... + 1/J  defines the Jth harmonic number.

!  Find an approximation to H(J) of the form:

!  ln(J) + c(1) + c(2)/J + ... + c(k)/J**(k-1)

!  Integrating 1/x from 1 to J gives ln(J) as a first approximation, and we generate N data
!  points (x(i),y(i)) where x(i) is J and y(i) is H(J) for various J values.  Then we do a
!  least squares fit of the model function c(1) + c(2)/J + ... + c(k)/J**(k-1) to the data
!  (x(i),y(i)-ln(i)).

!  Since this is a sample problem, we can compare the results of the fit to the "true"
!  asymptotic formula, where c(1) = 0.57721566..., Euler's constant, and for i > 1,
!  c(i) = -B(i-1)/(i-1).  The B values are Bernoulli numbers, and the first few are:
!  B(1) = -1/2, B(2) = 1/6, B(4) = -1/30, B(6) = 1/42, ..., with the others being zero:
!  B(3) = B(5) = B(7) = ... = 0.

!  The first c's in the list of fitted coefficients give the most agreement with the
!  theoretical values, and the last ones the least.  The linear system is ill-conditioned,
!  but by using high precision we can get good accuracy for several coefficients.
!  For example, using 400 digit precision, 60 data points at intervals of 100 (i.e.,
!  x(i) = 100, 200, 300, ..., 6000), and fitting 60 coefficients, we get at least 50
!  decimal agreement between the fitted c's and the theoretical ones for c(1), ..., c(29).
!  c(41) agrees to 16 decimals, and because the number is large this is 31 significant
!  digit agreement.

      INTEGER :: J, K, N, NGAP
      TYPE (FM) :: H_N, ONE, DET
      TYPE (FM), ALLOCATABLE :: A(:,:), B(:), C(:), X(:), Y(:)
      TYPE (FM), EXTERNAL :: F

!             This is not a good way to compute Euler's constant, but with 150 digit precision,
!             N = 40 data points at intervals of NGAP = 10, fitting K = 40 coefficients we get
!             c(1) = .57721566490153286060651209008240243104215933593992,
!             correct to 50 places.

!             Set FM precision.

      CALL FM_SET(150)

!             N is the number of harmonic data points.

      N = 40

!             NGAP is the gap between harmonic data points.

      NGAP = 10

!             K is the number of coefficients to fit.

      K = 40

      ALLOCATE(A(K,K),B(K),C(K),X(N),Y(N),STAT=J)
      IF (J /= 0) THEN
          WRITE (*,"(/' Error in HFIT.  Unable to allocate arrays with K,N = ',2I8/)") K,N
          STOP
      ENDIF

!             Generate the harmonic data points.
!             Since the coefficient of the first term in the model, ln(x), is assumed
!             to be 1 and is not being fitted, subtract that from the Y data points.

      H_N = 0
      ONE = 1
      WRITE (*,*) ' '
      WRITE (*,*) ' Data points:'
      WRITE (*,*) ' '
      DO J = 1, N*NGAP
         H_N = H_N + ONE/J
         IF (MOD(J,NGAP) == 0) THEN
             X(J/NGAP) = J
             Y(J/NGAP) = H_N - LOG(X(J/NGAP))
             WRITE (*,"(A,I4,A,I6,A,A)") ' I = ',J/NGAP,'  X = ',J,'  Y = ',  &
                                         TRIM(FM_FORMAT('F40.35',Y(J/NGAP)))
         ENDIF
      ENDDO

!             Generate the linear system for the normal equations.

      CALL FM_GENEQ(F,A,B,K,X,Y,N)

!             Solve the linear system for the normal equations.

      CALL FM_LIN_SOLVE(A,C,B,N,DET)

!             Print the solution.
!             When using F format, FM doesn't like to print 0.00000...0 showing no
!             significant digits when the actual number is too small for that format.
!             FM will shift to E format when possible, to avoid showing all zeroes.
!             In this example, all the even-numbered coefficients are zero in the
!             asymptotic series for the harmonic numbers, so any non-zero digits
!             found in the fit are not interesting.  Therefore the if statement
!             below prints exactly zero when C(J) is too small, making the output
!             look neater.

      WRITE (*,*) ' '
      WRITE (*,*) ' Fitted coefficients:'
      DO J = 1, K
         IF (ABS(C(J)) > 1.0D-50) THEN
             WRITE (*,"(A,I3,A,A)") ' J = ',J,' C(J) = ',TRIM(FM_FORMAT('F60.50',C(J)))
         ELSE
             WRITE (*,"(A,I3,A,A)") ' J = ',J,' C(J) = ',TRIM(FM_FORMAT('F60.50',TO_FM(0)))
         ENDIF
      ENDDO

      END PROGRAM TEST

      FUNCTION F(J,X)     RESULT (RETURN_VALUE)
      USE FMZM
      IMPLICIT NONE

!  This defines the model function being fitted to the data points.
!  For the harmonic number case, the model function is:

!  F(J,X) = 1/X**(J-1)

!  This will fit the terms  c1 + c2/n + c3/n**2 + ... to the harmonic model function
!  ln(x) + c1 + c2/n + c3/n**2 + ....

      INTEGER :: J
      TYPE (FM) :: RETURN_VALUE, X

      RETURN_VALUE = 1/X**(J-1)

      END FUNCTION F

max_stars_repo_name: andreypudov/projecteuler

max_issues_repo_name: andreypudov/projecteuler

max_forks_repo_name: andreypudov/projecteuler


! Chemfiles, an efficient IO library for chemistry file formats
! Copyright (C) 2015-2019 Guillaume Fraux -- BSD licence
module chemfiles_property
    use iso_c_binding
    use chemfiles_ffi
    use chemfiles_strings
    implicit none

    private
    public :: chfl_property

    type, extends(chfl_ptr) :: chfl_property
    contains
        generic :: init => bool_init, double_init, string_init, vector3d_init
        procedure, private :: bool_init
        procedure, private :: double_init
        procedure, private :: string_init
        procedure, private :: vector3d_init

        procedure :: kind => get_kind
        procedure :: bool
        procedure :: double
        procedure :: string
        procedure :: vector3d
    end type

contains
    subroutine bool_init(this, value, status)
        implicit none
        class(chfl_property), intent(out) :: this
        logical, intent(in) :: value
        integer(chfl_status), optional, intent(out) :: status

        call this%unsafe_set_ptr(c_chfl_property_bool(logical(value, c_bool)), status)
    end subroutine

    subroutine double_init(this, value, status)
        implicit none
        class(chfl_property), intent(inout) :: this
        real(kind=c_double), intent(in) :: value
        integer(chfl_status), intent(out), optional :: status

        call this%unsafe_set_ptr(c_chfl_property_double(value), status)
    end subroutine

    subroutine string_init(this, value, status)
        implicit none
        class(chfl_property), intent(inout) :: this
        character(len=*), intent(in) :: value
        integer(chfl_status), intent(out), optional :: status

        call this%unsafe_set_ptr(c_chfl_property_string(f_to_c_str(value)), status)
    end subroutine

    subroutine vector3d_init(this, value, status)
        implicit none
        class(chfl_property), intent(inout) :: this
        real(kind=c_double), dimension(3), intent(in) :: value
        integer(chfl_status), intent(out), optional :: status

        call this%unsafe_set_ptr(c_chfl_property_vector3d(value), status)
    end subroutine

    function get_kind(this, status)
        implicit none
        class(chfl_property), intent(in) :: this
        integer(chfl_status), intent(out), optional :: status
        integer(chfl_property_kind) :: get_kind
        integer(chfl_status) :: status_tmp

        status_tmp = c_chfl_property_get_kind(this%unsafe_const_ptr(), get_kind)

        if (present(status)) then
            status = status_tmp
        end if
    end function

    logical function bool(this, status)
        implicit none
        class(chfl_property), intent(in) :: this
        integer(chfl_status), intent(out), optional :: status

        integer(chfl_status) :: status_tmp
        logical(c_bool) :: value_tmp

        status_tmp = c_chfl_property_get_bool(this%unsafe_const_ptr(), value_tmp)
        bool = logical(value_tmp)

        if (present(status)) then
            status = status_tmp
        end if
    end function

    real(kind=c_double) function double(this, status)
        implicit none
        class(chfl_property), intent(in) :: this
        integer(chfl_status), intent(out), optional :: status
        integer(chfl_status) :: status_tmp

        status_tmp = c_chfl_property_get_double(this%unsafe_const_ptr(), double)

        if (present(status)) then
            status = status_tmp
        end if
    end function

    character(len=CHFL_STRING_LENGTH) function string(this, status)
        implicit none
        class(chfl_property), intent(in) :: this
        integer(chfl_status), intent(out), optional :: status
        integer(chfl_status) :: status_tmp

        status_tmp = c_chfl_property_get_string(this%unsafe_const_ptr(), string, int(CHFL_STRING_LENGTH, c_int64_t))
        string = rm_null_in_str(string)
        if (present(status)) then
            status = status_tmp
        end if
    end function

    function vector3d(this, status)
        implicit none
        class(chfl_property), intent(in) :: this
        integer(chfl_status), intent(out), optional :: status
        real(kind=c_double), dimension(3) :: vector3d
        integer(chfl_status) :: status_tmp

        status_tmp = c_chfl_property_get_vector3d(this%unsafe_const_ptr(), vector3d)

        if (present(status)) then
            status = status_tmp
        end if
    end function
end module


text:
\ surface
\
\ Copyright (C) 2016 David J Goehrig <dave@dloh.org>
\
\  This software is provided 'as-is', without any express or implied
\  warranty.  In no event will the authors be held liable for any damages
\  arising from the use of this software.
\
\  Permission is granted to anyone to use this software for any purpose,
\  including commercial applications, and to alter it and redistribute it
\  freely, subject to the following restrictions:
\
\  1. The origin of this software must not be misrepresented; you must not
\     claim that you wrote the original software. If you use this software
\     in a product, an acknowledgment in the product documentation would be
\     appreciated but is not required.
\  2. Altered source versions must be plainly marked as such, and must not be
\     misrepresented as being the original software.
\  3. This notice may not be removed or altered from any source distribution.
\

$00000000 constant SDL_SWSURFACE
$00000001 constant SDL_PREALLOC
$00000002 constant SDL_RLEACCEL
$00000004 constant SDL_DONTFREE

: SDL_Surface 0 value ;

: SDL_Surface:flags ;
: SDL_Surface:format 1 cells + ;
: SDL_Surface:w 2 cells + ;
: SDL_Surface:h 3 cells + ;
: SDL_Surface:pitch 4 cells + ;
: SDL_Surface:pixels 5 cells + ;
: SDL_Surface:userdata 6 cells + ;
: SDL_Surface:locked 7 cells + ;
: SDL_Surface:lock_data 8 cells + ;
: SDL_Surface:clip_rect 9 cells + ;
: SDL_Surface:map ." SDL_Surface:map not implemented " ;
: SDL_Surface:refcount ." SDL_Surface:refcount not implemented " ;

: SDL_MUSTLOCK @ SDL_RLEACCEL and 0= ;

FUNCTION: SDL_CreateRGBSurface ( flags w h d r g b a -- surface )
FUNCTION: SDL_CreateRGBSurfaceFrom ( pixels w h d p r g b a -- surface )
FUNCTION: SDL_FreeSurface ( surface -- )
FUNCTION: SDL_SetSurfacePalette ( surface palette -- flag )
FUNCTION: SDL_LockSurface ( surface -- flag )
FUNCTION: SDL_UnlockSurface ( surface -- )
FUNCTION: SDL_LoadBMP_RW ( src n -- surface )

\ need to include rwops first!
: SDL_LoadBMP ( z -- surface ) z" rb" SDL_RWFromFile 1 SDL_LoadBMP_RW ;

FUNCTION: SDL_SaveBMP_RW ( surface dst n -- flag )

: SDL_SaveBMP ( surface z -- flag ) z" wb" SDL_RWFromFile 1 SDL_SaveBMP_RW ;

FUNCTION: SDL_SetSurfaceRLE ( surface flag -- flag )
FUNCTION: SDL_SetColorKey ( surface flag key -- flag )
FUNCTION: SDL_GetColorKey ( surface key* -- flag )
FUNCTION: SDL_SetSurfaceColorMod ( surface r g b -- flag )
FUNCTION: SDL_GetSurfaceColorMod ( surface r* g* b* -- flag )
FUNCTION: SDL_SetSurfaceAlphaMod ( surface a -- flag )
FUNCTION: SDL_GetSurfaceAlphaMod ( surface a* -- flag )
FUNCTION: SDL_SetSurfaceBlendMode ( surface mode -- flag )
FUNCTION: SDL_GetSurfaceBlendMode ( surface mode* -- flag )
FUNCTION: SDL_SetClipRect ( surface rect -- flag )
FUNCTION: SDL_GetClipRect ( surface rect* -- )
FUNCTION: SDL_ConvertSurface ( surface fmt flags -- surface )
FUNCTION: SDL_ConvertSurfaceFormat ( surface n flags -- surface )
FUNCTION: SDL_ConvertPixels ( w h fmt src* pitch fmt dst* pitch -- flag )
FUNCTION: SDL_FillRect ( surface rect color -- flag )
FUNCTION: SDL_FillRects ( surface rect* count color -- flag )
FUNCTION: SDL_UpperBlit ( src rect dst rect -- flag )

: SDL_BlitSurface SDL_UpperBlit ;

FUNCTION: SDL_LowerBlit ( src rect dst rect -- flag )
FUNCTION: SDL_SoftStretch ( src rect dst rect -- flag )
FUNCTION: SDL_UpperBlitScaled ( src rect dst rect -- flag )

: SDL_BlitScaled SDL_UpperBlitScaled ;

FUNCTION: SDL_LowerBlitScaled ( src rect dst rect -- flag )

max_stars_repo_name: cthulhuology/sdlforth

max_issues_repo_name: cthulhuology/sdlforth

max_forks_repo_name: cthulhuology/sdlforth


!this subroutine substitutes the extract and the split routine in the specific case of Erin's light curves

subroutine extract_split_erin(filename, o)
  use dyn_lc
  implicit none

  character (len = 500)  :: filename
  integer, INTENT(IN)    :: o
  integer :: j,  read_file_length
  double precision :: temp, temp_time

  int_len_dim = read_file_length(filename) - 1
  int_number = 1

  ! write(*,*) 'length of the light curve', int_len_dim
  open(2, file = filename)

  read(2,*) !skip the first line which is the header
  do j = 1, int_len_dim
     ! read(2, *) temp, time_obs(1, j, obs_index),
     read(2,*)  temp_time, lc_en_obs(1,j,1,o), temp, lc_en_obs(1,j,2,o), temp, lc_en_obs(1,j,3,o), temp, lc_en_obs(1,j,4,o), temp, lc_en_obs(1,j,5,o), temp, lc_en_obs(1,j,6,o), temp, lc_en_obs(1,j,7,o), temp, lc_en_obs(1,j,8,o), temp, lc_en_obs(1,j,9,o), temp, lc_en_obs(1,j,10,o), temp, lc_en_obs(1,j,11,o), temp, lc_en_obs(1,j,12,o), temp, lc_en_obs(1,j,13,o), temp, lc_en_obs(1,j,14,o), temp, lc_en_obs(1,j,15,o), temp, lc_en_obs(1,j,16,o), temp, lc_en_obs(1,j,17,o), temp, lc_en_obs(1,j,18,o), temp, lc_en_obs(1,j,19,o), temp, lc_en_obs(1,j,20,o), temp, lc_en_obs(1,j,21,o), temp, lc_en_obs(1,j,22,o), temp, lc_en_obs(1,j,23,o), temp, lc_en_obs(1,j,24,o), temp, lc_en_obs(1,j,25,o), temp, lc_en_obs(1,j,26,o), temp, lc_en_obs(1,j,27,o), temp, lc_en_obs(1,j,28,o), temp, lc_en_obs(1,j,29,o), temp, lc_en_obs(1,j,30,o), temp, lc_en_obs(1,j,31,o), temp, lc_en_obs(1,j,32,o), temp, lc_en_obs(1,j,33,o), temp, lc_en_obs(1,j,34,o), temp, lc_en_obs(1,j,35,o), temp
     ! write(*,*) j, temp_time

  enddo
  close(2)

 end subroutine extract_split_erin