text:
!***************************************************************************************************************************
! This module is intended to calculate environmental factors
! 15SEP2017 added function TFAC5 to estimate cardinal temperatures with a plateau by @lpmorenoc
! Atributes:
!   - tMin minimmun temprature of the day
!   - tMax maximun temprature of the day
!   - dewPointTemp dew point temperature of the day
!   - dayRadiation registered during the day
! Type functions:
!        hourlyTemperature
!        hourlySVP
!        hourlyWHCAIR
!        hourlyRH
!        hourlyVPD
!        hourlyRadiation
!        hourlyBiomass
!        hourlyTranspiration
! Static functions:
!        calculateSVP
!        calculateWHCAIR
!        calculateRH
!        calculateVPD
!        calculatePortionOfRadiation
!        calculatateHourlyRadiation
!        calculateBiomass
!        calculateTranspiration
!        calculateTemperatureFactor
! Authors
! @danipilze
!*********

    Module YCA_Control_Environment         !Module of environment

    ! STATIC ATRIBUTES
    real, private :: PI =  4 * atan (1.0_8)
    real, private :: HOURS_OF_DAY = 24
    real, private :: LIGHT_HOURS = 12
    real, private :: DAWN_TIME = 5                                              !dawn time

    type DailyEnvironment_type

        ! OBJECT ATRIBUTES

        real, private :: tMin_ = 0
        real, private :: tMax_ = 0
        real, private :: dewPointTemp_ = 0
        real, private :: dayRadiation_ = 0                                      ! solar radiation

    contains

        ! type functions
        procedure, pass (this) :: hourlyTemperature
        procedure, pass (this) :: hourlySVP
        procedure, pass (this) :: hourlyWHCAIR
        procedure, pass (this) :: hourlyRH
        procedure, pass (this) :: hourlyVPD
        procedure, pass (this) :: hourlyRadiation
        procedure, pass (this) :: hourlyBiomass
        procedure, pass (this) :: hourlyTranspiration

        ! getters and setters
        procedure, pass (this) :: getTMin
        procedure, pass (this) :: setTMin
        procedure, pass (this) :: getTMax
        procedure, pass (this) :: setTMax
        procedure, pass (this) :: getDewPointTemp
        procedure, pass (this) :: setDewPointTemp
        procedure, pass (this) :: getDayRadiation
        procedure, pass (this) :: setDayRadiation

    end Type DailyEnvironment_type

    ! interface to reference the constructor
    interface DailyEnvironment_type
        module procedure DailyEnvironment_type_constructor
    end interface DailyEnvironment_type

    contains

    ! constructor for the type
    ! tMin          mandatory
    ! tMax          mandatory
    ! dewPointTemp
    ! dayRadiation
    type (DailyEnvironment_type) function DailyEnvironment_type_constructor(tMin, tMax, dewPointTemp, dayRadiation)
        implicit none
        real, intent (in) :: tMin, tMax, dewPointTemp,dayRadiation
        DailyEnvironment_type_constructor%TMin_ = tMin
        DailyEnvironment_type_constructor%TMax_ = tMax
        DailyEnvironment_type_constructor%dewPointTemp_ = dewPointTemp
        DailyEnvironment_type_constructor%dayRadiation_ = dayRadiation
    end function DailyEnvironment_type_constructor

    !-------------------------------------------
    ! TYPE FUNCTIONS
    !-------------------------------------------

    ! obtain the temperature for any given hour of the day
    ! T(Hour) = Amplitude*sin[w(t - a)] + C.
    ! Amplitude is called the amplitude the height of each peak above the baseline
    ! Hod is Hours Of Day, the period or wavelength (the length of each cycle)
    ! a  is the phase shift (the horizontal offset of the basepoint; where the curve crosses the baseline as it ascends)
    ! C is average temperature,  the vertical offset (height of the baseline)
    ! w is the angular frequency, given by w = 2PI/hod
    real function hourlyTemperature(this, hour)
        implicit none
        class (DailyEnvironment_type), intent(in) :: this
        integer, intent (in) :: hour
        real :: Amplitude, C, w, g

        w = (2*PI)/HOURS_OF_DAY
        Amplitude = ((this%TMax_ - this%TMin_)/2)                                               ! half distance between temperatures
        C = (this%TMin_ + this%TMax_)/2                                                         ! mean temperature
        g = w*(hour - 8)

        hourlyTemperature = Amplitude*SIN(g)+C                                                  ! calculate temperature acording to the current time

    end function hourlyTemperature

    !obtain the Saturation Vapour Pressure for any given hour                                   ! (pascals)
    real function hourlySVP(this, Hour)
        implicit none
        class (DailyEnvironment_type), intent(in) :: this
        integer, intent (in) :: hour

        hourlySVP = calculateSVP(hourlyTemperature(this, Hour))

    end function hourlySVP

    ! obtain the water holding capacity of the air for any given hour
    real function hourlyWHCAIR(this, hour) ! (kg/m3)
        implicit none
        class (DailyEnvironment_type), intent(in) :: this
        integer, intent (in) :: hour

        hourlyWHCAIR = calculateWHCAIR(hourlyTemperature(this, Hour))

    end function hourlyWHCAIR

    ! obtain the relative humidity for any given hour
    real function hourlyRH(this, hour)
        implicit none
        class (DailyEnvironment_type), intent(in) :: this
        integer, intent (in) :: hour

        hourlyRH =  calculateRH(hourlyTemperature(this, Hour), this%dewPointTemp_)

    end function hourlyRH

    ! obtain the Vapor Preasure Deficit VPD for any given hour
    real function hourlyVPD(this, hour)
        implicit none
        class (DailyEnvironment_type), intent(in) :: this
        integer, intent (in) :: hour

        hourlyVPD =  calculateVPD(hourlyTemperature(this, Hour), this%dewPointTemp_)

    end function hourlyVPD

    ! obtain the radiation VPD for any given hour
    real function hourlyRadiation(this, hour)
        implicit none
        class (DailyEnvironment_type), intent(in) :: this
        integer, intent (in) :: hour

        hourlyRadiation = calculatateHourlyRadiation(hour, DAWN_TIME, this%dayRadiation_ , LIGHT_HOURS)

    end function hourlyRadiation

    ! obtain the biomass produced at a given hour
    real function hourlyBiomass(this, hour, extinctionCoefficient,LAI, radiationUseEfficiency,stomatalConductance)
        implicit none
        class (DailyEnvironment_type), intent(in) :: this
        integer, intent (in) :: hour
        real, intent (in) :: extinctionCoefficient, LAI,radiationUseEfficiency,stomatalConductance

        hourlyBiomass = calculateBiomass(calculatateHourlyRadiation(hour, DAWN_TIME, this%dayRadiation_ , LIGHT_HOURS) , extinctionCoefficient , LAI, radiationUseEfficiency, stomatalConductance )

    end function hourlyBiomass

    ! obtain the transpiration
    real function hourlyTranspiration(this, hour, stomatalConductance)
        implicit none
        class (DailyEnvironment_type), intent(in) :: this
        integer, intent (in) :: hour
        real, intent (in) :: stomatalConductance

        hourlyTranspiration = calculateTranspiration(hourlyVPD(this, hour), stomatalConductance)

    end function hourlyTranspiration

    !-------------------------------------------
    ! STATIC FUNCTIONS
    !-------------------------------------------

    ! obtain the Saturation Vapour Pressure  for any given temperature                          ! (pascals)
    real function calculateSVP(temperature)
        implicit none
        real, intent (in) :: temperature

        !Saturation vapour pressure in pascals: svp = 610.78 *exp( t / ( t + 238.3 ) *17.2694 )
        calculateSVP = 610.78 * exp( temperature/ ( temperature  + 238.3 ) * 17.2694 )

    end function calculateSVP

    ! obtain the water holding capacity of the air for any given temperature                    ! (kg/m3)
    real function calculateWHCAIR(temperature)
        implicit none
        real, intent (in) :: temperature

        ! WHCAIR = 0.002166 *SVP / ( t + 273.16 )
        calculateWHCAIR = 0.002166 * calculateSVP(temperature) / ( temperature  + 273.16 )

    end function calculateWHCAIR

    ! obtain the relative humidity                                                              ! %
    real function calculateRH(temperature, dewPointTemp)
        implicit none
        real, intent (in) :: temperature, dewPointTemp

        !(RHt) = WHCAIRdp/WHCAIRt
        calculateRH =  calculateWHCAIR(dewPointTemp)/calculateWHCAIR(temperature)

    end function calculateRH

    ! obtain the Vapor Preasure Deficit VPD                                                     ! (KPa)
    real function calculateVPD(temperature, dewPointTemp)
        implicit none
        real, intent (in) :: temperature, dewPointTemp

        ! VPD = (1 - (RH/100)) * SVP
        calculateVPD =  ((1 - (calculateRH(temperature, dewPointTemp)/100)) * calculateSVP(temperature))/1000

    end function calculateVPD

    ! calculates proportion of radiation
    ! K is the extinction coefficient
    ! LAI is leaf area index
    real function calculatePortionOfRadiation(K, LAI)                                           ! fraction
        implicit none
        real, intent (in) :: K, LAI
        real :: value = 0

        calculatePortionOfRadiation = 1.0 - exp(-K*LAI)

    end function calculatePortionOfRadiation

    ! obtain the radiation at a given hour
    real function calculatateHourlyRadiation(hour, dawnTime, dayRadiation, lightHours)          !MJ/m2
        implicit none
        integer, intent (in) :: hour
        real :: Amplitude, C, w, g, value, dawnTime, dayRadiation, lightHours

        w = (2*PI)/24
        C = 0
        Amplitude = (dayRadiation * PI)/24
        g = w*(hour - dawnTime)

        value = Amplitude*SIN(g)+C           ! calculate hourly radiation
        if(value < 0.0 ) value = 0           ! no radiation on non-light hours

        calculatateHourlyRadiation=value

    end function calculatateHourlyRadiation

    ! obtain the biomass produced accoording to the radiation at a given moment             ! MJ/m2 * # * # * g/MJ * fraction = g/m2
    real function calculateBiomass(radiation, extinctionCoefficient , LAI, radiationUseEfficiency, stomatalConductance)
        implicit none
        real, intent (in) :: radiation, extinctionCoefficient, LAI, radiationUseEfficiency, stomatalConductance
        real :: value = 0

        calculateBiomass = calculatePortionOfRadiation(extinctionCoefficient, LAI) * radiation * radiationUseEfficiency  * stomatalConductance

    end function calculateBiomass

    ! obtain the transpiration
    real function calculateTranspiration(VPD, stomatalConductance)
        implicit none
        real, intent (in) :: stomatalConductance, VPD

        calculateTranspiration = stomatalConductance * VPD

    end function calculateTranspiration

    ! Estimate thermal time with a plateau after second cardinal temperature                        ! by @lpmorenoc
    ! Calculate temp factor and thermal units from cardinal temps                                   ! 0EQN 58
    ! ARGUMENTS
    ! tcard,temp,tunit
    ! MODIFES
    ! tunit
    ! RETURNS
    ! tfac5
    real function calculateTemperatureFactor(tcard,temp,tunit)

      implicit none
      real, intent (in) ::      tcard(4),temp
      real                      tfac5
      real                      tunit

      IF (temp <= tcard(1)) THEN
        tfac5 = 0.0
      ELSEIF (temp > tcard(1) .AND. temp <= tcard(2)) THEN
        tfac5 = (temp-tcard(1))/(tcard(2)-tcard(1))
      ELSEIF (temp > tcard(2) .AND. temp <= tcard(3)) THEN
        tfac5 = 1.0
      ELSEIF (temp > tcard(3) .AND. temp <= tcard(4) .AND. tcard(3) < tcard(4)) THEN
        tfac5 = 1.0 - ((temp-tcard(3))/(tcard(4)-tcard(3)))
      ELSEIF (temp > tcard(4) .AND. tcard(3) < tcard(4)) THEN
        tfac5 = 0.0
      ELSEIF (temp > tcard(3)) THEN
        tfac5 = 1.0
      ENDIF

      calculateTemperatureFactor = AMAX1(0.0,AMIN1(1.0,tfac5))
      tunit = tfac5 * (tcard(2)-tcard(1))

      end function calculateTemperatureFactor
    !-------------------------------------------
    ! GETTERS AND SETTERS
    !------------------------------------------

    ! get TMin
    real function getTMin(this)
        implicit none
        class (DailyEnvironment_type), intent(in) :: this

        getTMin = this%tMin_
    end function getTMin

    ! set TMin
    subroutine setTMin(this, tMin)
        implicit none
        class (DailyEnvironment_type), intent(inout) :: this
        real, intent (in) :: tMin

        this%tMin_ = tMin
    end subroutine setTMin

    ! get tMax
    real function getTMax(this)
        implicit none
        class (DailyEnvironment_type), intent(in) :: this

        getTMax = this%tMax_
    end function getTMax

    ! set tMax
    subroutine setTMax(this, tMax)
        implicit none
        class (DailyEnvironment_type), intent(inout) :: this
        real, intent (in) :: tMax

        this%tMax_ = tMax
    end subroutine setTMax

    ! get dewPointTemp
    real function getDewPointTemp(this)
        implicit none
        class (DailyEnvironment_type), intent(in) :: this

        getDewPointTemp = this%dewPointTemp_
    end function getDewPointTemp

    ! set dewPointTemp
    subroutine setDewPointTemp(this, dewPointTemp)
        implicit none
        class (DailyEnvironment_type), intent(inout) :: this
        real, intent (in) :: dewPointTemp

        this%dewPointTemp_ = dewPointTemp
    end subroutine setDewPointTemp

    ! get day radiation
    real function getDayRadiation(this)
        implicit none
        class (DailyEnvironment_type), intent(in) :: this

        getDayRadiation = this%dayRadiation_
    end function getDayRadiation

    ! set day radiation
    subroutine setDayRadiation(this, dayRadiation)
        implicit none
        class (DailyEnvironment_type), intent(inout) :: this
        real, intent (in) :: dayRadiation

        this%dayRadiation_ = dayRadiation
    end subroutine setDayRadiation

END Module YCA_Control_Environment

max_stars_repo_name: vstocca/dssat-csm-os

max_issues_repo_name: vstocca/dssat-csm-os

max_forks_repo_name: vstocca/dssat-csm-os


text:
MODULE BARDELL

CONTAINS

SUBROUTINE integral_ff(i, j, x1t, x1r, x2t, x2r, y1t, y1r, y2t, y2r, out)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: i, j
    REAL*8, INTENT(IN) :: x1t, x1r, x2t, x2r, y1t, y1r, y2t, y2r
    REAL*8, INTENT(OUT) :: out
    out = 0
    SELECT CASE (i)
    CASE (1)
        SELECT CASE (j)
        CASE (1)
            out = 0.742857142857143*x1t*y1t
            RETURN
        CASE (2)
            out = 0.104761904761905*x1t*y1r
            RETURN
        CASE (3)
            out = 0.257142857142857*x1t*y2t
            RETURN
        CASE (4)
            out = -0.0619047619047619*x1t*y2r
            RETURN
        CASE (5)
            out = 0.0666666666666667*x1t
            RETURN
        CASE (6)
            out = -0.0126984126984127*x1t
            RETURN
        CASE (8)
            out = 0.000288600288600289*x1t
            RETURN
        END SELECT
    CASE (2)
        SELECT CASE (j)
        CASE (1)
            out = 0.104761904761905*x1r*y1t
            RETURN
        CASE (2)
            out = 0.019047619047619*x1r*y1r
            RETURN
        CASE (3)
            out = 0.0619047619047619*x1r*y2t
            RETURN
        CASE (4)
            out = -0.0142857142857143*x1r*y2r
            RETURN
        CASE (5)
            out = 0.0142857142857143*x1r
            RETURN
        CASE (6)
            out = -0.00158730158730159*x1r
            RETURN
        CASE (7)
            out = -0.000529100529100529*x1r
            RETURN
        CASE (8)
            out = 0.000144300144300144*x1r
            RETURN
        END SELECT
    CASE (3)
        SELECT CASE (j)
        CASE (1)
            out = 0.257142857142857*x2t*y1t
            RETURN
        CASE (2)
            out = 0.0619047619047619*x2t*y1r
            RETURN
        CASE (3)
            out = 0.742857142857143*x2t*y2t
            RETURN
        CASE (4)
            out = -0.104761904761905*x2t*y2r
            RETURN
        CASE (5)
            out = 0.0666666666666667*x2t
            RETURN
        CASE (6)
            out = 0.0126984126984127*x2t
            RETURN
        CASE (8)
            out = -0.000288600288600289*x2t
            RETURN
        END SELECT
    CASE (4)
        SELECT CASE (j)
        CASE (1)
            out = -0.0619047619047619*x2r*y1t
            RETURN
        CASE (2)
            out = -0.0142857142857143*x2r*y1r
            RETURN
        CASE (3)
            out = -0.104761904761905*x2r*y2t
            RETURN
        CASE (4)
            out = 0.019047619047619*x2r*y2r
            RETURN
        CASE (5)
            out = -0.0142857142857143*x2r
            RETURN
        CASE (6)
            out = -0.00158730158730159*x2r
            RETURN
        CASE (7)
            out = 0.000529100529100529*x2r
            RETURN
        CASE (8)
            out = 0.000144300144300144*x2r
            RETURN
        END SELECT
    CASE (5)
        SELECT CASE (j)
        CASE (1)
            out = 0.0666666666666667*y1t
            RETURN
        CASE (2)
            out = 0.0142857142857143*y1r
            RETURN
        CASE (3)
            out = 0.0666666666666667*y2t
            RETURN
        CASE (4)
            out = -0.0142857142857143*y2r
            RETURN
        CASE (5)
            out = 0.0126984126984127
            RETURN
        CASE (7)
            out = -0.000769600769600770
            RETURN
        CASE (9)
            out = 4.44000444000444e-5
            RETURN
        END SELECT
    CASE (6)
        SELECT CASE (j)
        CASE (1)
            out = -0.0126984126984127*y1t
            RETURN
        CASE (2)
            out = -0.00158730158730159*y1r
            RETURN
        CASE (3)
            out = 0.0126984126984127*y2t
            RETURN
        CASE (4)
            out = -0.00158730158730159*y2r
            RETURN
        CASE (6)
            out = 0.00115440115440115
            RETURN
        CASE (8)
            out = -0.000177600177600178
            RETURN
        CASE (10)
            out = 1.48000148000148e-5
            RETURN
        END SELECT
    CASE (7)
        SELECT CASE (j)
        CASE (2)
            out = -0.000529100529100529*y1r
            RETURN
        CASE (4)
            out = 0.000529100529100529*y2r
            RETURN
        CASE (5)
            out = -0.000769600769600770
            RETURN
        CASE (7)
            out = 0.000266400266400266
            RETURN
        CASE (9)
            out = -5.92000592000592e-5
            RETURN
        CASE (11)
            out = 6.09412374118256e-6
            RETURN
        END SELECT
    CASE (8)
        SELECT CASE (j)
        CASE (1)
            out = 0.000288600288600289*y1t
            RETURN
        CASE (2)
            out = 0.000144300144300144*y1r
            RETURN
        CASE (3)
            out = -0.000288600288600289*y2t
            RETURN
        CASE (4)
            out = 0.000144300144300144*y2r
            RETURN
        CASE (6)
            out = -0.000177600177600178
            RETURN
        CASE (8)
            out = 8.88000888000888e-5
            RETURN
        CASE (10)
            out = -2.43764949647303e-5
            RETURN
        CASE (12)
            out = 2.88669019319174e-6
            RETURN
        END SELECT
    CASE (9)
        SELECT CASE (j)
        CASE (5)
            out = 4.44000444000444e-5
            RETURN
        CASE (7)
            out = -5.92000592000592e-5
            RETURN
        CASE (9)
            out = 3.65647424470954e-5
            RETURN
        CASE (11)
            out = -1.15467607727670e-5
            RETURN
        CASE (13)
            out = 1.51207581548139e-6
            RETURN
        END SELECT
    CASE (10)
        SELECT CASE (j)
        CASE (6)
            out = 1.48000148000148e-5
            RETURN
        CASE (8)
            out = -2.43764949647303e-5
            RETURN
        CASE (10)
            out = 1.73201411591504e-5
            RETURN
        CASE (12)
            out = -6.04830326192555e-6
            RETURN
        CASE (14)
            out = 8.54651547880785e-7
            RETURN
        END SELECT
    CASE (11)
        SELECT CASE (j)
        CASE (7)
            out = 6.09412374118256e-6
            RETURN
        CASE (9)
            out = -1.15467607727670e-5
            RETURN
        CASE (11)
            out = 9.07245489288833e-6
            RETURN
        CASE (13)
            out = -3.41860619152314e-6
            RETURN
        CASE (15)
            out = 5.12790928728471e-7
            RETURN
        END SELECT
    CASE (12)
        SELECT CASE (j)
        CASE (8)
            out = 2.88669019319174e-6
            RETURN
        CASE (10)
            out = -6.04830326192555e-6
            RETURN
        CASE (12)
            out = 5.12790928728471e-6
            RETURN
        CASE (14)
            out = -2.05116371491388e-6
            RETURN
        CASE (16)
            out = 3.22868362532741e-7
            RETURN
        END SELECT
    CASE (13)
        SELECT CASE (j)
        CASE (9)
            out = 1.51207581548139e-6
            RETURN
        CASE (11)
            out = -3.41860619152314e-6
            RETURN
        CASE (13)
            out = 3.07674557237082e-6
            RETURN
        CASE (15)
            out = -1.29147345013096e-6
            RETURN
        CASE (17)
            out = 2.11534444418003e-7
            RETURN
        END SELECT
    CASE (14)
        SELECT CASE (j)
        CASE (10)
            out = 8.54651547880785e-7
            RETURN
        CASE (12)
            out = -2.05116371491388e-6
            RETURN
        CASE (14)
            out = 1.93721017519645e-6
            RETURN
        CASE (16)
            out = -8.46137777672011e-7
            RETURN
        CASE (18)
            out = 1.43297526863808e-7
            RETURN
        END SELECT
    CASE (15)
        SELECT CASE (j)
        CASE (11)
            out = 5.12790928728471e-7
            RETURN
        CASE (13)
            out = -1.29147345013096e-6
            RETURN
        CASE (15)
            out = 1.26920666650802e-6
            RETURN
        CASE (17)
            out = -5.73190107455233e-7
            RETURN
        CASE (19)
            out = 9.98740338747754e-8
            RETURN
        END SELECT
    CASE (16)
        SELECT CASE (j)
        CASE (12)
            out = 3.22868362532741e-7
            RETURN
        CASE (14)
            out = -8.46137777672011e-7
            RETURN
        CASE (16)
            out = 8.59785161182849e-7
            RETURN
        CASE (18)
            out = -3.99496135499102e-7
            RETURN
        CASE (20)
            out = 7.13385956248396e-8
            RETURN
        END SELECT
    CASE (17)
        SELECT CASE (j)
        CASE (13)
            out = 2.11534444418003e-7
            RETURN
        CASE (15)
            out = -5.73190107455233e-7
            RETURN
        CASE (17)
            out = 5.99244203248653e-7
            RETURN
        CASE (19)
            out = -2.85354382499358e-7
            RETURN
        CASE (21)
            out = 5.20578941046127e-8
            RETURN
        END SELECT
    CASE (18)
        SELECT CASE (j)
        CASE (14)
            out = 1.43297526863808e-7
            RETURN
        CASE (16)
            out = -3.99496135499102e-7
            RETURN
        CASE (18)
            out = 4.28031573749038e-7
            RETURN
        CASE (20)
            out = -2.08231576418451e-7
            RETURN
        CASE (22)
            out = 3.87097161290710e-8
            RETURN
        END SELECT
    CASE (19)
        SELECT CASE (j)
        CASE (15)
            out = 9.98740338747754e-8
            RETURN
        CASE (17)
            out = -2.85354382499358e-7
            RETURN
        CASE (19)
            out = 3.12347364627676e-7
            RETURN
        CASE (21)
            out = -1.54838864516284e-7
            RETURN
        CASE (23)
            out = 2.92683219512488e-8
            RETURN
        END SELECT
    CASE (20)
        SELECT CASE (j)
        CASE (16)
            out = 7.13385956248396e-8
            RETURN
        CASE (18)
            out = -2.08231576418451e-7
            RETURN
        CASE (20)
            out = 2.32258296774426e-7
            RETURN
        CASE (22)
            out = -1.17073287804995e-7
            RETURN
        CASE (24)
            out = 2.24617354509584e-8
            RETURN
        END SELECT
    CASE (21)
        SELECT CASE (j)
        CASE (17)
            out = 5.20578941046127e-8
            RETURN
        CASE (19)
            out = -1.54838864516284e-7
            RETURN
        CASE (21)
            out = 1.75609931707493e-7
            RETURN
        CASE (23)
            out = -8.98469418038335e-8
            RETURN
        CASE (25)
            out = 1.74702386840787e-8
            RETURN
        END SELECT
    CASE (22)
        SELECT CASE (j)
        CASE (18)
            out = 3.87097161290710e-8
            RETURN
        CASE (20)
            out = -1.17073287804995e-7
            RETURN
        CASE (22)
            out = 1.34770412705750e-7
            RETURN
        CASE (24)
            out = -6.98809547363149e-8
            RETURN
        CASE (26)
            out = 1.37531666236364e-8
            RETURN
        END SELECT
    CASE (23)
        SELECT CASE (j)
        CASE (19)
            out = 2.92683219512488e-8
            RETURN
        CASE (21)
            out = -8.98469418038335e-8
            RETURN
        CASE (23)
            out = 1.04821432104472e-7
            RETURN
        CASE (25)
            out = -5.50126664945458e-8
            RETURN
        CASE (27)
            out = 1.09463979249351e-8
            RETURN
        END SELECT
    CASE (24)
        SELECT CASE (j)
        CASE (20)
            out = 2.24617354509584e-8
            RETURN
        CASE (22)
            out = -6.98809547363149e-8
            RETURN
        CASE (24)
            out = 8.25189997418187e-8
            RETURN
        CASE (26)
            out = -4.37855916997405e-8
            RETURN
        CASE (28)
            out = 8.80004539063412e-9
            RETURN
        END SELECT
    CASE (25)
        SELECT CASE (j)
        CASE (21)
            out = 1.74702386840787e-8
            RETURN
        CASE (23)
            out = -5.50126664945458e-8
            RETURN
        CASE (25)
            out = 6.56783875496108e-8
            RETURN
        CASE (27)
            out = -3.52001815625365e-8
            RETURN
        CASE (29)
            out = 7.13965946787297e-9
            RETURN
        END SELECT
    CASE (26)
        SELECT CASE (j)
        CASE (22)
            out = 1.37531666236364e-8
            RETURN
        CASE (24)
            out = -4.37855916997405e-8
            RETURN
        CASE (26)
            out = 5.28002723438047e-8
            RETURN
        CASE (28)
            out = -2.85586378714919e-8
            RETURN
        CASE (30)
            out = 5.84153956462334e-9
            RETURN
        END SELECT
    CASE (27)
        SELECT CASE (j)
        CASE (23)
            out = 1.09463979249351e-8
            RETURN
        CASE (25)
            out = -3.52001815625365e-8
            RETURN
        CASE (27)
            out = 4.28379568072378e-8
            RETURN
        CASE (29)
            out = -2.33661582584934e-8
            RETURN
        END SELECT
    CASE (28)
        SELECT CASE (j)
        CASE (24)
            out = 8.80004539063412e-9
            RETURN
        CASE (26)
            out = -2.85586378714919e-8
            RETURN
        CASE (28)
            out = 3.50492373877400e-8
            RETURN
        CASE (30)
            out = -1.92668322482314e-8
            RETURN
        END SELECT
    CASE (29)
        SELECT CASE (j)
        CASE (25)
            out = 7.13965946787297e-9
            RETURN
        CASE (27)
            out = -2.33661582584934e-8
            RETURN
        CASE (29)
            out = 2.89002483723470e-8
            RETURN
        END SELECT
    CASE (30)
        SELECT CASE (j)
        CASE (26)
            out = 5.84153956462334e-9
            RETURN
        CASE (28)
            out = -1.92668322482314e-8
            RETURN
        CASE (30)
            out = 2.40019011905933e-8
            RETURN
        END SELECT
    END SELECT
END SUBROUTINE integral_ff

SUBROUTINE integral_ffxi(i, j, x1t, x1r, x2t, x2r, y1t, y1r, y2t, y2r, out)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: i, j
    REAL*8, INTENT(IN) :: x1t, x1r, x2t, x2r, y1t, y1r, y2t, y2r
    REAL*8, INTENT(OUT) :: out
    out = 0
    SELECT CASE (i)
    CASE (1)
        SELECT CASE (j)
        CASE (1)
            out = -0.5*x1t*y1t
            RETURN
        CASE (2)
            out = 0.1*x1t*y1r
            RETURN
        CASE (3)
            out = 0.5*x1t*y2t
            RETURN
        CASE (4)
            out = -0.1*x1t*y2r
            RETURN
        CASE (5)
            out = 0.0857142857142857*x1t
            RETURN
        CASE (7)
            out = -0.00317460317460317*x1t
            RETURN
        END SELECT
    CASE (2)
        SELECT CASE (j)
        CASE (1)
            out = -0.1*x1r*y1t
            RETURN
        CASE (3)
            out = 0.1*x1r*y2t
            RETURN
        CASE (4)
            out = -0.0166666666666667*x1r*y2r
            RETURN
        CASE (5)
            out = 0.00952380952380952*x1r
            RETURN
        CASE (6)
            out = 0.00476190476190476*x1r
            RETURN
        CASE (7)
            out = -0.00158730158730159*x1r
            RETURN
        END SELECT
    CASE (3)
        SELECT CASE (j)
        CASE (1)
            out = -0.5*x2t*y1t
            RETURN
        CASE (2)
            out = -0.1*x2t*y1r
            RETURN
        CASE (3)
            out = 0.5*x2t*y2t
            RETURN
        CASE (4)
            out = 0.1*x2t*y2r
            RETURN
        CASE (5)
            out = -0.0857142857142857*x2t
            RETURN
        CASE (7)
            out = 0.00317460317460317*x2t
            RETURN
        END SELECT
    CASE (4)
        SELECT CASE (j)
        CASE (1)
            out = 0.1*x2r*y1t
            RETURN
        CASE (2)
            out = 0.0166666666666667*x2r*y1r
            RETURN
        CASE (3)
            out = -0.1*x2r*y2t
            RETURN
        CASE (5)
            out = 0.00952380952380952*x2r
            RETURN
        CASE (6)
            out = -0.00476190476190476*x2r
            RETURN
        CASE (7)
            out = -0.00158730158730159*x2r
            RETURN
        END SELECT
    CASE (5)
        SELECT CASE (j)
        CASE (1)
            out = -0.0857142857142857*y1t
            RETURN
        CASE (2)
            out = -0.00952380952380952*y1r
            RETURN
        CASE (3)
            out = 0.0857142857142857*y2t
            RETURN
        CASE (4)
            out = -0.00952380952380952*y2r
            RETURN
        CASE (6)
            out = 0.00634920634920635
            RETURN
        CASE (8)
            out = -0.000577200577200577
            RETURN
        END SELECT
    CASE (6)
        SELECT CASE (j)
        CASE (2)
            out = -0.00476190476190476*y1r
            RETURN
        CASE (4)
            out = 0.00476190476190476*y2r
            RETURN
        CASE (5)
            out = -0.00634920634920635
            RETURN
        CASE (7)
            out = 0.00173160173160173
            RETURN
        CASE (9)
            out = -0.000222000222000222
            RETURN
        END SELECT
    CASE (7)
        SELECT CASE (j)
        CASE (1)
            out = 0.00317460317460317*y1t
            RETURN
        CASE (2)
            out = 0.00158730158730159*y1r
            RETURN
        CASE (3)
            out = -0.00317460317460317*y2t
            RETURN
        CASE (4)
            out = 0.00158730158730159*y2r
            RETURN
        CASE (6)
            out = -0.00173160173160173
            RETURN
        CASE (8)
            out = 0.000666000666000666
            RETURN
        CASE (10)
            out = -0.000103600103600104
            RETURN
        END SELECT
    CASE (8)
        SELECT CASE (j)
        CASE (5)
            out = 0.000577200577200577
            RETURN
        CASE (7)
            out = -0.000666000666000666
            RETURN
        CASE (9)
            out = 0.000310800310800311
            RETURN
        CASE (11)
            out = -5.48471136706431e-5
            RETURN
        END SELECT
    CASE (9)
        SELECT CASE (j)
        CASE (6)
            out = 0.000222000222000222
            RETURN
        CASE (8)
            out = -0.000310800310800311
            RETURN
        CASE (10)
            out = 0.000164541341011929
            RETURN
        CASE (12)
            out = -3.17535921251092e-5
            RETURN
        END SELECT
    CASE (10)
        SELECT CASE (j)
        CASE (7)
            out = 0.000103600103600104
            RETURN
        CASE (9)
            out = -0.000164541341011929
            RETURN
        CASE (11)
            out = 9.52607763753275e-5
            RETURN
        CASE (13)
            out = -1.96569856012580e-5
            RETURN
        END SELECT
    CASE (11)
        SELECT CASE (j)
        CASE (8)
            out = 5.48471136706431e-5
            RETURN
        CASE (10)
            out = -9.52607763753275e-5
            RETURN
        CASE (12)
            out = 5.89709568037741e-5
            RETURN
        CASE (14)
            out = -1.28197732182118e-5
            RETURN
        END SELECT
    CASE (12)
        SELECT CASE (j)
        CASE (9)
            out = 3.17535921251092e-5
            RETURN
        CASE (11)
            out = -5.89709568037741e-5
            RETURN
        CASE (13)
            out = 3.84593196546353e-5
            RETURN
        CASE (15)
            out = -8.71744578838400e-6
            RETURN
        END SELECT
    CASE (13)
        SELECT CASE (j)
        CASE (10)
            out = 1.96569856012580e-5
            RETURN
        CASE (12)
            out = -3.84593196546353e-5
            RETURN
        CASE (14)
            out = 2.61523373651520e-5
            RETURN
        CASE (16)
            out = -6.13449888812208e-6
            RETURN
        END SELECT
    CASE (14)
        SELECT CASE (j)
        CASE (11)
            out = 1.28197732182118e-5
            RETURN
        CASE (13)
            out = -2.61523373651520e-5
            RETURN
        CASE (15)
            out = 1.84034966643662e-5
            RETURN
        CASE (17)
            out = -4.44222333277806e-6
            RETURN
        END SELECT
    CASE (15)
        SELECT CASE (j)
        CASE (12)
            out = 8.71744578838400e-6
            RETURN
        CASE (14)
            out = -1.84034966643662e-5
            RETURN
        CASE (16)
            out = 1.33266699983342e-5
            RETURN
        CASE (18)
            out = -3.29584311786759e-6
            RETURN
        END SELECT
    CASE (16)
        SELECT CASE (j)
        CASE (13)
            out = 6.13449888812208e-6
            RETURN
        CASE (15)
            out = -1.33266699983342e-5
            RETURN
        CASE (17)
            out = 9.88752935360277e-6
            RETURN
        CASE (19)
            out = -2.49685084686939e-6
            RETURN
        END SELECT
    CASE (17)
        SELECT CASE (j)
        CASE (14)
            out = 4.44222333277806e-6
            RETURN
        CASE (16)
            out = -9.88752935360277e-6
            RETURN
        CASE (18)
            out = 7.49055254060816e-6
            RETURN
        CASE (20)
            out = -1.92614208187067e-6
            RETURN
        END SELECT
    CASE (18)
        SELECT CASE (j)
        CASE (15)
            out = 3.29584311786759e-6
            RETURN
        CASE (17)
            out = -7.49055254060816e-6
            RETURN
        CASE (19)
            out = 5.77842624561201e-6
            RETURN
        CASE (21)
            out = -1.50967892903377e-6
            RETURN
        END SELECT
    CASE (19)
        SELECT CASE (j)
        CASE (16)
            out = 2.49685084686939e-6
            RETURN
        CASE (18)
            out = -5.77842624561201e-6
            RETURN
        CASE (20)
            out = 4.52903678710130e-6
            RETURN
        CASE (22)
            out = -1.20000120000120e-6
            RETURN
        END SELECT
    CASE (20)
        SELECT CASE (j)
        CASE (17)
            out = 1.92614208187067e-6
            RETURN
        CASE (19)
            out = -4.52903678710130e-6
            RETURN
        CASE (21)
            out = 3.60000360000360e-6
            RETURN
        CASE (23)
            out = -9.65854624391210e-7
            RETURN
        END SELECT
    CASE (21)
        SELECT CASE (j)
        CASE (18)
            out = 1.50967892903377e-6
            RETURN
        CASE (20)
            out = -3.60000360000360e-6
            RETURN
        CASE (22)
            out = 2.89756387317363e-6
            RETURN
        CASE (24)
            out = -7.86160740783543e-7
            RETURN
        END SELECT
    CASE (22)
        SELECT CASE (j)
        CASE (19)
            out = 1.20000120000120e-6
            RETURN
        CASE (21)
            out = -2.89756387317363e-6
            RETURN
        CASE (23)
            out = 2.35848222235063e-6
            RETURN
        CASE (25)
            out = -6.46398831310913e-7
            RETURN
        END SELECT
    CASE (23)
        SELECT CASE (j)
        CASE (20)
            out = 9.65854624391210e-7
            RETURN
        CASE (22)
            out = -2.35848222235063e-6
            RETURN
        CASE (24)
            out = 1.93919649393274e-6
            RETURN
        CASE (26)
            out = -5.36373498321821e-7
            RETURN
        END SELECT
    CASE (24)
        SELECT CASE (j)
        CASE (21)
            out = 7.86160740783543e-7
            RETURN
        CASE (23)
            out = -1.93919649393274e-6
            RETURN
        CASE (25)
            out = 1.60912049496546e-6
            RETURN
        CASE (27)
            out = -4.48802314922340e-7
            RETURN
        END SELECT
    CASE (25)
        SELECT CASE (j)
        CASE (22)
            out = 6.46398831310913e-7
            RETURN
        CASE (24)
            out = -1.60912049496546e-6
            RETURN
        CASE (26)
            out = 1.34640694476702e-6
            RETURN
        CASE (28)
            out = -3.78401951797267e-7
            RETURN
        END SELECT
    CASE (26)
        SELECT CASE (j)
        CASE (23)
            out = 5.36373498321821e-7
            RETURN
        CASE (25)
            out = -1.34640694476702e-6
            RETURN
        CASE (27)
            out = 1.13520585539180e-6
            RETURN
        CASE (29)
            out = -3.21284676054284e-7
            RETURN
        END SELECT
    CASE (27)
        SELECT CASE (j)
        CASE (24)
            out = 4.48802314922340e-7
            RETURN
        CASE (26)
            out = -1.13520585539180e-6
            RETURN
        CASE (28)
            out = 9.63854028162851e-7
            RETURN
        CASE (30)
            out = -2.74552359537297e-7
            RETURN
        END SELECT
    CASE (28)
        SELECT CASE (j)
        CASE (25)
            out = 3.78401951797267e-7
            RETURN
        CASE (27)
            out = -9.63854028162851e-7
            RETURN
        CASE (29)
            out = 8.23657078611891e-7
            RETURN
        END SELECT
    CASE (29)
        SELECT CASE (j)
        CASE (26)
            out = 3.21284676054284e-7
            RETURN
        CASE (28)
            out = -8.23657078611891e-7
            RETURN
        CASE (30)
            out = 7.08056085122503e-7
            RETURN
        END SELECT
    CASE (30)
        SELECT CASE (j)
        CASE (27)
            out = 2.74552359537297e-7
            RETURN
        CASE (29)
            out = -7.08056085122503e-7
            RETURN
        END SELECT
    END SELECT
END SUBROUTINE integral_ffxi

SUBROUTINE integral_ffxixi(i, j, x1t, x1r, x2t, x2r, y1t, y1r, y2t, y2r, out)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: i, j
    REAL*8, INTENT(IN) :: x1t, x1r, x2t, x2r, y1t, y1r, y2t, y2r
    REAL*8, INTENT(OUT) :: out
    out = 0
    SELECT CASE (i)
    CASE (1)
        SELECT CASE (j)
        CASE (1)
            out = -0.6*x1t*y1t
            RETURN
        CASE (2)
            out = -0.55*x1t*y1r
            RETURN
        CASE (3)
            out = 0.6*x1t*y2t
            RETURN
        CASE (4)
            out = -0.05*x1t*y2r
            RETURN
        CASE (6)
            out = 0.0285714285714286*x1t
            RETURN
        END SELECT
    CASE (2)
        SELECT CASE (j)
        CASE (1)
            out = -0.05*x1r*y1t
            RETURN
        CASE (2)
            out = -0.0666666666666667*x1r*y1r
            RETURN
        CASE (3)
            out = 0.05*x1r*y2t
            RETURN
        CASE (4)
            out = 0.0166666666666667*x1r*y2r
            RETURN
        CASE (5)
            out = -0.0333333333333333*x1r
            RETURN
        CASE (6)
            out = 0.0142857142857143*x1r
            RETURN
        END SELECT
    CASE (3)
        SELECT CASE (j)
        CASE (1)
            out = 0.6*x2t*y1t
            RETURN
        CASE (2)
            out = 0.05*x2t*y1r
            RETURN
        CASE (3)
            out = -0.6*x2t*y2t
            RETURN
        CASE (4)
            out = 0.55*x2t*y2r
            RETURN
        CASE (6)
            out = -0.0285714285714286*x2t
            RETURN
        END SELECT
    CASE (4)
        SELECT CASE (j)
        CASE (1)
            out = -0.05*x2r*y1t
            RETURN
        CASE (2)
            out = 0.0166666666666667*x2r*y1r
            RETURN
        CASE (3)
            out = 0.05*x2r*y2t
            RETURN
        CASE (4)
            out = -0.0666666666666667*x2r*y2r
            RETURN
        CASE (5)
            out = 0.0333333333333333*x2r
            RETURN
        CASE (6)
            out = 0.0142857142857143*x2r
            RETURN
        END SELECT
    CASE (5)
        SELECT CASE (j)
        CASE (2)
            out = -0.0333333333333333*y1r
            RETURN
        CASE (4)
            out = 0.0333333333333333*y2r
            RETURN
        CASE (5)
            out = -0.0380952380952381
            RETURN
        CASE (7)
            out = 0.00634920634920635
            RETURN
        END SELECT
    CASE (6)
        SELECT CASE (j)
        CASE (1)
            out = 0.0285714285714286*y1t
            RETURN
        CASE (2)
            out = 0.0142857142857143*y1r
            RETURN
        CASE (3)
            out = -0.0285714285714286*y2t
            RETURN
        CASE (4)
            out = 0.0142857142857143*y2r
            RETURN
        CASE (6)
            out = -0.0126984126984127
            RETURN
        CASE (8)
            out = 0.00288600288600289
            RETURN
        END SELECT
    CASE (7)
        SELECT CASE (j)
        CASE (5)
            out = 0.00634920634920635
            RETURN
        CASE (7)
            out = -0.00577200577200577
            RETURN
        CASE (9)
            out = 0.00155400155400155
            RETURN
        END SELECT
    CASE (8)
        SELECT CASE (j)
        CASE (6)
            out = 0.00288600288600289
            RETURN
        CASE (8)
            out = -0.00310800310800311
            RETURN
        CASE (10)
            out = 0.000932400932400932
            RETURN
        END SELECT
    CASE (9)
        SELECT CASE (j)
        CASE (7)
            out = 0.00155400155400155
            RETURN
        CASE (9)
            out = -0.00186480186480186
            RETURN
        CASE (11)
            out = 0.000603318250377074
            RETURN
        END SELECT
    CASE (10)
        SELECT CASE (j)
        CASE (8)
            out = 0.000932400932400932
            RETURN
        CASE (10)
            out = -0.00120663650075415
            RETURN
        CASE (12)
            out = 0.000412796697626419
            RETURN
        END SELECT
    CASE (11)
        SELECT CASE (j)
        CASE (9)
            out = 0.000603318250377074
            RETURN
        CASE (11)
            out = -0.000825593395252838
            RETURN
        CASE (13)
            out = 0.000294854784018871
            RETURN
        END SELECT
    CASE (12)
        SELECT CASE (j)
        CASE (10)
            out = 0.000412796697626419
            RETURN
        CASE (12)
            out = -0.000589709568037741
            RETURN
        CASE (14)
            out = 0.000217936144709600
            RETURN
        END SELECT
    CASE (13)
        SELECT CASE (j)
        CASE (11)
            out = 0.000294854784018871
            RETURN
        CASE (13)
            out = -0.000435872289419200
            RETURN
        CASE (15)
            out = 0.000165631469979296
            RETURN
        END SELECT
    CASE (14)
        SELECT CASE (j)
        CASE (12)
            out = 0.000217936144709600
            RETURN
        CASE (14)
            out = -0.000331262939958592
            RETURN
        CASE (16)
            out = 0.000128824476650564
            RETURN
        END SELECT
    CASE (15)
        SELECT CASE (j)
        CASE (13)
            out = 0.000165631469979296
            RETURN
        CASE (15)
            out = -0.000257648953301127
            RETURN
        CASE (17)
            out = 0.000102171136653895
            RETURN
        END SELECT
    CASE (16)
        SELECT CASE (j)
        CASE (14)
            out = 0.000128824476650564
            RETURN
        CASE (16)
            out = -0.000204342273307791
            RETURN
        CASE (18)
            out = 8.23960779466897e-5
            RETURN
        END SELECT
    CASE (17)
        SELECT CASE (j)
        CASE (15)
            out = 0.000102171136653895
            RETURN
        CASE (17)
            out = -0.000164792155893379
            RETURN
        CASE (19)
            out = 6.74149728654734e-5
            RETURN
        END SELECT
    CASE (18)
        SELECT CASE (j)
        CASE (16)
            out = 8.23960779466897e-5
            RETURN
        CASE (18)
            out = -0.000134829945730947
            RETURN
        CASE (20)
            out = 5.58581203742494e-5
            RETURN
        END SELECT
    CASE (19)
        SELECT CASE (j)
        CASE (17)
            out = 6.74149728654734e-5
            RETURN
        CASE (19)
            out = -0.000111716240748499
            RETURN
        CASE (21)
            out = 4.68000468000468e-5
            RETURN
        END SELECT
    CASE (20)
        SELECT CASE (j)
        CASE (18)
            out = 5.58581203742494e-5
            RETURN
        CASE (20)
            out = -9.36000936000936e-5
            RETURN
        CASE (22)
            out = 3.96000396000396e-5
            RETURN
        END SELECT
    CASE (21)
        SELECT CASE (j)
        CASE (19)
            out = 4.68000468000468e-5
            RETURN
        CASE (21)
            out = -7.92000792000792e-5
            RETURN
        CASE (23)
            out = 3.38049118536923e-5
            RETURN
        END SELECT
    CASE (22)
        SELECT CASE (j)
        CASE (20)
            out = 3.96000396000396e-5
            RETURN
        CASE (22)
            out = -6.76098237073847e-5
            RETURN
        CASE (24)
            out = 2.90879474089911e-5
            RETURN
        END SELECT
    CASE (23)
        SELECT CASE (j)
        CASE (21)
            out = 3.38049118536923e-5
            RETURN
        CASE (23)
            out = -5.81758948179822e-5
            RETURN
        CASE (25)
            out = 2.52095544211256e-5
            RETURN
        END SELECT
    CASE (24)
        SELECT CASE (j)
        CASE (22)
            out = 2.90879474089911e-5
            RETURN
        CASE (24)
            out = -5.04191088422512e-5
            RETURN
        CASE (26)
            out = 2.19913134311947e-5
            RETURN
        END SELECT
    CASE (25)
        SELECT CASE (j)
        CASE (23)
            out = 2.52095544211256e-5
            RETURN
        CASE (25)
            out = -4.39826268623894e-5
            RETURN
        CASE (27)
            out = 1.92984995416606e-5
            RETURN
        END SELECT
    CASE (26)
        SELECT CASE (j)
        CASE (24)
            out = 2.19913134311947e-5
            RETURN
        CASE (26)
            out = -3.85969990833213e-5
            RETURN
        CASE (28)
            out = 1.70280878308770e-5
            RETURN
        END SELECT
    CASE (27)
        SELECT CASE (j)
        CASE (25)
            out = 1.92984995416606e-5
            RETURN
        CASE (27)
            out = -3.40561756617541e-5
            RETURN
        CASE (29)
            out = 1.51003797745513e-5
            RETURN
        END SELECT
    CASE (28)
        SELECT CASE (j)
        CASE (26)
            out = 1.70280878308770e-5
            RETURN
        CASE (28)
            out = -3.02007595491027e-5
            RETURN
        CASE (30)
            out = 1.34530656173275e-5
            RETURN
        END SELECT
    CASE (29)
        SELECT CASE (j)
        CASE (27)
            out = 1.51003797745513e-5
            RETURN
        CASE (29)
            out = -2.69061312346551e-5
            RETURN
        END SELECT
    CASE (30)
        SELECT CASE (j)
        CASE (28)
            out = 1.34530656173275e-5
            RETURN
        CASE (30)
            out = -2.40739068941651e-5
            RETURN
        END SELECT
    END SELECT
END SUBROUTINE integral_ffxixi

SUBROUTINE integral_fxifxi(i, j, x1t, x1r, x2t, x2r, y1t, y1r, y2t, y2r, out)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: i, j
    REAL*8, INTENT(IN) :: x1t, x1r, x2t, x2r, y1t, y1r, y2t, y2r
    REAL*8, INTENT(OUT) :: out
    out = 0
    SELECT CASE (i)
    CASE (1)
        SELECT CASE (j)
        CASE (1)
            out = 0.6*x1t*y1t
            RETURN
        CASE (2)
            out = 0.05*x1t*y1r
            RETURN
        CASE (3)
            out = -0.6*x1t*y2t
            RETURN
        CASE (4)
            out = 0.05*x1t*y2r
            RETURN
        CASE (6)
            out = -0.0285714285714286*x1t
            RETURN
        END SELECT
    CASE (2)
        SELECT CASE (j)
        CASE (1)
            out = 0.05*x1r*y1t
            RETURN
        CASE (2)
            out = 0.0666666666666667*x1r*y1r
            RETURN
        CASE (3)
            out = -0.05*x1r*y2t
            RETURN
        CASE (4)
            out = -0.0166666666666667*x1r*y2r
            RETURN
        CASE (5)
            out = 0.0333333333333333*x1r
            RETURN
        CASE (6)
            out = -0.0142857142857143*x1r
            RETURN
        END SELECT
    CASE (3)
        SELECT CASE (j)
        CASE (1)
            out = -0.6*x2t*y1t
            RETURN
        CASE (2)
            out = -0.05*x2t*y1r
            RETURN
        CASE (3)
            out = 0.6*x2t*y2t
            RETURN
        CASE (4)
            out = -0.05*x2t*y2r
            RETURN
        CASE (6)
            out = 0.0285714285714286*x2t
            RETURN
        END SELECT
    CASE (4)
        SELECT CASE (j)
        CASE (1)
            out = 0.05*x2r*y1t
            RETURN
        CASE (2)
            out = -0.0166666666666667*x2r*y1r
            RETURN
        CASE (3)
            out = -0.05*x2r*y2t
            RETURN
        CASE (4)
            out = 0.0666666666666667*x2r*y2r
            RETURN
        CASE (5)
            out = -0.0333333333333333*x2r
            RETURN
        CASE (6)
            out = -0.0142857142857143*x2r
            RETURN
        END SELECT
    CASE (5)
        SELECT CASE (j)
        CASE (2)
            out = 0.0333333333333333*y1r
            RETURN
        CASE (4)
            out = -0.0333333333333333*y2r
            RETURN
        CASE (5)
            out = 0.0380952380952381
            RETURN
        CASE (7)
            out = -0.00634920634920635
            RETURN
        END SELECT
    CASE (6)
        SELECT CASE (j)
        CASE (1)
            out = -0.0285714285714286*y1t
            RETURN
        CASE (2)
            out = -0.0142857142857143*y1r
            RETURN
        CASE (3)
            out = 0.0285714285714286*y2t
            RETURN
        CASE (4)
            out = -0.0142857142857143*y2r
            RETURN
        CASE (6)
            out = 0.0126984126984127
            RETURN
        CASE (8)
            out = -0.00288600288600289
            RETURN
        END SELECT
    CASE (7)
        SELECT CASE (j)
        CASE (5)
            out = -0.00634920634920635
            RETURN
        CASE (7)
            out = 0.00577200577200577
            RETURN
        CASE (9)
            out = -0.00155400155400155
            RETURN
        END SELECT
    CASE (8)
        SELECT CASE (j)
        CASE (6)
            out = -0.00288600288600289
            RETURN
        CASE (8)
            out = 0.00310800310800311
            RETURN
        CASE (10)
            out = -0.000932400932400932
            RETURN
        END SELECT
    CASE (9)
        SELECT CASE (j)
        CASE (7)
            out = -0.00155400155400155
            RETURN
        CASE (9)
            out = 0.00186480186480186
            RETURN
        CASE (11)
            out = -0.000603318250377074
            RETURN
        END SELECT
    CASE (10)
        SELECT CASE (j)
        CASE (8)
            out = -0.000932400932400932
            RETURN
        CASE (10)
            out = 0.00120663650075415
            RETURN
        CASE (12)
            out = -0.000412796697626419
            RETURN
        END SELECT
    CASE (11)
        SELECT CASE (j)
        CASE (9)
            out = -0.000603318250377074
            RETURN
        CASE (11)
            out = 0.000825593395252838
            RETURN
        CASE (13)
            out = -0.000294854784018871
            RETURN
        END SELECT
    CASE (12)
        SELECT CASE (j)
        CASE (10)
            out = -0.000412796697626419
            RETURN
        CASE (12)
            out = 0.000589709568037741
            RETURN
        CASE (14)
            out = -0.000217936144709600
            RETURN
        END SELECT
    CASE (13)
        SELECT CASE (j)
        CASE (11)
            out = -0.000294854784018871
            RETURN
        CASE (13)
            out = 0.000435872289419200
            RETURN
        CASE (15)
            out = -0.000165631469979296
            RETURN
        END SELECT
    CASE (14)
        SELECT CASE (j)
        CASE (12)
            out = -0.000217936144709600
            RETURN
        CASE (14)
            out = 0.000331262939958592
            RETURN
        CASE (16)
            out = -0.000128824476650564
            RETURN
        END SELECT
    CASE (15)
        SELECT CASE (j)
        CASE (13)
            out = -0.000165631469979296
            RETURN
        CASE (15)
            out = 0.000257648953301127
            RETURN
        CASE (17)
            out = -0.000102171136653895
            RETURN
        END SELECT
    CASE (16)
        SELECT CASE (j)
        CASE (14)
            out = -0.000128824476650564
            RETURN
        CASE (16)
            out = 0.000204342273307791
            RETURN
        CASE (18)
            out = -8.23960779466897e-5
            RETURN
        END SELECT
    CASE (17)
        SELECT CASE (j)
        CASE (15)
            out = -0.000102171136653895
            RETURN
        CASE (17)
            out = 0.000164792155893379
            RETURN
        CASE (19)
            out = -6.74149728654734e-5
            RETURN
        END SELECT
    CASE (18)
        SELECT CASE (j)
        CASE (16)
            out = -8.23960779466897e-5
            RETURN
        CASE (18)
            out = 0.000134829945730947
            RETURN
        CASE (20)
            out = -5.58581203742494e-5
            RETURN
        END SELECT
    CASE (19)
        SELECT CASE (j)
        CASE (17)
            out = -6.74149728654734e-5
            RETURN
        CASE (19)
            out = 0.000111716240748499
            RETURN
        CASE (21)
            out = -4.68000468000468e-5
            RETURN
        END SELECT
    CASE (20)
        SELECT CASE (j)
        CASE (18)
            out = -5.58581203742494e-5
            RETURN
        CASE (20)
            out = 9.36000936000936e-5
            RETURN
        CASE (22)
            out = -3.96000396000396e-5
            RETURN
        END SELECT
    CASE (21)
        SELECT CASE (j)
        CASE (19)
            out = -4.68000468000468e-5
            RETURN
        CASE (21)
            out = 7.92000792000792e-5
            RETURN
        CASE (23)
            out = -3.38049118536923e-5
            RETURN
        END SELECT
    CASE (22)
        SELECT CASE (j)
        CASE (20)
            out = -3.96000396000396e-5
            RETURN
        CASE (22)
            out = 6.76098237073847e-5
            RETURN
        CASE (24)
            out = -2.90879474089911e-5
            RETURN
        END SELECT
    CASE (23)
        SELECT CASE (j)
        CASE (21)
            out = -3.38049118536923e-5
            RETURN
        CASE (23)
            out = 5.81758948179822e-5
            RETURN
        CASE (25)
            out = -2.52095544211256e-5
            RETURN
        END SELECT
    CASE (24)
        SELECT CASE (j)
        CASE (22)
            out = -2.90879474089911e-5
            RETURN
        CASE (24)
            out = 5.04191088422512e-5
            RETURN
        CASE (26)
            out = -2.19913134311947e-5
            RETURN
        END SELECT
    CASE (25)
        SELECT CASE (j)
        CASE (23)
            out = -2.52095544211256e-5
            RETURN
        CASE (25)
            out = 4.39826268623894e-5
            RETURN
        CASE (27)
            out = -1.92984995416606e-5
            RETURN
        END SELECT
    CASE (26)
        SELECT CASE (j)
        CASE (24)
            out = -2.19913134311947e-5
            RETURN
        CASE (26)
            out = 3.85969990833213e-5
            RETURN
        CASE (28)
            out = -1.70280878308770e-5
            RETURN
        END SELECT
    CASE (27)
        SELECT CASE (j)
        CASE (25)
            out = -1.92984995416606e-5
            RETURN
        CASE (27)
            out = 3.40561756617541e-5
            RETURN
        CASE (29)
            out = -1.51003797745513e-5
            RETURN
        END SELECT
    CASE (28)
        SELECT CASE (j)
        CASE (26)
            out = -1.70280878308770e-5
            RETURN
        CASE (28)
            out = 3.02007595491027e-5
            RETURN
        CASE (30)
            out = -1.34530656173275e-5
            RETURN
        END SELECT
    CASE (29)
        SELECT CASE (j)
        CASE (27)
            out = -1.51003797745513e-5
            RETURN
        CASE (29)
            out = 2.69061312346551e-5
            RETURN
        END SELECT
    CASE (30)
        SELECT CASE (j)
        CASE (28)
            out = -1.34530656173275e-5
            RETURN
        CASE (30)
            out = 2.40739068941651e-5
            RETURN
        END SELECT
    END SELECT
END SUBROUTINE integral_fxifxi

SUBROUTINE integral_fxifxixi(i, j, x1t, x1r, x2t, x2r, y1t, y1r, y2t, y2r, out)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: i, j
    REAL*8, INTENT(IN) :: x1t, x1r, x2t, x2r, y1t, y1r, y2t, y2r
    REAL*8, INTENT(OUT) :: out
    out = 0
    SELECT CASE (i)
    CASE (1)
        SELECT CASE (j)
        CASE (2)
            out = 0.25*x1t*y1r
            RETURN
        CASE (4)
            out = -0.25*x1t*y2r
            RETURN
        CASE (5)
            out = 0.2*x1t
            RETURN
        END SELECT
    CASE (2)
        SELECT CASE (j)
        CASE (1)
            out = -0.25*x1r*y1t
            RETURN
        CASE (2)
            out = -0.125*x1r*y1r
            RETURN
        CASE (3)
            out = 0.25*x1r*y2t
            RETURN
        CASE (4)
            out = -0.125*x1r*y2r
            RETURN
        CASE (5)
            out = 0.1*x1r
            RETURN
        END SELECT
    CASE (3)
        SELECT CASE (j)
        CASE (2)
            out = -0.25*x2t*y1r
            RETURN
        CASE (4)
            out = 0.25*x2t*y2r
            RETURN
        CASE (5)
            out = -0.2*x2t
            RETURN
        END SELECT
    CASE (4)
        SELECT CASE (j)
        CASE (1)
            out = 0.25*x2r*y1t
            RETURN
        CASE (2)
            out = 0.125*x2r*y1r
            RETURN
        CASE (3)
            out = -0.25*x2r*y2t
            RETURN
        CASE (4)
            out = 0.125*x2r*y2r
            RETURN
        CASE (5)
            out = 0.1*x2r
            RETURN
        END SELECT
    CASE (5)
        SELECT CASE (j)
        CASE (1)
            out = -0.2*y1t
            RETURN
        CASE (2)
            out = -0.1*y1r
            RETURN
        CASE (3)
            out = 0.2*y2t
            RETURN
        CASE (4)
            out = -0.1*y2r
            RETURN
        CASE (6)
            out = 0.0571428571428571
            RETURN
        END SELECT
    CASE (6)
        SELECT CASE (j)
        CASE (5)
            out = -0.0571428571428571
            RETURN
        CASE (7)
            out = 0.0317460317460317
            RETURN
        END SELECT
    CASE (7)
        SELECT CASE (j)
        CASE (6)
            out = -0.0317460317460317
            RETURN
        CASE (8)
            out = 0.0202020202020202
            RETURN
        END SELECT
    CASE (8)
        SELECT CASE (j)
        CASE (7)
            out = -0.0202020202020202
            RETURN
        CASE (9)
            out = 0.0139860139860140
            RETURN
        END SELECT
    CASE (9)
        SELECT CASE (j)
        CASE (8)
            out = -0.0139860139860140
            RETURN
        CASE (10)
            out = 0.0102564102564103
            RETURN
        END SELECT
    CASE (10)
        SELECT CASE (j)
        CASE (9)
            out = -0.0102564102564103
            RETURN
        CASE (11)
            out = 0.00784313725490196
            RETURN
        END SELECT
    CASE (11)
        SELECT CASE (j)
        CASE (10)
            out = -0.00784313725490196
            RETURN
        CASE (12)
            out = 0.00619195046439629
            RETURN
        END SELECT
    CASE (12)
        SELECT CASE (j)
        CASE (11)
            out = -0.00619195046439629
            RETURN
        CASE (13)
            out = 0.00501253132832080
            RETURN
        END SELECT
    CASE (13)
        SELECT CASE (j)
        CASE (12)
            out = -0.00501253132832080
            RETURN
        CASE (14)
            out = 0.00414078674948240
            RETURN
        END SELECT
    CASE (14)
        SELECT CASE (j)
        CASE (13)
            out = -0.00414078674948240
            RETURN
        CASE (15)
            out = 0.00347826086956522
            RETURN
        END SELECT
    CASE (15)
        SELECT CASE (j)
        CASE (14)
            out = -0.00347826086956522
            RETURN
        CASE (16)
            out = 0.00296296296296296
            RETURN
        END SELECT
    CASE (16)
        SELECT CASE (j)
        CASE (15)
            out = -0.00296296296296296
            RETURN
        CASE (17)
            out = 0.00255427841634738
            RETURN
        END SELECT
    CASE (17)
        SELECT CASE (j)
        CASE (16)
            out = -0.00255427841634738
            RETURN
        CASE (18)
            out = 0.00222469410456062
            RETURN
        END SELECT
    CASE (18)
        SELECT CASE (j)
        CASE (17)
            out = -0.00222469410456062
            RETURN
        CASE (19)
            out = 0.00195503421309873
            RETURN
        END SELECT
    CASE (19)
        SELECT CASE (j)
        CASE (18)
            out = -0.00195503421309873
            RETURN
        CASE (20)
            out = 0.00173160173160173
            RETURN
        END SELECT
    CASE (20)
        SELECT CASE (j)
        CASE (19)
            out = -0.00173160173160173
            RETURN
        CASE (21)
            out = 0.00154440154440154
            RETURN
        END SELECT
    CASE (21)
        SELECT CASE (j)
        CASE (20)
            out = -0.00154440154440154
            RETURN
        CASE (22)
            out = 0.00138600138600139
            RETURN
        END SELECT
    CASE (22)
        SELECT CASE (j)
        CASE (21)
            out = -0.00138600138600139
            RETURN
        CASE (23)
            out = 0.00125078173858662
            RETURN
        END SELECT
    CASE (23)
        SELECT CASE (j)
        CASE (22)
            out = -0.00125078173858662
            RETURN
        CASE (24)
            out = 0.00113442994895065
            RETURN
        END SELECT
    CASE (24)
        SELECT CASE (j)
        CASE (23)
            out = -0.00113442994895065
            RETURN
        CASE (25)
            out = 0.00103359173126615
            RETURN
        END SELECT
    CASE (25)
        SELECT CASE (j)
        CASE (24)
            out = -0.00103359173126615
            RETURN
        CASE (26)
            out = 0.000945626477541371
            RETURN
        END SELECT
    CASE (26)
        SELECT CASE (j)
        CASE (25)
            out = -0.000945626477541371
            RETURN
        CASE (27)
            out = 0.000868432479374729
            RETURN
        END SELECT
    CASE (27)
        SELECT CASE (j)
        CASE (26)
            out = -0.000868432479374729
            RETURN
        CASE (28)
            out = 0.000800320128051221
            RETURN
        END SELECT
    CASE (28)
        SELECT CASE (j)
        CASE (27)
            out = -0.000800320128051221
            RETURN
        CASE (29)
            out = 0.000739918608953015
            RETURN
        END SELECT
    CASE (29)
        SELECT CASE (j)
        CASE (28)
            out = -0.000739918608953015
            RETURN
        CASE (30)
            out = 0.000686106346483705
            RETURN
        END SELECT
    CASE (30)
        SELECT CASE (j)
        CASE (29)
            out = -0.000686106346483705
            RETURN
        END SELECT
    END SELECT
END SUBROUTINE integral_fxifxixi

SUBROUTINE integral_fxixifxixi(i, j, x1t, x1r, x2t, x2r, y1t, y1r, y2t, y2r, out)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: i, j
    REAL*8, INTENT(IN) :: x1t, x1r, x2t, x2r, y1t, y1r, y2t, y2r
    REAL*8, INTENT(OUT) :: out
    out = 0
    SELECT CASE (i)
    CASE (1)
        SELECT CASE (j)
        CASE (1)
            out = 1.5*x1t*y1t
            RETURN
        CASE (2)
            out = 0.75*x1t*y1r
            RETURN
        CASE (3)
            out = -1.5*x1t*y2t
            RETURN
        CASE (4)
            out = 0.75*x1t*y2r
            RETURN
        END SELECT
    CASE (2)
        SELECT CASE (j)
        CASE (1)
            out = 0.75*x1r*y1t
            RETURN
        CASE (2)
            out = 0.5*x1r*y1r
            RETURN
        CASE (3)
            out = -0.75*x1r*y2t
            RETURN
        CASE (4)
            out = 0.25*x1r*y2r
            RETURN
        END SELECT
    CASE (3)
        SELECT CASE (j)
        CASE (1)
            out = -1.5*x2t*y1t
            RETURN
        CASE (2)
            out = -0.75*x2t*y1r
            RETURN
        CASE (3)
            out = 1.5*x2t*y2t
            RETURN
        CASE (4)
            out = -0.75*x2t*y2r
            RETURN
        END SELECT
    CASE (4)
        SELECT CASE (j)
        CASE (1)
            out = 0.75*x2r*y1t
            RETURN
        CASE (2)
            out = 0.25*x2r*y1r
            RETURN
        CASE (3)
            out = -0.75*x2r*y2t
            RETURN
        CASE (4)
            out = 0.5*x2r*y2r
            RETURN
        END SELECT
    CASE (5)
        SELECT CASE (j)
        CASE (5)
            out = 0.400000000000000
            RETURN
        END SELECT
    CASE (6)
        SELECT CASE (j)
        CASE (6)
            out = 0.285714285714286
            RETURN
        END SELECT
    CASE (7)
        SELECT CASE (j)
        CASE (7)
            out = 0.222222222222222
            RETURN
        END SELECT
    CASE (8)
        SELECT CASE (j)
        CASE (8)
            out = 0.181818181818182
            RETURN
        END SELECT
    CASE (9)
        SELECT CASE (j)
        CASE (9)
            out = 0.153846153846154
            RETURN
        END SELECT
    CASE (10)
        SELECT CASE (j)
        CASE (10)
            out = 0.133333333333333
            RETURN
        END SELECT
    CASE (11)
        SELECT CASE (j)
        CASE (11)
            out = 0.117647058823529
            RETURN
        END SELECT
    CASE (12)
        SELECT CASE (j)
        CASE (12)
            out = 0.105263157894737
            RETURN
        END SELECT
    CASE (13)
        SELECT CASE (j)
        CASE (13)
            out = 0.0952380952380952
            RETURN
        END SELECT
    CASE (14)
        SELECT CASE (j)
        CASE (14)
            out = 0.0869565217391304
            RETURN
        END SELECT
    CASE (15)
        SELECT CASE (j)
        CASE (15)
            out = 0.0800000000000000
            RETURN
        END SELECT
    CASE (16)
        SELECT CASE (j)
        CASE (16)
            out = 0.0740740740740741
            RETURN
        END SELECT
    CASE (17)
        SELECT CASE (j)
        CASE (17)
            out = 0.0689655172413793
            RETURN
        END SELECT
    CASE (18)
        SELECT CASE (j)
        CASE (18)
            out = 0.0645161290322581
            RETURN
        END SELECT
    CASE (19)
        SELECT CASE (j)
        CASE (19)
            out = 0.0606060606060606
            RETURN
        END SELECT
    CASE (20)
        SELECT CASE (j)
        CASE (20)
            out = 0.0571428571428571
            RETURN
        END SELECT
    CASE (21)
        SELECT CASE (j)
        CASE (21)
            out = 0.0540540540540541
            RETURN
        END SELECT
    CASE (22)
        SELECT CASE (j)
        CASE (22)
            out = 0.0512820512820513
            RETURN
        END SELECT
    CASE (23)
        SELECT CASE (j)
        CASE (23)
            out = 0.0487804878048781
            RETURN
        END SELECT
    CASE (24)
        SELECT CASE (j)
        CASE (24)
            out = 0.0465116279069767
            RETURN
        END SELECT
    CASE (25)
        SELECT CASE (j)
        CASE (25)
            out = 0.0444444444444444
            RETURN
        END SELECT
    CASE (26)
        SELECT CASE (j)
        CASE (26)
            out = 0.0425531914893617
            RETURN
        END SELECT
    CASE (27)
        SELECT CASE (j)
        CASE (27)
            out = 0.0408163265306122
            RETURN
        END SELECT
    CASE (28)
        SELECT CASE (j)
        CASE (28)
            out = 0.0392156862745098
            RETURN
        END SELECT
    CASE (29)
        SELECT CASE (j)
        CASE (29)
            out = 0.0377358490566038
            RETURN
        END SELECT
    CASE (30)
        SELECT CASE (j)
        CASE (30)
            out = 0.0363636363636364
            RETURN
        END SELECT
    END SELECT
END SUBROUTINE integral_fxixifxixi

END MODULE BARDELL

max_stars_repo_name: saullocastro/structmanager

max_issues_repo_name: saullocastro/structmanager

max_forks_repo_name: saullocastro/structmanager


text:
      FUNCTION PIMACH (DUM)
C
C     THIS SUBPROGRAM SUPPLIES THE VALUE OF THE CONSTANT PI CORRECT TO
C     MACHINE PRECISION WHERE
C
C     PI=3.1415926535897932384626433832795028841971693993751058209749446
C
      PIMACH = 3.14159265358979
      RETURN
      END

ext: f

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
subroutine rupclo(ndime,ngaus,posgp,weigp)

!-----------------------------------------------------------------------
!
!     This routine sets up the integration constants of closed rules for
!     PRISMS
!
!             NDIME = 3
!
!          NGAUS  EXACT POL.
!          -----  ----------
!            3       p1
!
!-----------------------------------------------------------------------

  use      typre
  implicit none

  integer(ip), intent(in)  :: ndime,ngaus
  real(rp),    intent(out) :: posgp(ndime,ngaus),weigp(ngaus)
  integer(ip)              :: istop

!
! Area integral
!
  if(ndime==3) then
     if(ngaus==6) then
        posgp(1,1)= 0.0_rp
        posgp(2,1)= 0.0_rp
        posgp(3,1)= 0.0_rp
        posgp(1,2)= 1.0_rp
        posgp(2,2)= 0.0_rp
        posgp(3,2)= 0.0_rp
        posgp(1,3)= 0.0_rp
        posgp(2,3)= 1.0_rp
        posgp(3,3)= 0.0_rp
        posgp(1,4)= 0.0_rp
        posgp(2,4)= 0.0_rp
        posgp(3,4)= 1.0_rp
        posgp(1,5)= 1.0_rp
        posgp(2,5)= 0.0_rp
        posgp(3,5)= 1.0_rp
        posgp(1,6)= 0.0_rp
        posgp(2,6)= 1.0_rp
        posgp(3,6)= 1.0_rp
        weigp(  1)= 1.0_rp/12.0_rp
        weigp(  2)= 1.0_rp/12.0_rp
        weigp(  3)= 1.0_rp/12.0_rp
        weigp(  4)= 1.0_rp/12.0_rp
        weigp(  5)= 1.0_rp/12.0_rp
        weigp(  6)= 1.0_rp/12.0_rp
     else
        istop=1
     end if
  else if(ndime==2 .and. ngaus==0) then
  else
     istop=1
  end if
!
! Errors
!
  if(istop==1) call runend('RUTCLO: NOT AVAILABLE INTEGRATION RULE')

end subroutine rupclo

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
subroutine binmodel(npt,wv,nmodel,wmod,fmod,fmodbin,rv)
!Jason Rowe 2015 - jasonfrowe@gmail.com
use precision
implicit none
integer :: npt,nmodel,i,nbin
integer, allocatable, dimension(:) :: q,nfmodbin
real(double), parameter :: c=2.99792458d8
real(double) :: wvdiffmed,zpt,dl,rv
real(double), dimension(:) :: wv,wmod,fmod,fmodbin
real(double), allocatable, dimension(:) :: wvdiff

!calculate spacing between observed spectra.. assumes spectra is sorted short to long
allocate(q(npt-1),wvdiff(npt-1)) !allocate work arrays
do i=1,npt-1
   wvdiff(i)=wv(i+1)-wv(i) !calculate difference
enddo
call rqsort(npt-1,wvdiff,q) !sort differences
wvdiffmed=wvdiff(q((npt-1)/2)) !get median value
!write(0,*) "wvdiffmed: ",wvdiffmed
deallocate(q,wvdiff) !do not need work arrays
!write(0,*) "deallocated q,wvdiff"
zpt=wv(1)-wvdiffmed/2.0d0
!write(0,*) "zpt: ",zpt

!write(0,*) "allocating"
allocate(nfmodbin(npt))
!write(0,*) "done allocating"
nfmodbin=0
 fmodbin=0.0d0
do i=1,nmodel
!   write(0,*) i
   dl=wmod(i)*rv/c
!   write(0,*) dl
   nbin=int((wmod(i)-zpt+dl)/wvdiffmed)+1
!   write(0,*) dl,nbin
   if((nbin.gt.0).and.(nbin.le.npt))then
      nfmodbin(nbin)=nfmodbin(nbin)+1
       fmodbin(nbin)= fmodbin(nbin)+fmod(i)
   endif
enddo
do i=1,npt
   if(nfmodbin(i).gt.0)then
      fmodbin(i)=fmodbin(i)/dble(nfmodbin(i))
   else
      fmodbin(i)=0.0d0
   endif
enddo
!deallocate(nfmodbin)

return
end subroutine binmodel

max_stars_repo_name: AntoineDarveau/jwst-mtl

max_issues_repo_name: AntoineDarveau/jwst-mtl

max_forks_repo_name: AntoineDarveau/jwst-mtl


text:
\ CORE EXT

there

0 equ false                             \ 2.1000 -- false
-1 equ true                             \ 2.1010 -- true

: within over - >r - r> u< ;            \ 2.1020 x xlo xhi -- flag
: /string >r swap r@ + swap r> - ;      \ 2.1030 addr1 u1 n -- addr2 u2
: 0<>    0= 0= ; macro                  \ 2.1040 x y -- f
: 0>     negate 0< ;                    \ 2.1050 n -- f
: u>     swap u< ;                      \ 2.1060 u1 u2 -- flag
: 2>r        swap r> swap >r swap >r >r \ 2.1070 d -- | -- d
; no-tail-recursion
: 2r>        r> r> swap r> swap >r swap \ 2.1080 -- d | d --
; no-tail-recursion
: 2r@  r> r> r@ swap >r swap r@ swap >r \ 2.1090 -- d | d -- d
; no-tail-recursion
: third  >r >r dup r> swap r> swap ;    \ 2.1100 x1 x2 x3 -- x1 x2 x3 x1

\ lshift and rshift are hardly ever used so they are slow

: lshift  								\ 2.1110 x1 u -- x2
    dup if  for  2*  next  exit
    then drop
;

: rshift  								\ 2.1120 x1 u -- x2
    dup if  for  0 >carry 2/c  next  exit
    then drop
;

there swap - . .( instructions used by core ext) cr

max_stars_repo_name: uho/chad

max_issues_repo_name: uho/chad

max_forks_repo_name: uho/chad


      subroutine sd_channel_output (ichan)

      use sd_channel_module
      use basin_module
      use time_module
      use hydrograph_module
      use water_body_module

      implicit none

      integer, intent (in) :: ichan         !             |
      integer :: iob                        !             |
      real :: const                         !             |

      iob = sp_ob1%chandeg + ichan - 1

      !ch_stor_m(ichan) = ch_stor_m(ichan) + ch_stor(ichan)
      ch_in_m(ichan) = ch_in_m(ichan) + ch_in_d(ichan)
      ch_out_m(ichan) = ch_out_m(ichan) + ch_out_d(ichan)
      ch_wat_m(ichan) = ch_wat_m(ichan) + ch_wat_d(ichan)

!!!!! daily print
       if (pco%day_print == "y" .and. pco%int_day_cur == pco%int_day) then
        if (pco%sd_chan%d == "y") then
          write (2500,100) time%day, time%mo, time%day_mo, time%yrc, ichan, ob(iob)%gis_id, ob(iob)%name, ch_wat_d(ichan), &
            ch_stor(ichan), ch_in_d(ichan), ch_out_d(ichan), wtemp
           if (pco%csvout == "y") then
             write (2504,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, ichan, ob(iob)%gis_id, ob(iob)%name,  &
            ch_wat_d(ichan), ch_stor(ichan), ch_in_d(ichan), ch_out_d(ichan), wtemp
           end if
        end if
      end if

!!!!! monthly print
        if (time%end_mo == 1) then
          !ch_stor_y(ichan) = ch_stor_y(ichan) + ch_stor_m(ichan)
          ch_in_y(ichan) = ch_in_y(ichan) + ch_in_m(ichan)
          ch_out_y(ichan) = ch_out_y(ichan) + ch_out_m(ichan)
          ch_wat_y(ichan) = ch_wat_y(ichan) + ch_wat_m(ichan)

          const = float (ndays(time%mo + 1) - ndays(time%mo))
          !ch_stor_m(ichan) = ch_stor_m(ichan) / const           !! all storage variables are averages
          ch_wat_m(ichan) = ch_wat_m(ichan) // const            !! // only divides area (daily average values)

          if (pco%sd_chan%m == "y") then
          write (2501,100) time%day, time%mo, time%day_mo, time%yrc, ichan, ob(iob)%gis_id, ob(iob)%name, ch_wat_m(ichan), &
            ch_stor(ichan), ch_in_m(ichan), ch_out_m(ichan), wtemp
          if (pco%csvout == "y") then
            write (2505,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, ichan, ob(iob)%gis_id, ob(iob)%name,   &
            ch_wat_m(ichan), ch_stor(ichan), ch_in_m(ichan), ch_out_m(ichan), wtemp
          end if
        end if
        !ch_stor_m(ichan) = chaz
        ch_in_m(ichan) = chaz
        ch_out_m(ichan) = chaz
        ch_wat_m(ichan) = wbodz
        end if

!!!!! yearly print
      if (time%end_yr == 1) then
        !ch_stor_a(ichan) = ch_stor_a(ichan) + ch_stor_y(ichan)
        ch_in_a(ichan) = ch_in_a(ichan) + ch_in_y(ichan)
        ch_out_a(ichan) = ch_out_a(ichan) + ch_out_y(ichan)
        ch_wat_a(ichan) = ch_wat_a(ichan) + ch_wat_y(ichan)

        const = time%day_end_yr
        !ch_stor_y(ichan) = ch_stor_y(ichan) / const     !! all storage variables are averages
        ch_wat_y(ichan) = ch_wat_y(ichan) // const      !! // only divides area (daily average values)

        if (pco%sd_chan%y == "y") then
          write (2502,100) time%day, time%mo, time%day_mo, time%yrc, ichan, ob(iob)%gis_id, ob(iob)%name, ch_wat_y(ichan), &
            ch_stor(ichan), ch_in_y(ichan), ch_out_y(ichan), wtemp
          if (pco%csvout == "y") then
           write (2506,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, ichan, ob(iob)%gis_id, ob(iob)%name,    &
            ch_wat_y(ichan), ch_stor(ichan), ch_in_y(ichan), ch_out_y(ichan), wtemp
          end if
        end if
      end if

!!!!! average annual print
      if (time%end_sim == 1) then
        !ch_stor_a(ichan) = ch_stor_a(ichan) / time%yrs_prt      !! all storage variables (averaged) must be divided by years
        ch_in_a(ichan) = ch_in_a(ichan) / time%yrs_prt          !! all inflow and outflow varaibles (summed) are divided by years
        ch_out_a(ichan) = ch_out_a(ichan) / time%yrs_prt
        ch_wat_a(ichan) = ch_wat_a(ichan) / time%yrs_prt        !! all summed variables divided by years
        ch_wat_a(ichan) = ch_wat_a(ichan) // time%yrs_prt       !! all averaged variables divided by years

        if (pco%sd_chan%a == "y") then
        write (2503,100) time%day, time%mo, time%day_mo, time%yrc, ichan, ob(iob)%gis_id, ob(iob)%name, ch_wat_a(ichan), &
          ch_stor(ichan), ch_in_a(ichan), ch_out_a(ichan), wtemp
        if (pco%csvout == "y") then
          write (2507,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, ichan, ob(iob)%gis_id, ob(iob)%name,   &
          ch_wat_a(ichan), ch_stor(ichan), ch_in_a(ichan), ch_out_a(ichan), wtemp
        end if
       end if
     end if

      return

100   format (4i6,2i8,2x,a,61e15.4)

      end subroutine sd_channel_output


text:
C$Procedure      WNFETD ( Fetch an interval from a DP window )

      SUBROUTINE WNFETD ( WINDOW, N, LEFT, RIGHT )

C$ Abstract
C
C     Fetch a particular interval from a double precision window.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C      WINDOWS
C
C$ Keywords
C
C     WINDOWS
C
C$ Declarations

      INTEGER               LBCELL
      PARAMETER           ( LBCELL = -5 )

      DOUBLE PRECISION      WINDOW      ( LBCELL:* )
      INTEGER               N
      DOUBLE PRECISION      LEFT
      DOUBLE PRECISION      RIGHT

C$ Brief_I/O
C
C      VARIABLE  I/O  DESCRIPTION
C      --------  ---  --------------------------------------------------
C      WINDOW     I   Input window.
C      N          I   Index of interval to be fetched.
C      LEFT,
C      RIGHT      O   Left, right endpoints of the Nth interval.
C
C$ Detailed_Input
C
C      WINDOW      is a window containing zero or more intervals.
C
C      N           is the index of a particular interval within the
C                  window. Indices range from 1 to CARD(WINDOW)/2.
C
C$ Detailed_Output
C
C      LEFT,
C      RIGHT       are the left and right endpoints of the Nth interval
C                  in the input window. If the interval is not found,
C                  LEFT and RIGHT are not defined.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1) If N is less than one, the error SPICE(NOINTERVAL) signals.
C
C     2) If the interval does not exist, i.e. N > CARD(WINDOW)/2, the
C         error SPICE(NOINTERVAL) signals.
C
C$ Files
C
C      None.
C
C$ Particulars
C
C      None.
C
C$ Examples
C
C      Let A contain the intervals
C
C            [ 1, 3 ]  [ 7, 11 ]  [ 23, 27 ]
C
C      This window has a cardinality of 6, so N may have
C      value 1, 2, or 3 ( N =< CARD(WINDOW)/2 ).
C
C      Then the following calls
C
C            CALL WNFETD ( A,  1, LEFT, RIGHT )       [1]
C            CALL WNFETD ( A,  2, LEFT, RIGHT )       [2]
C            CALL WNFETD ( A,  3, LEFT, RIGHT )       [3]
C
C      yield the following values of LEFT and RIGHT
C
C            LEFT         RIGHT
C            ---------    ---------
C            1            3
C            7            11
C            23           27
C
C$ Restrictions
C
C      None.
C
C$ Literature_References
C
C      None.
C
C$ Author_and_Institution
C
C      W.L. Taber      (JPL)
C      I.M. Underwood  (JPL)
C
C$ Version
C
C-     SPICELIB Version 1.0.3, 30-JUL-2007 (EDW)
C
C         Removed erroneous description in the Examples section
C         indicating "Undefined" as a return state after an error
C         event caused by an invalid value of N.
C
C-     SPICELIB Version 1.0.1, 10-MAR-1992 (WLT)
C
C         Comment section for permuted index source lines was added
C         following the header.
C
C-     SPICELIB Version 1.0.0, 31-JAN-1990 (WLT) (IMU)
C
C-&

C$ Index_Entries
C
C     fetch an interval from a d.p. window
C
C-&

C
C     SPICELIB functions
C
      INTEGER               CARDD
      LOGICAL               RETURN

C
C     Local variables
C
      INTEGER               CARD
      INTEGER               END

C
C     Set up the error processing.
C
      IF ( RETURN () ) RETURN
      CALL CHKIN ( 'WNFETD' )

C
C
C     How many endpoints in the window? Enough? Normally, endpoints
C     of the Nth interval are stored in elements 2N and 2N-1.
C
      CARD = CARDD ( WINDOW )
      END  = 2*N

      IF ( N .LT. 1  .OR.  CARD .LT. END ) THEN
         CALL SETMSG ( 'WNFETD: No such interval.' )
         CALL SIGERR ( 'SPICE(NOINTERVAL)' )
      ELSE
         LEFT  = WINDOW(END-1)
         RIGHT = WINDOW(END  )
      END IF

      CALL CHKOUT ( 'WNFETD' )

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
!! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=83146
!! Fixed in 8.1.1 (20180521)
!!
!! ICE on select case statement with associate name
!!
!! $ gfortran --version
!! GNU Fortran (GCC) 8.0.0 20171123 (experimental)
!!
!! $ gfortran -c gfortran-20171124.f90
!! gfortran-20171124.f90:9:0:
!!
!!    select case (n_array(1))
!!
!! internal compiler error: in gfc_get_element_type, at fortran/trans-types.c:1231
!! 0x5a5da8 gfc_get_element_type(tree_node*)
!! 	../../gcc/fortran/trans-types.c:1231
!! 0x95d247 trans_associate_var
!! 	../../gcc/fortran/trans-stmt.c:1632
!! 0x95d247 gfc_trans_block_construct(gfc_code*)
!! 	../../gcc/fortran/trans-stmt.c:1890
!! 0x8e48c7 trans_code
!! 	../../gcc/fortran/trans.c:1924
!! 0x90e7a8 gfc_generate_function_code(gfc_namespace*)
!! 	../../gcc/fortran/trans-decl.c:6437
!! 0x89d036 translate_all_program_units
!! 	../../gcc/fortran/parse.c:6091
!! 0x89d036 gfc_parse_file()
!! 	../../gcc/fortran/parse.c:6294
!! 0x8e0eaf gfc_be_parse_file
!! 	../../gcc/fortran/f95-lang.c:204
!!

type foo
  integer n
end type
type bar
  type(foo) array(2)
end type
type(bar) b
associate (n_array => b%array%n)
  select case (n_array(1))
  case default
  end select
end associate
end

!! Turns out the associate block can be empty. This gives an ICE too
!type foo
!  integer n
!end type
!type bar
!  type(foo) array(2)
!end type
!type(bar) b
!associate (n_array => b%array%n)
!end associate
!end

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests