text:
! Program that blurs image data, by making every pixel the average of it's neighbours.
! Does not handle any actual image formats, so it just applies the bluring to dummy data.

module img_blur_module
    implicit none

    ! Custom pixel type, just has RGB values.
    type pixel
        integer :: r, g, b
    end type pixel

    ! Function declarations.
    public :: blur, set_pixel

contains

    function blur(img, n) result (blurred)
        integer, intent (in) :: n
        type (pixel), dimension (n, n), intent (in) :: img
        type (pixel), dimension (n, n) :: blurred

        integer :: i, j, cnt, r, g, b
        type (pixel) :: p

        do i = 1, n
            do j = 1, n
                cnt = 1
                r = img(i, j)%r
                g = img(i, j)%g
                b = img(i, j)%b

                ! Get average RGB values.
                if (i > 1) then
                    cnt = cnt + 1
                    r = r + img(i - 1, j)%r
                    g = g + img(i - 1, j)%g
                    b = b + img(i - 1, j)%b
                end if
                if (i < n) then
                    cnt = cnt + 1
                    r = r + img(i + 1, j)%r
                    g = g + img(i + 1, j)%g
                    b = b + img(i + 1, j)%b
                end if
                if (j > 1) then
                    cnt = cnt + 1
                    r = r + img(i, j - 1)%r
                    g = g + img(i, j - 1)%g
                    b = b + img(i, j - 1)%b
                end if
                if (j < n) then
                    cnt = cnt + 1
                    r = r + img(i, j + 1)%r
                    g = g + img(i, j + 1)%g
                    b = b + img(i, j + 1)%b
                end if

                ! Set the new, averaged value.
                call set_pixel(r / cnt, g / cnt, b / cnt, p)
                blurred(i, j) = p
            end do
        end do

    end function blur

    ! Sets a given pixel's RGB values.
    subroutine set_pixel(r, g, b, p)
        integer, intent (in) :: r, g, b
        type (pixel), intent (out) :: p

        p%r = r
        p%g = g
        p%b = b
    end subroutine set_pixel

end module img_blur_module

! Driver code.
program main
    use img_blur_module
    implicit none

    ! Set up a dummy image that is just half green and half red.
    integer, parameter :: n = 6
    type (pixel), dimension (n, n) :: img, blurred

    integer :: i, j
    type (pixel) :: green, red

    ! Set green and red pixel types.
    call set_pixel(0, 255, 0, green)
    call set_pixel(255, 0, 0, red)

    ! Color the dummy image.
    do i = 1, n / 2
        do j = 1, n
            img(i, j) = green
        end do
    end do
    do i = (n / 2) + 1, n
        do j = 1, n
            img(i, j) = red
        end do
    end do

    ! Blur it and print the results.
    blurred = blur(img, n)

    print '(a)', 'Original image:'
    do i = 1, n
        do j = 1, n
            write (*, '(a,i0,a,i0,a,i0,a,x)', advance='no') '(', img(j, i)%r, ',', img(j,i)%g, ',', img(j,i)%b, ')'
        end do
        print *, ''
    end do

    print '(a)', 'Blurred image:'
    do i = 1, n
        do j = 1, n
            write (*, '(a,i0,a,i0,a,i0,a,x)', advance='no') '(', blurred(j, i)%r, ',', blurred(j,i)%g, ',', blurred(j,i)%b, ')'
        end do
        print *, ''
    end do

end program main

max_stars_repo_name: ngergel/fortran-examples

max_issues_repo_name: ngergel/fortran-examples

max_forks_repo_name: ngergel/fortran-examples


text:
!
! COMMON_r3p3b_Module
!
! Module containing procedures for COMMON Record r3p3b.
!
!
! CREATION HISTORY:
!       Written by:   Paul van Delst, 24-Dec-2012
!                     paul.vandelst@noaa.gov
!

MODULE COMMON_r3p3b_Module

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module usage
  USE Type_Kinds     , ONLY: fp
  USE File_Utility   , ONLY: File_Open
  USE Message_Handler, ONLY: SUCCESS, FAILURE, INFORMATION, Display_Message
  ! Line-by-line model parameters
  USE LBL_Parameters
  ! Disable implicit typing
  IMPLICIT NONE

  ! ----------
  ! Visibility
  ! ----------
  ! Everything private by default
  PRIVATE
  ! Datatypes
  PUBLIC :: COMMON_r3p3b_type
  ! Procedures
  PUBLIC :: COMMON_r3p3b_Associated
  PUBLIC :: COMMON_r3p3b_Destroy
  PUBLIC :: COMMON_r3p3b_Create
  PUBLIC :: COMMON_r3p3b_SetValue
  PUBLIC :: COMMON_r3p3b_Write

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
  ! Message string length
  INTEGER, PARAMETER :: ML = 256
  ! The record I/O format
  CHARACTER(*), PARAMETER :: COMMON_R3P3B_FMT = '(8f10.3)'

  ! -------------
  ! Derived types
  ! -------------
  TYPE :: COMMON_r3p3b_type
    PRIVATE
    LOGICAL               :: is_allocated         = .FALSE.  ! Allocation indicator
    LOGICAL               :: is_pressure          = .FALSE.  ! Data type flag
    INTEGER               :: n_calculation_levels = 0        ! No. of calculation levels
    REAL(fp), ALLOCATABLE :: bnd(:)                          ! Altitudes/pressures of layer boundaries
  END TYPE COMMON_r3p3b_type

CONTAINS

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       COMMON_r3p3b_Associated
!
! PURPOSE:
!       Pure function to test the status of the allocatable components
!       of an Common record 3.3b object.
!
! CALLING SEQUENCE:
!       Status = COMMON_r3p3b_Associated( record )
!
! OBJECTS:
!       record:    Common record 3.3b structure which is to have its member's
!                  status tested.
!                  UNITS:      N/A
!                  TYPE:       COMMON_r3p3b_type
!                  DIMENSION:  Scalar
!                  ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       Status:    The return value is a logical value indicating the
!                  status of the object members.
!                    .TRUE.  - if the array components are allocated.
!                    .FALSE. - if the array components are not allocated.
!                  UNITS:      N/A
!                  TYPE:       LOGICAL
!                  DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  PURE FUNCTION COMMON_r3p3b_Associated( self ) RESULT( Status )
    TYPE(COMMON_r3p3b_type), INTENT(IN) :: self
    LOGICAL :: Status
    Status = self%is_allocated
  END FUNCTION COMMON_r3p3b_Associated

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       COMMON_r3p3b_Destroy
!
! PURPOSE:
!       Pure subroutine to re-initialize Common record 3.3b objects.
!
! CALLING SEQUENCE:
!       CALL COMMON_r3p3b_Destroy( record )
!
! OBJECTS:
!       record:    Re-initialized COMMON_r3p3b structure.
!                  UNITS:      N/A
!                  TYPE:       COMMON_r3p3b_type
!                  DIMENSION:  Scalar
!                  ATTRIBUTES: INTENT(OUT)
!
!:sdoc-:
!--------------------------------------------------------------------------------

  PUREL SUBROUTINE COMMON_r3p3b_Destroy( self )
    TYPE(COMMON_r3p3b_type), INTENT(OUT) :: self
    self%Is_Allocated = .FALSE.
  END SUBROUTINE COMMON_r3p3b_Destroy

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       COMMON_r3p3b_Create
!
! PURPOSE:
!       Pure subroutine to create an instance of the Common record 3.3b object.
!
! CALLING SEQUENCE:
!       CALL COMMON_r3p3b_Create( &
!              record                   , &
!              n_calculation_levels     , &
!              is_pressure = is_pressure  )
!
! OBJECTS:
!       record:                 Common record 3.3b structure.
!                               UNITS:      N/A
!                               TYPE:       COMMON_r3p3b_type
!                               DIMENSION:  Scalar
!                               ATTRIBUTES: INTENT(OUT)
!
! INPUTS:
!       n_calculation_levels:   The number of LBLRTM calculation LEVELS, i.e. the
!                               number of calculation level boundaries.
!                               If not specified, a default layering configuration
!                               is used.
!                               Must be > 1.
!                               UNITS:      N/A
!                               TYPE:       INTEGER
!                               DIMENSION:  Scalar
!                               ATTRIBUTES: INTENT(IN)
!
! OPTIONAL INPUTS:
!       is_pressure:            Logical flag to indicate the boundary values are
!                               specified as pressure levels.
!                               If == .FALSE., boundaries are specified as altitudes [DEFAULT]
!                                  == .TRUE.,  boundaries are specified as pressures.
!                               UNITS:      N/A
!                               TYPE:       LOGICAL
!                               DIMENSION:  Scalar
!                               ATTRIBUTES: INTENT(IN), OPTIONAL
!
!:sdoc-:
!--------------------------------------------------------------------------------

  PURE SUBROUTINE COMMON_r3p3b_Create( &
    self                , &  ! Output
    n_calculation_levels, &  ! Input
    is_pressure           )  ! Optional input
    ! Arguments
    TYPE(COMMON_r3p3b_type), INTENT(OUT) :: self
    INTEGER,                 INTENT(IN)  :: n_calculation_levels
    LOGICAL,       OPTIONAL, INTENT(IN)  :: is_pressure
    ! Local variables
    INTEGER :: alloc_stat

    ! Check input
    IF ( n_calculation_levels < 2 ) RETURN

    ! Perform the allocation
    ALLOCATE( self%bnd( n_calculation_levels ), &
              STAT = alloc_stat )
    IF ( alloc_stat /= 0 ) RETURN

    ! Initialise
    IF ( PRESENT(is_pressure) ) self%is_pressure = is_pressure
    self%n_calculation_levels = n_calculation_levels
    self%bnd                  = ZERO

    ! Set allocation indicator
    self%is_allocated = .TRUE.

  END SUBROUTINE COMMON_r3p3b_Create

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       COMMON_r3p3b_SetValue
!
! PURPOSE:
!       Subroutine to set the calculation layer boundaries in the
!       Common record 3.3b object.
!
! CALLING SEQUENCE:
!       CALL COMMON_r3p3b_SetValue( record, bnd )
!
! OBJECTS:
!       record:     Common record 3.3b structure which is to have its data set.
!                   - If unallocated, the object is allocated to the size of
!                     the input array.
!                   - If already allocated, the dimensions of the input array
!                     must be the same as that of the record allocated component
!                   UNITS:      N/A
!                   TYPE:       COMMON_r3p3b_type
!                   DIMENSION:  Scalar
!                   ATTRIBUTES: INTENT(IN OUT)
!
! INPUTS:
!       bnd:        Array of calculation level boundary altitudes or pressures,
!                   ordered from surface to top-of-atmosphere.
!                   UNITS:      km, or hPa
!                   TYPE:       REAL(fp)
!                   DIMENSION:  Rank-1. Size must be > 1
!                   ATTRIBUTES: INTENT(IN)
!
! COMMENTS:
!   If an error occurs, the record object is deallocated.
!
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE COMMON_r3p3b_SetValue( &
    self, &  ! In/output
    bnd   )  ! Input
    ! Arguments
    TYPE(COMMON_r3p3b_1D_type), INTENT(IN OUT) :: self
    REAL(fp)                  , INTENT(IN)     :: bnd(:)
    ! Local parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'COMMON_r3p3b_SetValue'
    ! Local variables
    CHARACTER(ML) :: msg
    INTEGER :: n_calculation_levels

    ! Set up
    n_calculation_levels = SIZE(bnd)

    ! Allocate structure if necessary
    IF ( .NOT. COMMON_r3p3b_Associated(self) ) THEN
      IF ( n_calculation_levels < 2 ) THEN
        msg = 'Must specify at least two level boundaries'
        CALL Display_Message( ROUTINE_NAME, msg, FAILURE ); RETURN
      END IF
      CALL COMMON_r3p3b_Create( self, n_calculation_levels, is_pressure = bnd(2) < bnd(1) )
      IF ( .NOT. COMMON_r3p3b_Associated(self) ) THEN
        msg = 'Allocation of COMMON_r3p3b structure failed'
        CALL Display_Message( ROUTINE_NAME, msg, FAILURE ); RETURN
      END IF
    END IF

    ! Check dimensions are consistent
    IF ( self%n_calculation_levels /= n_calculation_levels ) THEN
      WRITE( msg,'("Different size between structure (",i0,") and array (",i0,")")' ) &
                 self%n_calculation_levels, n_calculation_levels
      CALL Display_Message( ROUTINE_NAME, msg, FAILURE )
      CALL COMMON_r3p3b_Destroy( self )
      RETURN
    END IF

    ! Assign the boundary level data
    self%bnd = bnd

  END SUBROUTINE COMMON_r3p3b_SetValue

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       COMMON_r3p3b_Write
!
! PURPOSE:
!       Function to write an instance of the Common record 3.3b object.
!
! CALLING SEQUENCE:
!       error_status = COMMON_r3p3b_Write( record, file_id )
!
! OBJECTS:
!       record:       Common record 3.3b structure.
!                     UNITS:      N/A
!                     TYPE:       COMMON_r3p3b_type
!                     DIMENSION:  Scalar
!                     ATTRIBUTES: INTENT(IN)
!
! INPUTS:
!       file_id:      The unit number of the already opened file to which the
!                     record is to be written.
!                     UNITS:      N/A
!                     TYPE:       INTEGER
!                     DIMENSION:  Scalar
!                     ATTRIBUTES: INTENT(IN)
!
! FUNCTION RESULT:
!       Error_Status: The return value is an integer defining the error status.
!                     The error codes are defined in the Message_Handler module.
!                     If == SUCCESS, the write was successful
!                        == FAILURE, an unrecoverable error occurred.
!                     UNITS:      N/A
!                     TYPE:       INTEGER
!                     DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  FUNCTION COMMON_r3p3b_Write(self,fid) RESULT(err_stat)
    ! Arguments
    TYPE(COMMON_r3p3b_type), INTENT(IN) :: self
    INTEGER                , INTENT(IN) :: fid
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'COMMON_r3p3b_Write'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    INTEGER :: io_stat

    ! Setup
    err_stat = SUCCESS
    ! ...Check there is something to write
    IF ( .NOT. COMMON_r3p3b_Associated(self) ) THEN
      msg = 'Input data record is not allocated'
      CALL Cleanup(); RETURN
    END IF
    ! ...Check unit is open
    IF ( .NOT. File_Open(fid) ) THEN
      msg = 'File unit is not connected'
      CALL Cleanup(); RETURN
    END IF

    ! Write the record
    WRITE( fid,FMT=COMMON_R3P3B_FMT,IOSTAT=io_stat,IOMSG=io_msg) self%bnd
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing record - '//TRIM(io_msg)
      CALL Cleanup(); RETURN
    END IF

  CONTAINS

    SUBROUTINE CleanUp()
      err_stat = FAILURE
      CALL Display_Message( ROUTINE_NAME,msg,err_stat )
    END SUBROUTINE CleanUp

  END FUNCTION COMMON_r3p3b_Write

END MODULE COMMON_r3p3b_Module

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
C Copyright(C) 2011-2017 National Technology & Engineering Solutions of
C Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C * Redistributions of source code must retain the above copyright
C    notice, this list of conditions and the following disclaimer.
C
C * Redistributions in binary form must reproduce the above
C   copyright notice, this list of conditions and the following
C   disclaimer in the documentation and/or other materials provided
C   with the distribution.
C
C * Neither the name of NTESS nor the names of its
C   contributors may be used to endorse or promote products derived
C   from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C=======================================================================
      SUBROUTINE EXPXYZ (XN, YN, ZN, ICONOD, XEXPL, YEXPL, ZEXPL,
     &   MATMAP, NELBLK, IDELB, NUMELB, NUMLNK, LINK, NUMNP, NDIM,
     $     MODE)
C=======================================================================

C   --*** EXPXYZ *** (GREPOS) Modify coordinates for each block separately
C   --   Written by Amy Gilkey - revised 05/09/88
C   --   Modified by Greg Sjaardema - 02/06/89
C   --
C   --EXPXYZ modifies the coordinate array for the database.
C   --       each block is treated separately if not connected
C   --
C   --Parameters:
C   --   XN, YN, ZN - OUT - the coordinates
C   --   IDELB - IN - the element block IDs for each block
C   --   NUMELB - IN - the number of elements in each block
C   --   NUMLNK - IN - the number of nodes per element in each block
C   --   LINK -   IN - the connectivity for each block
C   --   NUMNP - IN - Number of nodes
C   --   NDIM  - IN - Number of spatial dimensions
C   --   MODE  - 1 = Explode
C   --           2 = Scale
C   --           3 = Randomize

      CHARACTER*80 STRING
      REAL XN(*), YN(*), ZN(*)
      INTEGER ICONOD(NUMNP, NELBLK)
      REAL XEXPL(*), YEXPL(*), ZEXPL(*)
      INTEGER MATMAP(NELBLK,NELBLK)
      INTEGER IDELB(*)
      INTEGER NUMELB(*)
      INTEGER NUMLNK(*)
      INTEGER LINK(*)
      INTEGER MODE

C   --Determine material block numbers of each node

      call iniint(numnp*nelblk, 0, iconod)

      IELNK = 0
      DO 10 IBLK = 1, NELBLK
         ISLNK = IELNK + 1
         IELNK = IELNK + NUMLNK(IBLK) * NUMELB(IBLK)

         CALL EXFCON (IBLK, NUMELB(IBLK), NUMLNK(IBLK), LINK(ISLNK),
     &      ICONOD(1,IBLK))
   10 CONTINUE

C ... Setup MATMAP which identifies connected material blocks.

      DO 30 IMAT = 1, NELBLK
         DO 20 IMAT2 = 1, NELBLK
            MATMAP(IMAT,IMAT2) = 0
   20    CONTINUE
         MATMAP(IMAT,IMAT) = 1
   30 CONTINUE

      DO 60 INOD = 1, NUMNP
         DO 50 IBLK = 1, NELBLK
            IF (ICONOD(INOD, IBLK) .NE. 0) THEN
               DO 40 IRST = IBLK+1, NELBLK
                  IF (ICONOD(INOD, IRST) .NE. 0) THEN
                     MATMAP(IBLK, IRST) = 1
                  END IF
   40          CONTINUE
            END IF
   50    CONTINUE
   60 CONTINUE

C ... If column J of row I of MATMAP is nonzero, then set
C     _EXPL(J) equal to _EXPL(I) and print message that material blocks
C     are connected.

      DO 80 IROW = 1, NELBLK
         DO 70 ICOL = IROW+1, NELBLK
            IF (MATMAP(IROW,ICOL) .NE. 0) THEN
               XEXPL(ICOL) = XEXPL(IROW)
               YEXPL(ICOL) = YEXPL(IROW)
               IF (NDIM .EQ. 3) THEN
                  ZEXPL(ICOL) = ZEXPL(IROW)
               END IF
               WRITE (STRING, 90) IDELB(IROW), IDELB(ICOL)
               CALL SQZSTR (STRING, LSTR)
               CALL PRTERR ('CMDSPEC', STRING(:LSTR))
            END IF
   70    CONTINUE
   80 CONTINUE
   90 FORMAT ('EXPXYZ -- Material ',I10,
     &        ' is connected to material ',I10)

C ... Collapse ICONOD array so first column is material ---
C      if multiple materials attached to one node, use first

      DO 120 INOD = 1, NUMNP
         IF (ICONOD(INOD, 1) .EQ. 0) THEN
            DO 100 IBLK = 2, NELBLK
               IF (ICONOD(INOD, IBLK) .NE. 0) THEN
                  ICONOD(INOD, 1) = IBLK
                  GO TO 110
               END IF
  100       CONTINUE
            CALL PRTERR ('ERROR', 'Node not connected to any elements')
  110       CONTINUE
         END IF
  120 CONTINUE

C ... Offset or scale each node
      IF (MODE .EQ. 1) THEN
         DO 130 INOD = 1, NUMNP
            IBLK = ICONOD(INOD,1)
            XN(INOD) = XN(INOD) + XEXPL(IBLK)
            YN(INOD) = YN(INOD) + YEXPL(IBLK)
            IF (NDIM .EQ. 3) THEN
               ZN(INOD) = ZN(INOD) + ZEXPL(IBLK)
            END IF
  130    CONTINUE
C ... Scale Each Node
      ELSE IF (MODE .EQ. 2) THEN
         DO 140 INOD = 1, NUMNP
            IBLK = ICONOD(INOD,1)
            XN(INOD) = XN(INOD) * XEXPL(IBLK)
            YN(INOD) = YN(INOD) * YEXPL(IBLK)
            IF (NDIM .EQ. 3) THEN
               ZN(INOD) = ZN(INOD) * ZEXPL(IBLK)
            END IF
  140    CONTINUE
C ... Randomize Each Node
      ELSE IF (MODE .EQ. 3) THEN
         IDUM = 1
         DO 150 INOD = 1, NUMNP
            IBLK = ICONOD(INOD,1)
            XN(INOD) = (2.0*RAN1(IDUM)-1.0) * XEXPL(IBLK) + XN(INOD)
            YN(INOD) = (2.0*RAN1(IDUM)-1.0) * YEXPL(IBLK) + YN(INOD)
            IF (NDIM .EQ. 3) THEN
               ZN(INOD) = (2.0*RAN1(IDUM)-1.0) * ZEXPL(IBLK) + ZN(INOD)
            END IF
  150    CONTINUE
      END IF
      RETURN
      END

ext: f

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


text:
! Calculates the integral defined in the paper by Drake and Yan from 1995.
!  Started conversion to Fortran on 06/23/2009
!  Updated 08/02/2009 to use HYGFX subroutine for calculating 2F1 hypergeometric function.
!  Modified 08/02/2009 to be used in Ps-H integral calculations.
!  Modified 08/16/2009 to use precomputed values for the W function through WMatrix.
!  Modified 08/19/2009 to skip asymptotic expansion for any j_{mu nu} even.
!  Modified 09/05/2009 to move 1/(2q+1)**2 out of inner loop in T.
!  Somewhere in here, I got this to produce the same as PVR's code and then made some updates.
!  Updated 09/22/2009 to use ModifiedRiemannZeta from Drake and Yan '97 codebase.

! TODO: Rewrite array limits (TValues, etc.)
! TODO: Do we really need to reallocate TValues every time?
! TODO: Do we need pmax in HylleraasIntegral?

module FactPreComp
   implicit none
   real*16, dimension(0:401) :: FactArray  !@TODO: Any possible need to even go this high?
end module FactPreComp

real*16 function HylleraasIntegral(j1, j2, j3, j12, j23, j31, alpha, beta, gamma, qmax, WMatrix, Method)
	use WLimits
	implicit none
	common /globals/ PI, HyperB, HyperC, HyperVal
	real*16 PI, HyperB, HyperC, HyperVal
	external T, AsymptoticPart
	real*16 T, AsymptoticPart
	integer j1, j2, j3, j12, j23, j31, qmax, Method
	real*16 alpha, beta, gamma
	real*16, dimension(lmin:lmax, mmin:mmax, nmin:nmax, 6) :: WMatrix
	integer q, pmax, LambdaUpper, ZetaUpperBound
	real*16 LambdaLower
	real*16 Sum
	real*16, allocatable, dimension(:) :: TValues  ! We do not know the size of TValues yet (N+1).
	real*16 HylleraasIntegralRecursion

	! This is for the recursion relations from Pachucki's paper (an analytic solution).
	if (Method == 2) then
		HylleraasIntegral = HylleraasIntegralRecursion(j1, j2, j3, j12, j23, j31, &
				real(alpha,16), real(beta,16), real(gamma,16), qmax, Method)
		return
	endif

	! Initialize global variables
	PI = 3.141592653589793238462643383279502884197q0
	HyperB = 0.0q0
	HyperC = 0.0q0
	HyperVal = 0.0q0

	! Initialize local variables
	pmax = 75
	!pmax = 65
	LambdaUpper = 6
	ZetaUpperBound = 500
	Sum = 0

	! Allocate the TValues array
	allocate(TValues(qmax+1)) ! We have N+1 possibilities (0 to N).

	do q = 0, qmax, 1
		! We store these to calculate the asymptotic part.
		TValues(q+1) = T(q, j1, j2, j3, j12, j23, j31, alpha, beta, gamma, pmax, WMatrix)
		! If the q summation terminates, then every term past this will be 0 as well.
		if (TValues(q+1) == 0.0q0) then
			exit
		endif
		!write (*,*) 'T(', q, '): ', TValues(q+1)
	enddo

	! Do the direct summation
	Sum = 0.0q0
	do q = 0, qmax, 1
		! No need to sum more terms if this summation was finite.  Exiting early also
		!  allows us to do the check right after this do loop.
		if (TValues(q+1) == 0.0q0) then
			exit
		endif
		Sum = Sum + TValues(q+1)
	enddo
	!write (*,*) 'N =', N, ':'
	! Output the direct sum.
	!write (*,*) 'Direct sum:', Sum

	! If the summation was finite, there is no need to do the asymptotic approximation.
	!  AsymptoticPart will return exactly 0 if it is called for a finite summation.
	if (q < qmax) then
		HylleraasIntegral = Sum
		return
	endif

	! Equation (13) for lowercase lambda
	LambdaLower = 0.5q0 * real(j12 + 1,16) + 0.5q0 * real(j23 + 1,16) + 0.5q0 * real(j31 + 1,16) + 4.0q0

	! Add on the asymptotic part
	if (Method == 1) then
		Sum = Sum + AsymptoticPart(TValues, LambdaLower, qmax+1, LambdaUpper)
	endif

	! Output the final result, including the asymptotic part.
	!write (*,*) 'Final sum:', Sum
	!write (*,*)

	! Deallocate the TValues array before we end this routine.
	deallocate(TValues)

	HylleraasIntegral = Sum

	return
end

! This is a standard factorial function (not using recursion).
!  Only integer arguments are accepted.
! To keep the output of Calc3j correct, this needs to return 0 whenever we have a negative value for n.
! Note: Arguments past 400 will cause a crash.
real*16 function Factorial(n)
	use FactPreComp
	implicit none
	integer n
	integer i

	if (n < 0) then
		Factorial = 0.0q0
		!write (*,*) 'n < 0 in factorial function'
		return
	endif

	Factorial = FactArray(n)
	return

	!Factorial = 1.0q0
	!do i = 1, n
	!	Factorial = Factorial * real(i, 16)
	!enddo
	!return
end

subroutine GenFactorial(MaxN)
	use FactPreComp
	implicit none
	real*16 Fact
	integer i, MaxN

	Fact = 1.0_16
	FactArray(0) = Fact
	do i = 1, MaxN
		Fact = Fact * real(i,16)
		FactArray(i) = Fact
	enddo

	return
end

! This is your standard "choose" function, or the binomial coefficient.
real*16 function Binomial(n, k)
	integer n, k
	real*16 Factorial

	Binomial = Factorial(n) / (Factorial(k) * Factorial(n-k))
	return
end

! S is given in equation (4) - just a simple minimum function.
! TODO: Rewrite as returning an integer as in equation (13) of Perkins.
real*16 function S(q, j)
	integer q, j
	S = MIN(real(q-1, 16), 0.5q0 * real(j+1, 16))
	return
end

! C is given in equation (4).
real*16 function C(j, q, k)
	integer j, q, k
	real*16 S, Binomial
	real*16 Multiplication, SVal, Bin
	integer t

	Bin = Binomial(j+2, 2*k + 1)
	Multiplication = 1.0q0
	SVal = S(q, j)

	do t = 0, int(SVal), 1
		Multiplication = Multiplication * (2*k + 2*t - j) / real(2*k + 2*q - 2*t + 1, 16)
	enddo

	C = (2*q + 1) / real(j + 2, 16) * Multiplication * Bin
	return
end

! This is the W function from equation (8).
real*16 function W(l, m, n, alpha, beta, gamma, pmax)
	implicit none
	integer l, m, n, b, c
	real*16 alpha, beta, gamma
	real*16 z, Hypergeometric2F1
	integer pmax
	real*16 Factorial, WPartial
	real*16 LeadingFactor, Summation
	real*16 Hyper, HyperTest, Diff
	integer p
	integer, save :: LMin1 = 1000000, LMax1 = -1000000, MMin1 = 1000000, MMax1 = -1000000, NMin1 = 1000000, NMax1 = -1000000

	if (l < LMin1) LMin1 = l
	if (l > LMax1) LMax1 = l
	if (m < MMin1) MMin1 = m
	if (m > MMax1) MMax1 = m
	if (n < NMin1) NMin1 = n
	if (n > NMax1) NMax1 = n
	write (*,*) LMin1, LMax1, MMin1, MMax1, NMin1, NMax1

	LeadingFactor = Factorial(l) / ((alpha + beta + gamma) ** (l+m+n+3))
	Summation = 0.0q0

	! We do the summation backwards so that we can use the backwards recurrence relation in (9).
	!do p = 0, pmax, 1
	b = l+m+n+pmax+3
	c = l+m+pmax+3
	z = real((alpha + beta) / (alpha + beta + gamma),16)
	Hyper = Hypergeometric2F1(1.0_16, real(b,16), real(c,16), real(z,16))
	do p = pmax, 0, -1
		!write (*,*) p
		Summation = Summation + WPartial(l, m, n, alpha, beta, gamma, p, Hyper)
		b = b - 1;  c = c - 1;
		Hyper = 1.0_16 + real(b,16) / real(c,16) * z * Hyper
		!HyperTest = Hypergeometric2F1(1.0_16, real(l+m+n+p-1+3,16), real(l+m+p-1+3,16), real(z,16))
		!!HyperTest = Hypergeometric2F1(1.0_16, real(b+1,16), real(c+1,16), z)
		!Diff = abs((Hyper-HyperTest)/(Hyper+HyperTest)*2)
		!if (Diff > 1e-13) then
		!	write (*,*) Hyper, HyperTest, abs(Hyper-HyperTest), ":", b, c, z
		!	Hyper = HyperTest
		!end if
	enddo

	W = LeadingFactor * Summation

	return
end

! This is what's inside the summation of the W function from equation (8).
real*16 function WPartial(l, m, n, alpha, beta, gamma, p, hyper)
	integer l, m, n
	real*16 alpha, beta, gamma
	integer p
	real*16 Factorial
	real*16 Hypergeometric2F1, hyper

	WPartial = Factorial(l+m+n+p+2) / (Factorial(l+1+p) * (l+m+2+p)) * (alpha / (alpha + beta + gamma)) ** p  * &
		hyper
		!Hypergeometric2F1(1.0_16, real(l+m+n+p+3,16), real(l+m+p+3,16), real((alpha + beta) / (alpha + beta + gamma),16) )

	return
end

! The T function from equation (6)
real*16 function T(q, j1, j2, j3, j12, j23, j31, alpha, beta, gamma, pmax, WMatrix)
	use WLimits
	implicit none
	integer q, j1, j2, j3, j12, j23, j31
	real*16 alpha, beta, gamma
	integer pmax
	real*16, dimension(lmin:lmax, mmin:mmax, nmin:nmax, 6) :: WMatrix

	common /globals/ PI, HyperB, HyperC, HyperVal
	real*16 PI, HyperB, HyperC, HyperVal
	real*16 W, C
	real*16 Summation
	real*16 L12, L23, L31
	integer k12, k23, k31

	real*16 RemoveMe1, RemoveMe2
	real*16 Res1,Res2,Res3,Res4,Res5,Res6,Res
	real*16 Temp
real*16 Rqt, Piterm
real*16 Rcoeff, C12, C23, C31

	Summation = 0.0q0
	! TODO: L12, L23 and L31 are always even, because of the if statements following.
!	L12 = .5 * (j12 + 1)  ! Need to convert each parentheses to real?
!	L23 = .5 * (j23 + 1)
!	L31 = .5 * (j31 + 1)
	L12 = (j12 + 1) / 2  ! Need to convert each parentheses to real?
	L23 = (j23 + 1) / 2
	L31 = (j31 + 1) / 2

	! If any of j12, j23 or j31 are even, then it changes the upper limit on its
	!  corresponding k12, k23 or k31 summation.  This also terminates the q summation
	!  as mentioned on page 1988 of the paper by J.F. Perkins (1967).
	if (mod(j12,2) == 0) then
		L12 = j12/2 - q
		if (q > j12 / 2) then
			T = 0.0q0
			return
		endif
	endif
	if (mod(j23,2) == 0) then
		L23 = j23/2 - q
		if (q > j23 / 2) then
			T = 0.0q0
			return
		endif
	endif
	if (mod(j31,2) == 0) then
		L31 = j31/2 - q
		if (q > j31 / 2) then
			T = 0.0q0
			return
		endif
	endif

	Rqt = (4.0q0*PI)**3 / ((2.0q0*q + 1.0q0)**2.0q0)

	do k12 = 0, int(L12), 1  ! TODO: Determine if L12, L23 and L31 can be real.
		do k31 = 0, int(L31), 1
			do k23 = 0, int(L23), 1

!		do k31 = 0, int(L31), 1
!			do k23 = 0, int(L23), 1

!					Summation = Summation + 1.0q0 / dble(2*q + 1) ** 2.0q0 * C(j12, q, k12) * C(j23, q, k23) * C(j31, q, k31) * &
!					( W(j1 + 2*q + 2*k12 + 2*k31 + 2, j2 + j12 - 2*k12 + 2*k23 + 2, j3 + j23 - 2*q - 2*k23 + j31 - 2*k31 + 2, &
!						alpha, beta, gamma, pmax) &
!					+ W(j1 + 2*q + 2*k12 + 2*k31 + 2, j3 + j31 - 2*k31 + 2*k23 + 2, j2 + j12 - 2*q - 2*k12 + j23 - 2*k23 + 2, &
!						alpha, gamma, beta, pmax) &
!					+ W(j2 + 2*q + 2*k12 + 2*k23 + 2, j1 + j12 - 2*k12 + 2*k31 + 2, j3 + j23 - 2*q - 2*k23 + j31 - 2*k31 + 2, &
!						beta, alpha, gamma, pmax) &
!					+ W(j2 + 2*q + 2*k12 + 2*k23 + 2, j3 + j23 - 2*k23 + 2*k31 + 2, j1 + j12 - 2*q - 2*k12 + j31 - 2*k31 + 2, &
!						beta, gamma, alpha, pmax) &
!					+ W(j3 + 2*q + 2*k23 + 2*k31 + 2, j1 + j31 - 2*k31 + 2*k12 + 2, j2 + j12 - 2*q - 2*k12 + j23 - 2*k23 + 2, &
!						gamma, alpha, beta, pmax) &
!					+ W(j3 + 2*q + 2*k23 + 2*k31 + 2, j2 + j23 - 2*k23 + 2*k12 + 2, j1 + j12 - 2*q - 2*k12 + j31 - 2*k31 + 2, &
!						gamma, beta, alpha, pmax))

				! The original version is preserved above.  Now this uses precalculated values for the W function via WMatrix.
				!  The ordering is the same as above, so that the 3 corresponds to beta, alpha, gamma and so on.
				Res1 = WMatrix(j1 + 2*q + 2*k12 + 2*k31 + 2, j2 + j12 - 2*k12 + 2*k23 + 2, j3 + j23 - 2*q - 2*k23 + j31 - 2*k31 + 2, 1)
				Res2 = WMatrix(j1 + 2*q + 2*k12 + 2*k31 + 2, j3 + j31 - 2*k31 + 2*k23 + 2, j2 + j12 - 2*q - 2*k12 + j23 - 2*k23 + 2, 2)
				Res3 = WMatrix(j2 + 2*q + 2*k12 + 2*k23 + 2, j1 + j12 - 2*k12 + 2*k31 + 2, j3 + j23 - 2*q - 2*k23 + j31 - 2*k31 + 2, 3)
				Res4 = WMatrix(j2 + 2*q + 2*k12 + 2*k23 + 2, j3 + j23 - 2*k23 + 2*k31 + 2, j1 + j12 - 2*q - 2*k12 + j31 - 2*k31 + 2, 4)
				Res5 = WMatrix(j3 + 2*q + 2*k23 + 2*k31 + 2, j1 + j31 - 2*k31 + 2*k12 + 2, j2 + j12 - 2*q - 2*k12 + j23 - 2*k23 + 2, 5)
				Res6 = WMatrix(j3 + 2*q + 2*k23 + 2*k31 + 2, j2 + j23 - 2*k23 + 2*k12 + 2, j1 + j12 - 2*q - 2*k12 + j31 - 2*k31 + 2, 6)

				C12 = C(j12, q, k12)
				C23 = C(j23, q, k23)
				C31 = C(j31, q, k31)
				Rcoeff = C12 * C31 * C23 !* Rqt
				Summation = Summation + Rcoeff * (Res1+Res2+Res3+Res4+Res5+Res6)

				! For some reason, this does not give the same results as above.  If I cast each of the WMatrix elements to real*16,
				!  I do get the correct answer.  I think I may be running out of registers when it is accessing the elements of the
				!  WMatrix, so it is not keeping the intermediate values in a larger-sized (80/96-bit) register.
!				Summation = Summation + &! C(j12, q, k12) * C(j23, q, k23) * C(j31, q, k31) * &
!				( WMatrix(j1 + 2*q + 2*k12 + 2*k31 + 2, j2 + j12 - 2*k12 + 2*k23 + 2, j3 + j23 - 2*q - 2*k23 + j31 - 2*k31 + 2, 1) &
!				+ WMatrix(j1 + 2*q + 2*k12 + 2*k31 + 2, j3 + j31 - 2*k31 + 2*k23 + 2, j2 + j12 - 2*q - 2*k12 + j23 - 2*k23 + 2, 2) &
!				+ WMatrix(j2 + 2*q + 2*k12 + 2*k23 + 2, j1 + j12 - 2*k12 + 2*k31 + 2, j3 + j23 - 2*q - 2*k23 + j31 - 2*k31 + 2, 3) &
!				+ WMatrix(j2 + 2*q + 2*k12 + 2*k23 + 2, j3 + j23 - 2*k23 + 2*k31 + 2, j1 + j12 - 2*q - 2*k12 + j31 - 2*k31 + 2, 4) &
!				+ WMatrix(j3 + 2*q + 2*k23 + 2*k31 + 2, j1 + j31 - 2*k31 + 2*k12 + 2, j2 + j12 - 2*q - 2*k12 + j23 - 2*k23 + 2, 5) &
!				+ WMatrix(j3 + 2*q + 2*k23 + 2*k31 + 2, j2 + j23 - 2*k23 + 2*k12 + 2, j1 + j12 - 2*q - 2*k12 + j31 - 2*k31 + 2, 6))
			enddo
		enddo
	enddo

	T = ((4.0q0 * PI) ** 3) / real(2*q + 1, 16) ** 2 * Summation
	!T = Summation
	!T = Summation / dble(2*q + 1)
	!T = Summation
	return
end

!! TODO: Description here
!real*16 function AsymptoticPart(TValues, LambdaLower, N, LambdaUpper, ZetaUpperBound)
!	integer N
!	real*16 TValues(N)
!	real*16 LambdaLower
!	integer LambdaUpper, ZetaUpperBound
!
!	external DGESV, ModifiedRiemannZeta
!	real*16 ModifiedRiemannZeta
!
!	real*16 AsymptoticSum
!	integer q, i, j, Info
!	real*16 A(LambdaUpper+1, LambdaUpper+1)
!	real*16 BX(LambdaUpper+1);
!	real*16 IPIV(N)
!
!
!	! Initialize the summation
!	AsymptoticSum = 0.0
!
!	!// Fill in x vector  with the values of T(q) for q = N-LambdaUpper to N.
!	do i = 1, LambdaUpper+1, 1
!		BX(i) = TValues(i + (N - (LambdaUpper+1)))
!	enddo
!
!	!// This is our lower value for q, q = N-LambdaUpper.
!	q = N - (LambdaUpper+1)
!
!	! Fill in our matrix, A, with the coefficients.  These coefficients are the q^(i+LambdaLower) in equation (19).
!	!  Also notice that i and j are swapped from how they normally are in nested for loops.
!	do j = 1, LambdaUpper+1, 1
!		do i = 1, LambdaUpper+1, 1
!			A(j, i) = 1.0q0 / dble(q) ** ((i-1)+dble(LambdaLower))
!		enddo
!
!		q = q + 1  ! We want all integral values for q up to N
!	enddo
!
!	! Uses LAPACK to solve our system of LambdaUpper+1 linear equations for x.
!	!LaLinearSolve(A, x, b);
!	CALL dgesv(LambdaUpper+1, 1, A, LambdaUpper+1, IPIV, BX, LambdaUpper+1, Info)
!
!	do i = 1, LambdaUpper+1, 1
!		AsymptoticSum = AsymptoticSum + BX(i) * ModifiedRiemannZeta(N-1, (i-1)+LambdaLower, ZetaUpperBound)
!	enddo
!
!	!write (*,*) 'Asymptotic part:', AsymptoticSum
!
!	AsymptoticPart = AsymptoticSum
!	return
!end
!
!
!! This is the Riemann zeta function shown in (18) with its first N terms removed.
!real*16 function ModifiedRiemannZeta(N, i, ZetaUpperBound)
!	integer N
!	real*16 i
!	integer ZetaUpperBound
!
!	real*16 Sum, Term, Zeta2
!	integer j
!
!	Sum = 0.0q0
!
!	! For the special case of i == 2, we have a p-series with a definite sum (Zeta2 below).
!	!  This series converges *VERY* slowly, so we take it as a special case.
!	if (int(i) == 2) then
!		Zeta2 = 1.6449340668482264365q0
!		do j = 1, N, 1
!			Zeta2 = Zeta2 - 1.0q0 / dble(j**i)
!		enddo
!		ModifiedRiemannZeta = Zeta2
!		return
!	endif
!
!	j = N+1
!	!do j = N+1, ZetaUpperBound, 1
!	do
!		Term = 1.0q0 / (dble(j) ** dble(i))
!		Sum = Sum + Term
!		! TODO: How small do we really need to get this?
!		if (Term < 1d-21) exit
!		j = j + 1
!	enddo
!
!	ModifiedRiemannZeta = Sum
!	return
!end

! 2F1 hypergeometric function (or Gauss's hypergeometric function) definition from
!  http://mathworld.wolfram.com/HypergeometricFunction.html
! @TODO: Do we need to worry about convergence for our purposes?
!real*16 function Hypergeometric2F1(a, b, c, z)
!	real*16 a, b, c, z
!
!	real*16 Sum
!	real*16 Term
!	real*16 s, t
!	integer nUpper
!	!integer n
!	!external Factorial, PochhammerSymbol
!	real*16 Hypergeometric
!	real*16 Factorial
!	real*16 PochhammerSymbol
!
!	common /globals/ PI, HyperB, HyperC, HyperVal
!	real*16 PI, HyperB, HyperC, HyperVal
!
!	! @TODO: Error checking all in this section
!	! @TODO: Check how stable this is.
!	if (b == HyperB - 1.0q0) then
!		if (c == HyperC - 1.0q0) then
!			!write (*,*) 'Using recursion relation'
!			s = c
!			t = b - c
!			Hypergeometric2F1 = 1.0q0 + (s + t)/s * z * HyperVal
!			HyperVal = Hypergeometric2F1
!			HyperB = b
!			HyperC = c
!			return
!		endif
!	endif
!
!	Sum = 0.0_16
!	nUpper = 200 !100
!
!!	do n = 0, nUpper, 1
!!		Term = PochhammerSymbol(a, n) * PochhammerSymbol(b, n) / PochhammerSymbol(c, n) &
!!			* z ** n / Factorial(n)
!!		if (Term <= 1d-20) then
!!			Term = 0.0_16
!!			exit
!!		endif
!!		Sum = Sum + Term
!!	enddo
!
!	!call HYGFX(a, b, c, z, Hypergeometric2F1)
!
!!	Hypergeometric2F1 = Sum
!	Hypergeometric2F1 = Hypergeometric(a,b,c,z)
!	HyperVal = Hypergeometric2F1
!	HyperB = b
!	HyperC = c
!	return
!end

real*16 function Hypergeometric2F1Sum(a, b, c, z)
	implicit none
	real*16 a, b, c, z
	!integer a, b, c
	!real*16 z
	real*16 Sum, Term, Prev, nd
	integer n, nUpper

	Sum = 1.0_16
	Prev = 0.0_16
	nUpper = 10000  ! Only for extremely singular 2F1 will it even reach close to this number.
	!nUpper = 65

	Term = 1.0_16
	do n = 1, nUpper, 1
		nd = n
		Term = Term * (a+nd-1.0_16) * (b+nd-1.0_16) / ((c+nd-1.0_16) * nd) * z
		Sum = Sum + Term
		if (abs(Sum - Prev) <= 1d-30) then
		!	write (*,*) 'Number of terms: ', n
			exit
		endif
		Prev = Sum
	enddo

	Hypergeometric2F1Sum = Sum
end function Hypergeometric2F1Sum

!  Arguments b and c are assumed to be integers, and a is 1.
recursive function Hypergeometric2F1(a, b, c, z) result (hyper_result)
	implicit none
	real*16 a, b, c, z, s, t
	real*16 Hypergeometric2F1Sum, hyper_result
	real*16 Hypergeometric2F1First

	if (c > b) then
		!write (*,*) "This must have c <= b for the recursion relation to work properly.", b, c
		hyper_result = Hypergeometric2F1Sum(a, b, c, z)
		!!!!hyper_result = Hypergeometric2F1Second(a, b, c, z)  !!! DO NOT USE !!!
		return
	end if

	!hyper_result = Hypergeometric2F1Sum(a, b, c, z)
	hyper_result = Hypergeometric2F1First(real(a,16), real(b,16), real(c,16), real(z,16))
	return
end function Hypergeometric2F1

! This is used when c <= b.
recursive function Hypergeometric2F1First(a, b, c, z) result (hyper_result)
	implicit none
	real*16 a, b, c, z, s, t, hyper_result
    real*16 temp16
    real*8 temp8

	if (c == 1) then
		hyper_result = (1-z)**(-int(b))
		return
	endif

	s = c-1
	t = b-1-s
	hyper_result = (Hypergeometric2F1First(a,b-1,c-1,z) - 1) * s/(z*(s+t))
end function Hypergeometric2F1First

! Pochhammer symbol definition from http://mathworld.wolfram.com/PochhammerSymbol.html
!  Returns x(x+1)...(x+n-1)
real*16 function PochhammerSymbol(x, n)
	real*16 x
	integer n
	real*16 p

	p = 1.0_16;

	do i = 0, n - 1, 1
		p = p * (x + i)  ! Need to convert i to a double?
	enddo

	PochhammerSymbol = p  ! @TODO: Not any real need to use the intermediate p
	return
end

max_stars_repo_name: DentonW/Ps-H-Scattering

max_issues_repo_name: DentonW/Ps-H-Scattering

max_forks_repo_name: DentonW/Ps-H-Scattering


! louvain.f90

!Algoritmo para deteccao de comunidades em rede, baseado no artigo:
!------------------------------------------------------------------
!Paper: Fast unfolding of communities in large networks
!Authors: Vincent D Blondel, Jean-Loup Guillaume, Renaud Lambiotte and Etienne Lefebvre
!doi:10.1088/1742-5468/2008/10/P10008
!------------------------------------------------------------------
!Codigo escrito por: Daniel Prado (pradofisica@gmail.com) e Derick
!Este programa eh composto de tres partes: louvain.f90; move_nodes.f90; aggregate.f90

      !variaveis comuns as 3 partes do programa
      module variaveis
      integer*8, allocatable :: C1(:), C2(:)
      integer*8, allocatable :: G1(:,:), G2(:,:)
      integer*8  m, numb_comun
      end module variaveis

      !variaveis comuns somente a louvain.f90 e suas funcoes
      module variaveis_louvain
      integer*8 n
      integer*8, allocatable :: v(:,:), C0(:)
      integer*8, allocatable :: ki(:)
      end module variaveis_louvain

      program louvain
      use variaveis
      use variaveis_louvain
      implicit none
      real*8 Q1, Q2, modularity
      integer*8 n1
      integer*8 i, j
      integer*8, allocatable ::  lar(:)
      character(len=80) arq_graph, file_out

      open(50, file = 'entra.dat')
       read(50,*) arq_graph   !Nome do arquivo com a matriz de adjacencia
       read(50,*) n           !numero de nos do grafo
      close(50)

      allocate (v(n,n), lar(n))

!--------leitura do grafo----------
!      open(100, file = "grafo5space.dat")
       open(100, file = arq_graph)

        do i=1,n
           read(100,500) (lar(j),j=1,n)
           do j=1,n
           v(i,j) = lar(j)
           end do
        end do

!---------escrita da matriz de adjacencia na tela----------
        write(*,*) 'Grafo:', arq_graph
        do i=1,n
           write(*,500) (v(i,j), j=1,n)
        end do

      close(100)
!---------verifica a simetria na matriz de adjacencia---------------------
      do i=1,n
        do j=1,n
         if(v(i,j).ne.v(j,i)) then
           write(*,*) 'sem simetria','i=', i,'j=',j
         end if
        end do
      end do
!---------------------------------------
500   format (10000I2)
!-------grau dos nos (ki) do grafo v--------
      allocate (ki(n))
      write(*,*)'grau do no i:'
      do i=1,n
            ki(i) = 0
         do j =1,n
            ki(i) = ki(i)+v(i,j)
         end do
      end do
!------- soma das arestas--------------------
      m = sum(v)/2
      write(*,*)'Numero de arestas(m) =', m

!************Louvain*************************
      allocate (C0(n), C1(n), C2(n), G1(n,n))

      do i=1,n !comunidade inicial do no i
         C0(i)=i
      end do

      numb_comun = n
      G1 = v
      Q1 = -1000.0
      Q2 = modularity()

      do while(Q1.lt.Q2)
         Q1=Q2
         call move_nodes()       !primero passo do metodo Louvain
         call aggregate()        !segundo passo do metodo Louvain

          !coloca o no do grafo original na nova comunidade
          n1 = maxval(C0)
          do i=1, n1
             where(C0.eq.i)
               C0 = C1(i)
             end where
          end do

          !calcula modularidade
          Q2 = modularity()

      end do

!****************Escrever dados de saida****************
      file_out =  'modularity_'//arq_graph
      open(51,file= file_out)
      !Escreve os dados de saida
      write(51,*) arq_graph

      write(51,*) 'Matriz de adjacencia entre as comunidades:'
      do i=1,numb_comun
        write(51,*) (G1(i,j), j=1,numb_comun)
      end do
      write(51,fmt='(/A)') 'Comunidades:'
      do i=1, numb_comun
        write(51,fmt='(A,I4,A)', advance='no') 'C', i, '['
        do j=1,n
           if(C0(j).eq.i) then
             write(51,fmt='(I5,A1)', advance='no') j,','
           end if
        end do
        write(51,fmt='(A)')']'
      end do
      write(51,fmt='(/A,I6)')'numero de comunidades = ', numb_comun
      write(51,*)'modularidade = ', Q2

      close(51)
      write(*,*)'O arquivo de saida ja esta na pasta ', file_out
      end program louvain

!********************************************************
!--------Calcula modularidade-----------------
      function modularity()
      use variaveis
      use variaveis_louvain
      real*8 modularity
      integer*8 i, j

      modularity = 0.d0
      do i=1,n
         do j=1,n
           if(C0(i).eq.C0(j)) then
             modularity = modularity+v(i,j)-ki(i)*ki(j)/(2.0*m)
           end if
         end do
      end do

      modularity = modularity/(2*m)
      return
      end function modularity

!------------sub-rotinas--------------------
      include 'move_nodes.f90'
      include 'aggregate.f90'


text:
module base_state_module
  ! init_base_state is used to initialize the base state arrays from the
  ! model file.  The actual reading of the model file is handled by the
  ! model_parser_module in Util/
  !
  ! Note: The initial base state quantities returned from this routine
  ! are only a temporary base state.  These quantities are mapped onto
  ! the full 2- or 3-d state in initscaldata.f90 and a new base state is
  ! created after initialization by averaging the density and calling
  ! enforce_HSE in initialize.f90.

  use model_parser_module
  use eos_type_module
  use eos_module
  use amrex_constants_module
  use simple_log_module
  use inlet_bc_module
  use fundamental_constants_module, only: Gconst
  use amrex_fort_module, only: amrex_spacedim
  use network, only: nspec
  use meth_params_module, only: nscal, model_file, spherical, base_cutoff_density, &
                                do_2d_planar_octant, do_planar_invsq_grav, rho_comp, &
                                rhoh_comp, spec_comp, temp_comp, grav_const, &
                                planar_invsq_mass, print_init_hse_diag, prob_lo
  use base_state_geometry_module, only: nr_fine, dr, nr, max_radial_level
  use model_util_module, only : set_species, fv, dUdy
  use probin_module, only: rho_0, p_0

  implicit none

  private

contains

  subroutine init_base_state(s0_init,p0_init,rho0,rhoh0,p0,tempbar,tempbar_init) &
       bind(C, name="init_base_state")
    ! Binds to C function ``init_base_state``

    double precision, intent(inout) :: s0_init(0:max_radial_level,0:nr_fine-1,1:nscal)
    double precision, intent(inout) :: p0_init(0:max_radial_level,0:nr_fine-1)
    double precision, intent(inout) ::    rho0(0:max_radial_level,0:nr_fine-1)
    double precision, intent(inout) ::   rhoh0(0:max_radial_level,0:nr_fine-1)
    double precision, intent(inout) ::      p0(0:max_radial_level,0:nr_fine-1)
    double precision, intent(inout) :: tempbar(0:max_radial_level,0:nr_fine-1)
    double precision, intent(inout) :: tempbar_init(0:max_radial_level,0:nr_fine-1)

    ! local variables

    double precision :: x, y, z, fheat, rhopert, U_old(2), U_new(2), h
    double precision :: k1(2), k2(2), k3(2), k4(2)
    double precision, allocatable :: pres(:), dens(:)
    double precision :: xn(nspec)
    integer :: i, j, k, n, r, n_dy
    type (eos_t) :: eos_state

    do n=0,max_radial_level

        allocate(pres(0:nr(n)-1))
        allocate(dens(0:nr(n)-1))

        U_old(1) = log(rho_0)
        U_old(2) = log(p_0)

        do r=0,nr(n)-1

            y = prob_lo(amrex_spacedim) + (dble(r) + HALF)*dr(n)

            ! do HSE using RK2

            ! our integration starts at y - h
            if (r .eq. 0) then
                h = dr(n) * HALF
            else
                h = dr(n)
            endif

            k1(:) = dUdy(y - h, U_old)
            U_new(:) = U_old(:) + h * dUdy(y - HALF*h, U_old + HALF*h * k1)

            dens(r) = exp(U_new(1))
            pres(r) = exp(U_new(2))

            U_old(:) = U_new(:)

        end do

        do r=0,nr(n)-1

            y = prob_lo(amrex_spacedim) + (dble(r) + HALF)*dr(n)

            if (r < 0) then
                eos_state % rho = rho_0
                eos_state % p = p_0
            else
                eos_state%rho = dens(r)
                eos_state%p = pres(r)
            endif

            eos_state%xn(:) = set_species(y)

            call eos(eos_input_rp, eos_state)

            s0_init(n, r, rho_comp) = eos_state % rho
            s0_init(n, r, rhoh_comp) = eos_state % rho * eos_state % h
            s0_init(n, r, spec_comp:spec_comp+nspec-1) = eos_state % rho * eos_state % xn(:)
            p0_init(n, r) = eos_state % p
            s0_init(n, r, temp_comp) = eos_state % t

        end do

       deallocate(pres)
       deallocate(dens)

       ! initialize any inlet BC parameters
       call set_inlet_bcs()

    end do ! end loop over levels

    ! copy s0_init and p0_init into rho0, rhoh0, p0, and tempbar
    rho0(:,:) = s0_init(:,:,rho_comp)
    rhoh0(:,:) = s0_init(:,:,rhoh_comp)
    tempbar(:,:) = s0_init(:,:,temp_comp)
    tempbar_init(:,:) = s0_init(:,:,temp_comp)
    p0(:,:) = p0_init(:,:)

  end subroutine init_base_state

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine init_base_state_irreg(s0_init,p0_init,rho0,rhoh0,p0,tempbar,tempbar_init, &
                                     r_cc_loc, r_edge_loc) &
       bind(C, name="init_base_state_irreg")
    ! Binds to C function ``init_base_state_irreg``

    double precision, intent(inout) :: s0_init(0:max_radial_level,0:nr_fine-1,1:nscal)
    double precision, intent(inout) :: p0_init(0:max_radial_level,0:nr_fine-1)
    double precision, intent(inout) ::    rho0(0:max_radial_level,0:nr_fine-1)
    double precision, intent(inout) ::   rhoh0(0:max_radial_level,0:nr_fine-1)
    double precision, intent(inout) ::      p0(0:max_radial_level,0:nr_fine-1)
    double precision, intent(inout) :: tempbar(0:max_radial_level,0:nr_fine-1)
    double precision, intent(inout) :: tempbar_init(0:max_radial_level,0:nr_fine-1)
    double precision, intent(in   ) ::   r_cc_loc(0:max_radial_level,0:nr_fine-1)
    double precision, intent(in   ) :: r_edge_loc(0:max_radial_level,0:nr_fine  )

  end subroutine init_base_state_irreg
end module base_state_module

max_stars_repo_name: pgrete/MAESTROeX

max_issues_repo_name: pgrete/MAESTROeX

max_forks_repo_name: pgrete/MAESTROeX


text:
!! This should be enabled only after #888
!! check if conststrings are stored correctly for special characters
!RUN: %flang -g -S -emit-llvm %s -o - | FileCheck %s

!! \0 = 0 , " = 34 = 0x22
!CHECK: call void @llvm.dbg.value(metadata [10 x i8] c"a\00d\22g     ", metadata [[CONSTR1:![0-9a-f]+]], metadata !DIExpression())
!CHECK: call void @llvm.dbg.value(metadata [10 x i8] c"h i~j     ", metadata [[CONSTR2:![0-9a-f]+]], metadata !DIExpression())
!! \n = 10 \r = 13
!CHECK: call void @llvm.dbg.value(metadata [10 x i8] c"k\0Al\0Dm     ", metadata [[CONSTR3:![0-9a-f]+]], metadata !DIExpression())
!CHECK-LABEL: distinct !DISubprogram(name: "main"
!CHECK: [[CONSTR1]] = !DILocalVariable(name: "constr1",
!CHECK: [[CONSTR2]] = !DILocalVariable(name: "constr2",
!CHECK: [[CONSTR3]] = !DILocalVariable(name: "constr3",

program main
  character(10),parameter :: constr1 = "a"//achar(0)//"d"//achar(34)//"g"
  character(10),parameter :: constr2 = "h i~j"
  character(10),parameter :: constr3 = "k"//achar(10)//"l"//achar(13)//"m"

  print *,constr1
  print *,constr2
  print *,constr3
end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
MODULE m_hubbard1Distance

   USE m_types
   USE m_constants

   IMPLICIT NONE

   CONTAINS
   SUBROUTINE hubbard1Distance(n_mmp_in,n_mmp_out,results)

      COMPLEX,             INTENT(IN)     :: n_mmp_in (-lmaxU_const:,-lmaxU_const:,:)
      COMPLEX,             INTENT(IN)     :: n_mmp_out(-lmaxU_const:,-lmaxU_const:,:)
      TYPE(t_results),     INTENT(INOUT)  :: results

      REAL    :: elementDistance(SIZE(n_mmp_in,3))
      INTEGER :: ispin,m,mp
      REAL    :: n_in,n_out

      !Calculates the distance for two density matrices (maximum distance between two elements)
      n_out = 0.0
      n_in = 0.0
      elementDistance = 0.0
      DO ispin = 1, SIZE(n_mmp_in,3)
         DO m = -lmaxU_const,lmaxU_const
            DO mp = -lmaxU_const,lmaxU_const
               !Distance between individual elements
               IF((ABS(n_mmp_out(m,mp,ispin) - n_mmp_in(m,mp,ispin))).GT.elementDistance(ispin)) THEN
                  elementDistance(ispin) = ABS(n_mmp_out(m,mp,ispin) - n_mmp_in(m,mp,ispin))
               ENDIF
               !Distance of the spin diagonal trace
               IF(m.EQ.mp.AND.ispin<3) THEN
                  n_out = n_out + REAL(n_mmp_out(m,m,ispin))
                  n_in  = n_in  + REAL(n_mmp_in (m,m,ispin))
               ENDIF
            END DO
         END DO
      ENDDO
      results%last_occdistance    = MAX(results%last_occdistance   ,ABS(n_out-n_in))
      results%last_mmpMatdistance = MAX(results%last_mmpMatdistance,MAXVAL(elementDistance))

      !IO to out file
      WRITE(oUnit,'(A)') "Hubbard 1 Distances:"
      WRITE(oUnit,'(A)') "-------------------------------"
      WRITE(oUnit,*) "  Occupation distance:     ", ABS(n_out-n_in)
      WRITE(oUnit,*) "  Density matrix distance: ", MAXVAL(elementDistance)
      DO ispin = 1, SIZE(n_mmp_in,3)
         WRITE(oUnit,'(/,A,I1,A,f14.8)') "Delta Spin ", ispin, ": ", elementDistance(ispin)
         WRITE(oUnit,'(7f12.7)') REAL(  n_mmp_out(-lmaxU_const:,-lmaxU_const:,ispin) &
                                      - n_mmp_in (-lmaxU_const:,-lmaxU_const:,ispin)  )
      ENDDO

   END SUBROUTINE hubbard1Distance
END MODULE m_hubbard1Distance

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
!========================================
! NEED TO ALLOW FOR JUST ONE HEMISPHERE !!!! - marc 27/3/03
! Marc 25/3/03
! Reads a real and imagniary parts of 1D array from a netcdf file
!========================================
      subroutine read1d_comp_nc(filein, varname, arraysize, arrayout)
      implicit none

!----------------------------------------
! For precision
!----------------------------------------
      integer :: realkind
!----------------------------------------
! Define variables in the subroutine header
!----------------------------------------
      character(*), intent(in) :: filein    ! Name of file to read from
      character(*), intent(in) :: varname   ! name of variable to collect
      integer, intent(in)      :: arraysize ! Size of array
      complex, intent(out) :: arrayout(arraysize)  ! the array which is read in
!----------------------------------------
! Define the netcdf variables
!----------------------------------------
      integer :: ncid                       ! netCDF dataset ID
     &     , status                         ! return code
     &     , var_reid                       ! ID of variable (real component)
     &     , var_imid                       ! ID of variable (imag. component)
!----------------------------------------
! Local variables
!----------------------------------------
      complex, parameter :: ic=(0.0, 1.0)   ! -i in complex number notation
      character(300)     :: varname_re      ! name for real comp. of variable
     &     , varname_im                     ! name for imag comp. of variable
      real :: array_re(arraysize)           ! real part of arrayout
     &     , array_im(arraysize)            ! imaginary part of arrayout
      integer :: length_of_varname_re  ! length of real (and imag) var. string
      integer, external :: length_of_char_nc ! funct. to calc. length of string
!----------------------------------------
! Include files
!----------------------------------------
      include 'netcdf.inc'

      realkind=kind(arrayout)

!----------------------------------------
! open the netcdf file
!----------------------------------------
      status=nf_open(filein, nf_nowrite, ncid)
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: opening ', filein
         stop
      endif
!----------------------------------------
! assume that real and imaginary part of varname are varname'_re' and
! varname'_im'.
!----------------------------------------
      varname_re=varname//'_re'
      varname_im=varname//'_im'
      length_of_varname_re=length_of_char_nc(varname_re)
!----------------------------------------
! find the variable ID
!----------------------------------------
      status=nf_inq_varid(ncid,
     &     varname_re(1:length_of_varname_re), var_reid)
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: could not find real part of ', varname
         stop
      endif
      status=nf_inq_varid(ncid,
     &     varname_im(1:length_of_varname_re), var_imid)
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: could not find imaginary part of ',
     &        varname
         stop
      endif
!----------------------------------------
! get variable
!----------------------------------------

      if (realkind.eq.4) then
        status=nf_get_var_real(ncid, var_reid, array_re)
      else if (realkind.eq.8) then
        status=nf_get_var_double(ncid, var_reid, array_re)
      else
        print*,'precision problem in read1d_comp_nc'
        stop
      endif
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: getting real part of variable'
         stop
      endif

      if (realkind.eq.4) then
        status=nf_get_var_real(ncid, var_imid, array_im)
      else if (realkind.eq.8) then
        status=nf_get_var_double(ncid, var_imid, array_im)
      else
        print*,'precision problem in read1d_comp_nc'
        stop
      endif
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: getting imaginary part of variable'
         stop
      endif
!----------------------------------------
! put real and imaginary parts together
!----------------------------------------
      arrayout=array_re+ic*array_im
!----------------------------------------
! close file
!----------------------------------------
      status=nf_close(ncid)
      if (status .ne. nf_noerr) then
         write(6, *) 'ERROR: closing ', filein
         stop
      endif
      end


program int_example
    implicit none

    integer,parameter :: nmax = 1000000
    integer :: npts, nerr
    real*4 :: data(nmax)
    real*4 :: beg, dt

    integer sac_compare

    ! Read in the data file
    call rsac1('raw.sac', data, npts, beg, dt, nmax, nerr)

    call int_trap(data, npts, dble(dt))

    ! write the seismogram with trend removed back to disk
    call wsac0('int.sac', data, data, nerr)

    if(sac_compare('int_sac.sac', data, npts-1, beg + 0.5 * dt, dt) .ne. 1) then
       write(*,*)'data does not match file'
    endif

end program int_example


text:
function func(a)
    implicit none
    integer func
    integer, intent(in) :: a
    func = a + 10
end function

max_stars_repo_name: wilsonify/fwrap

max_issues_repo_name: fwrap/fwrap

max_forks_repo_name: fwrap/fwrap


text:
c
c           user_routines_warp3d.f   Distribution version
c
c           Updated:  2/6/2015
c
c     ****************************************************************
c     *                                                              *
c     *                 subroutine user_solution_parmeters           *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 11/7/2015 rhd              *
c     *                                                              *
c     *     user routine to potentially alter solution parameters    *
c     *     and loading before the next load (time) step             *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine user_solution_parameters( now_step, iout,
     &    sol_parms, next_step_loading, convergence_history )
      implicit none
c
      integer now_step, iout ! others are user-defined types
c
c           declaration of data structure used to communicate between
c           WARP3D and the user_solution_params routine
c
c           See Section 2.10 of User Manual for more details. The
c           variables below have 1-to-1 correspondence with
c           user input commands.
c
      include 'include_usr_parm'
      type :: step_convergence_data
        logical :: step_converged
        logical :: adaptive_used
        integer :: iterations_for_convergence
        integer :: adapt_substeps
      end type
c
      type(solution_parameters) :: sol_parms
      type(step_definition) :: next_step_loading
      type(step_convergence_data), dimension(5) ::
     &                         convergence_history
c
c            locals
c
      integer i, num_patts, step_to_write
      logical debug_local
      data debug_local  / .false. /
c
      if( debug_local ) then ! values provided by WARP3D
        write(iout,*)
     &   '... called user_solution_parameters before step: ',now_step
        write(iout,9000) sol_parms%time_step, sol_parms%newmark_beta,
     &   sol_parms%maximum_iterations, sol_parms%minimum_iterations,
     &   sol_parms%nonconvergent_solutions_flag,
     &   sol_parms%adaptive_solution,
     &   sol_parms%extrapolate_solution
        write(iout,9010)
     &   sol_parms%convergence_test_1,
     &   sol_parms%convergence_test_1_tolerance,
     &   sol_parms%convergence_test_2,
     &   sol_parms%convergence_test_2_tolerance,
     &   sol_parms%convergence_test_3,
     &   sol_parms%convergence_test_3_tolerance,
     &   sol_parms%convergence_test_4,
     &   sol_parms%convergence_test_4_tolerance
        write(iout,9020)
     &   sol_parms%batch_messages,
     &   sol_parms%wall_time_limit_seconds,
     &   sol_parms%material_messages,
     &   sol_parms%trace_solution,
     &   sol_parms%bbar_stabiliation_factor,
     &   sol_parms%consistent_q_matrix,
     &   sol_parms%reset_load_reduction_factor
c
        num_patts =  next_step_loading%number_load_patt
        write(iout,9150) next_step_loading%constraints_multiplier,
     &             next_step_loading%number_load_patt
        if( num_patts .ne. 0 ) then
          do i = 1, num_patts
            write(iout,9160) next_step_loading%load_patt_nums(i),
     &      next_step_loading%load_patt_ids(i),
     &      next_step_loading%load_patt_multipliers(i)
          end do
        end if
c
c            convergence history over last 5 steps. step n is in 5th
c            row, step n-1 in row 4, etc
c
        if( now_step .gt. 1 ) write(iout,9300) now_step
        step_to_write = now_step - 1
        do i = 5, 1, -1
           if( step_to_write .lt. 1 ) exit
           write(iout,9310) step_to_write,
     &        convergence_history(i)%step_converged,
     &        convergence_history(i)%adaptive_used,
     &        convergence_history(i)%iterations_for_convergence
           if( convergence_history(i)%adaptive_used ) then
             write(iout,9320) convergence_history(i)%adapt_substeps
           end if
          step_to_write = step_to_write - 1
          end do
      end if
c
c            modify nonlinear solution parameters as desired
c            and loading/constraints for next step.
c
      if( debug_local ) then
        write(iout,*) ' '
        write(iout,*) '... leaving user_solution_parameters ...'
      end if

      return
 9000 format(/,6x,'Nonlinear solution parameters: ',
     & /,8x,'time step:                   ',e14.6,
     & /,8x,'newmark beta:                ',f8.3,
     & /,8x,'max, min global iterations:  ',2i3,
     & /,8x,'stop on nonconvergent iters: ',l2,
     & /,8x,'adaptive solution:           ',l2,
     & /,8x,'extrapolate displacements:   ',l2 )
 9010  format(
     &   8x,'convergence test 1, tol:     ',l2,1x,f10.4,
     & /,8x,'convergence test 2, tol:     ',l2,1x,f10.4,
     & /,8x,'convergence test 3, tol:     ',l2,1x,f10.4,
     & /,8x,'convergence test 4, tol:     ',l2,1x,f10.4 )
 9020 format(
     & /,8x,'batch messages:              ',l2,
     & /,8x,'wall time limits (secs):     ',f10.0,
     & /,8x,'material messages:           ',l2,
     & /,8x,'trace solution:              ',l2,
     & /,8x,'bbar stabilize factor:       ',f5.2,
     & /,8x,'consistent [Q] matrix:       ',l2,
     & /,8x,'reset load reduction factor: ',f5.2)
 9150 format(/,6x,'Next step definition:',
     & /,8x,'constraint multiplier:       ',f10.5
     & /,8x,'number of loading patterns:  ',i2 )
 9160 format(15x,'pattern #: ',i5,'   id: ',a8,3x,'multiplier: ',
     & f10.5 )
 9200 format(/ )
 9300 format(/,6x,'Convergence history. Next step is: ',i8)
 9310 format(8x,'Step:', i8,'  Converged: ',l2,
     &  '  Adaptive used: ',l2,
     &  ' global iterations: ', i4 )
 9320 format(23x,'Adaptive substeps used: ',i4)
      end
c
c
c *******************************************************************
c *                                                                 *
c *  User nodal loads routine (optional)                            *
c *                                                                 *
c *******************************************************************
c
c
      subroutine user_nodal_loads(
     &  load_name, user_file_name, pattern_values, initial_temps,
     &  nnode, step_np1, time_n, dtime,
     &  node_coords, forces_set, temps_set, use_linear_k, kout )
      implicit none
c
      character :: load_name*8, user_file_name*80
      integer  nnode, step_np1, kout
      double precision pattern_values(nnode,4), initial_temps(*),
     &                 node_coords(nnode,3),
     &                 time_n, dtime
      logical forces_set, temps_set, use_linear_k
c
c       Routine for direct specification of incremental nodal
c       forces and temperatures for a loading condition (pattern).
c
c       In the input file...
c
c           loading right_end   <--  loading name is right_end
c              nodal loads
c               user_routine (<file name>)   <-- call this routine
c                                                for loads. file name
c                                                is optional
c           loading left
c              element loads
c                  .
c                  .
c
c  Arguments:  floating points are double precision
c  ----------
c    load_name       --  character * 8 (contains string right_end in
c                        above example)
c    user_file_name  --  the character string/label (if any) provided
c                        by user in input after keyword user_routine
c    pattern_values  --  number of model nodes x 4 array. Put
c                        values of force x, y, z and
c                        temperature in cols 1-4.
c                        array was zeroed before call to this routine
c    initial_temps   --  user specified initial temperatures
c    nnode           --  number of nodes in structural model
c    step_np1        --  value of n+1 where solution is being advanced
c                        from step n to n+1 (value = 1 for first
c                        step of analysis)
c    time_n          --  model simulation time at start of step (=0 for
c                        step 1)
c    dtime           --  increment of model simulation time from
c                        n -> n+1
c    node_coords     --  number of model nodes x 3 array. contains
c                        coordinates of all structure nodes (col 1
c                        is x, 2 is y, 3 is z).
c                        These are undeformed coordinates.
c    forces_set      --  (logical) set to .true. if any nodal forces
c                        are defined. initialized to .false.
c                        before call
c    temps_set       --  (logical) set to .true. if any nodal
c                        temperatures are defined.  initialized to .false.
c                        before call
c    use_linear_k    --  return as false or true. return as true
c                        if new step loads are sufficiently non-
c                        proportional that (1) linear stiffness should
c                        be used for iteration 1 of step and (2) no
c                        extrapolation of displacements at start of
c                        step.
c    kout            --  Fortran device number for error messages,
c                        debugging output
c
c               local variables
c               ---------------
c
      double precision x
      integer node
c
c                code below is just a simple example for illustration
c
      forces_set = .true.
      temps_set = .false.
      use_linear_k = .false.
c
      do node = 1, nnode
       x = node_coords(node,1)
       if( abs(x-5.0d0) .lt. 0.001d0 ) then
         pattern_values(node,4) = 0.0d00
       end if
      end do
      return
c
      end


text:
program fibonacci_recur
    implicit none
    INTEGER:: n1
    real :: start, finish

    call cpu_time(start)
    n1 = fibonacci_recursive(38)
    call cpu_time(finish)

    print*, 'The number is: ',n1
    print*, 'Elapsed time is: ', finish - start

    contains !Separate the program from the function or sub-routine

    recursive function fibonacci_recursive(n) result(fbn) !
        integer, intent(in) :: n !Intent (in) that is the variable recieves info from outside the function
        integer :: fbn

        select case (n)

            case (:2);      fbn = 1
            case default;   fbn = fibonacci_recursive(n-1) + fibonacci_recursive(n-2)

        end select
    end function fibonacci_recursive

end program fibonacci_recur

max_stars_repo_name: arturofburgos/Comparing-Languages


text:
!==============================================================================
! Module AQUIFER_MODULE                                           (28-Jun-2017)
!
! Written by:
! 	   Dr. Randal J. Barnes
!     Department of Civil, Environmental, and Geo- Engineering
!     University of Minnesota
!     <barne003@umn.edu>
!==============================================================================

!==============================================================================
MODULE AQUIFER_MODULE
!==============================================================================
   USE CONSTANTS_MODULE
   IMPLICIT NONE

   !===========================================================================
   ! Type T_AQUIFER
   !===========================================================================
   TYPE T_AQUIFER
      REAL(8) :: Base         = 0         ! base elevation of aquifer
      REAL(8) :: Thickness    = INFINITY  ! aquifer thickness
      REAL(8) :: Porosity     = -1        ! aquifer porosity
      REAL(8) :: Conductivity = 1         ! hydraulic conductivity k
   END TYPE T_AQUIFER

   !===========================================================================
   ! INTERFACES
   !===========================================================================
   INTERFACE DischargeToVelocity
      MODULE PROCEDURE DischargeToVelocity_Aquifer
   END INTERFACE

   INTERFACE ElevationToPotential
      MODULE PROCEDURE ElevationToPotential_Aquifer
   END INTERFACE

   INTERFACE ElevationToPotentialStats
      MODULE PROCEDURE ElevationToPotentialStats_Aquifer
   END INTERFACE

   INTERFACE HeadToPotential
      MODULE PROCEDURE HeadToPotential_Aquifer
   END INTERFACE

   INTERFACE PotentialToElevation
      MODULE PROCEDURE PotentialToElevation_Aquifer
   END INTERFACE

   INTERFACE PotentialToHead
      MODULE PROCEDURE PotentialToHead_Aquifer
   END INTERFACE

   INTERFACE Reset
      MODULE PROCEDURE Reset_Aquifer
   END INTERFACE

   INTERFACE SetBase
      MODULE PROCEDURE SetBase_Aquifer
   END INTERFACE

   INTERFACE SetConductivity
      MODULE PROCEDURE SetConductivity_Aquifer
   END INTERFACE

   INTERFACE SetPorosity
      MODULE PROCEDURE SetPorosity_Aquifer
   END INTERFACE

   INTERFACE SetThickness
      MODULE PROCEDURE SetThickness_Aquifer
   END INTERFACE

CONTAINS

   !---------------------------------------------------------------------------
   ! DischargeToVelocity_Aquifer
   !---------------------------------------------------------------------------
   FUNCTION DischargeToVelocity_Aquifer( Aquifer, Head, Q ) RESULT( V )

   ! Declare the arguments.
   TYPE(T_AQUIFER),       INTENT(IN)  :: Aquifer
   REAL(8),               INTENT(IN)  :: Head
   REAL(8), DIMENSION(2), INTENT(IN)  :: Q

   REAL(8), DIMENSION(2)              :: V

   ! Compute the velocity
   IF( Head > Aquifer%Thickness ) THEN
      V = Q / ( Aquifer%Thickness * Aquifer%Porosity )
   ELSE
      V = Q / ( Head * Aquifer%Porosity )
   END IF

   END FUNCTION DischargeToVelocity_Aquifer

   !---------------------------------------------------------------------------
   ! ElevationToPotential_Aquifer
   !---------------------------------------------------------------------------
   FUNCTION ElevationToPotential_Aquifer( Aquifer, Z ) RESULT( Potential )

   ! Declare the arguments.
   TYPE(T_AQUIFER), INTENT(IN) :: Aquifer
   REAL(8),         INTENT(IN) :: Z

   REAL(8)                     :: Potential

   ! Compute the potential.
   Potential = HeadToPotential( Aquifer, Z-Aquifer%Base )

   END FUNCTION ElevationToPotential_Aquifer

   !---------------------------------------------------------------------------
   ! ElevationToPotentialStats_Aquifer
   !---------------------------------------------------------------------------
   SUBROUTINE ElevationToPotentialStats_Aquifer( Aquifer, Zev, Zstd, Pev, Pstd )

   ! Declare the arguments.
   TYPE(T_AQUIFER), INTENT(IN)  :: Aquifer
   REAL(8),         INTENT(IN)  :: Zev       ! expected value of piezometric elevation
   REAL(8),         INTENT(IN)  :: Zstd      ! standard deviation of piezometric elevation
   REAL(8),         INTENT(OUT) :: Pev       ! expected value of discharge potential
   REAL(8),         INTENT(OUT) :: Pstd      ! standard deviation of discharge potential

   ! Declare local variables.
   REAL(8) :: Head

   ! Compute the expected value and the standard deviation of the discharge
   ! potential from the expected value and the standard deviation of the piezometric
   ! elevation, using a first-order, second moment approximation.
   Head = Zev - Aquifer%Base

   IF( Head < Aquifer%Thickness ) THEN
      Pev = 0.5*Aquifer%Conductivity * ( Head**2 + Zstd**2 )
      Pstd = Aquifer%Conductivity * ABS(Head) * Zstd
   ELSE
      Pev = Aquifer%Conductivity * Aquifer%Thickness * ( Head - 0.5*Aquifer%Thickness )
      Pstd = Aquifer%Conductivity * Aquifer%Thickness * Zstd
   END IF

   END SUBROUTINE ElevationToPotentialStats_Aquifer

   !---------------------------------------------------------------------------
   ! HeadToPotential_Aquifer
   !---------------------------------------------------------------------------
   FUNCTION HeadToPotential_Aquifer( Aquifer, Head ) RESULT( Potential )

   ! Declare the arguments.
   TYPE(T_AQUIFER), INTENT(IN) :: Aquifer
   REAL(8),         INTENT(IN) :: Head

   REAL(8)                     :: Potential

   ! Compute the discharge potential from the given head.
   IF( Head < Aquifer%Thickness ) THEN
      Potential = 0.5*Aquifer%Conductivity * Head * Head
   ELSE
      Potential = Aquifer%Conductivity * Aquifer%Thickness * ( Head - 0.5*Aquifer%Thickness )
   END IF

   END FUNCTION HeadToPotential_Aquifer

   !---------------------------------------------------------------------------
   ! PotentialToElevation_Aquifer
   !---------------------------------------------------------------------------
   FUNCTION PotentialToElevation_Aquifer( Aquifer, Potential ) RESULT( Elevation )

   ! Declare the arguments.
   TYPE(T_AQUIFER), INTENT(IN) :: Aquifer
   REAL(8),         INTENT(IN) :: Potential

   REAL(8)                     :: Elevation

   ! Compute the elevation.
   Elevation = PotentialToHead( Aquifer, Potential ) + Aquifer%Base

   END FUNCTION PotentialToElevation_Aquifer

   !---------------------------------------------------------------------------
   ! PotentialToHead_Aquifer
   !
   ! Notes:
   !
   ! o   A negative discharge potential is returned as a 0 head.
   !---------------------------------------------------------------------------
   FUNCTION PotentialToHead_Aquifer( Aquifer, Potential ) RESULT( Head )

   ! Declare the arguments.
   TYPE(T_AQUIFER), INTENT(IN) :: Aquifer
   REAL(8),         INTENT(IN) :: Potential

   REAL(8)                     :: Head

   ! Declare local variables.
   REAL(8) :: Term

   ! Compute the head from the given discharge potential.
   Term = 0.5 * Aquifer%Conductivity * Aquifer%Thickness**2

   IF( Potential < EPSILON(1.0_8) ) THEN
      Head = -1
   ELSE IF( Potential < Term ) THEN
      Head = SQRT( 2*Potential / Aquifer%Conductivity )
   ELSE
      Head  = ( Potential + Term ) / ( Aquifer%Conductivity * Aquifer%Thickness )
   END IF

   END FUNCTION PotentialToHead_Aquifer

   !---------------------------------------------------------------------------
   ! Reset_Aquifer
   !---------------------------------------------------------------------------
   SUBROUTINE Reset_Aquifer( Aquifer )

   ! Declare the arguments.
   TYPE(T_AQUIFER), INTENT(INOUT) :: Aquifer

   ! Reset the parameters the defaults.
   Aquifer%Base         = 0
   Aquifer%Thickness    = INFINITY
   Aquifer%Porosity     = -1
   Aquifer%Conductivity = 1

   END SUBROUTINE Reset_Aquifer

   !---------------------------------------------------------------------------
   ! SetBase_Aquifer
   !---------------------------------------------------------------------------
   SUBROUTINE SetBase_Aquifer( Aquifer, Base )

   ! Declare the arguments.
   TYPE(T_AQUIFER), INTENT(INOUT) :: Aquifer
   REAL(8),         INTENT(IN)    :: Base

   ! Set the conductivity.
   Aquifer%Base = Base

   END SUBROUTINE SetBase_Aquifer

   !---------------------------------------------------------------------------
   ! SetConductivity_Aquifer
   !---------------------------------------------------------------------------
   SUBROUTINE SetConductivity_Aquifer( Aquifer, Conductivity )

   ! Declare the arguments.
   TYPE(T_AQUIFER), INTENT(INOUT) :: Aquifer
   REAL(8),         INTENT(IN)    :: Conductivity

   ! Set the conductivity.
   Aquifer%Conductivity = Conductivity

   END SUBROUTINE SetConductivity_Aquifer

   !---------------------------------------------------------------------------
   ! SetPorosity_Aquifer
   !---------------------------------------------------------------------------
   SUBROUTINE SetPorosity_Aquifer( Aquifer, Porosity )

   ! Declare the arguments.
   TYPE(T_AQUIFER), INTENT(INOUT) :: Aquifer
   REAL(8),         INTENT(IN)    :: Porosity

   ! Set the conductivity.
   Aquifer%Porosity = Porosity

   END SUBROUTINE SetPorosity_Aquifer

   !---------------------------------------------------------------------------
   ! SetThickness_Aquifer
   !---------------------------------------------------------------------------
   SUBROUTINE SetThickness_Aquifer( Aquifer, Thickness )

   ! Declare the arguments.
   TYPE(T_AQUIFER), INTENT(INOUT) :: Aquifer
   REAL(8),         INTENT(IN)    :: Thickness

   ! Set the conductivity.
   Aquifer%Thickness = Thickness

   END SUBROUTINE SetThickness_Aquifer

!==============================================================================
END MODULE AQUIFER_MODULE
!==============================================================================

max_stars_repo_name: RandalJBarnes/Oneka

max_issues_repo_name: RandalJBarnes/Oneka

max_forks_repo_name: RandalJBarnes/Oneka


text:
module mpi_core

  use mpi

  implicit none
  save

  private
  public :: mp_initialize
  public :: mp_stop
  public :: mp_join
  public :: main_process

  ! MPI errors
  integer :: ierr

  ! Whether MPI is enabled
  logical :: mpi_enabled = .true.

  ! Number of active processes
  integer,public :: nproc

  ! Rank of current process
  integer,public :: rank

  ! The rank of the main process
  integer,parameter,public :: rank_main = 0

contains

  logical function main_process()
    main_process = rank == rank_main
  end function main_process

  subroutine mp_initialize()
    implicit none
    call mpi_init(ierr)
    call mpi_comm_rank(mpi_comm_world, rank, ierr)
    call mpi_comm_size(mpi_comm_world, nproc, ierr)
  end subroutine mp_initialize

  subroutine mp_stop()
    implicit none
    call mpi_finalize(ierr)
  end subroutine mp_stop

  subroutine mp_join()
    implicit none
    call mpi_barrier(mpi_comm_world, ierr)
  end subroutine mp_join

end module mpi_core

max_stars_repo_name: bluescarni/hyperion

max_issues_repo_name: bluescarni/hyperion

max_forks_repo_name: bluescarni/hyperion


text:
Pine Trails Ranch is a 20+ acre equestrian facility offering boarding, lessons, breeding and horse sales/rescue. Lessons include both English and Western style ridinng including hunter, jumper, dressage, western horsemanship and trail riding. They also offer summer camps and birthday events.  They are dedicated to giving riders of all ages and abilities a fun and high quality riding experience.

To learn more about equestrian resources in or near Davis, please visit our Horses page.

20080912 11:56:26 nbsp I ride at Pine Trails Ranch and find it to be very nice there.  Every one is really nice and its easy to feel at home.  they give lessons that are about 2 hours long for good prices.  They have 5 pastures that are big and are watered.  There are 2 trails on the ranch, the perimeter trail and the half perimeter trail.  Pine Trails lives up to its name, there are lots of Pine trees.  3 arenas, one is covered and the other two are HUGE and outdoor.  The barns have really clean, nice stalls and the horses will sick their headsout for a treat.  There are two ranch dogs that are very friendly and wont get into the way.  They take the lessons slowly and their horses are trained fairly well. Although they can be a bit touchy if you are late, its a place I would love to live FOREVER! Users/lilitraw

20080922 22:07:23 nbsp I used to ride here in elementary and High school. It is a decent stable for learning to ride, but not if you want to be competitive. The emphasis is NOT on showing or competing, and students move very slowly. After a certain level, you really cant learn any more from the lessons, which are little more than pony rides. If you want to show or compete, look elsewhere. Also, it is primarily a lesson stable for children, so its not an ideal boarding stable. Users/quarabaloosa

20081018 17:22:45 nbsp I have been riding there for five years and wouldnt consider going anywhere else despite the long drive I have to get there.  I took this sport up as a 45 year old, learned to ride, and bought an amazing horse from the ranch.  The instructors are terrific and the environment is very supportive and encouraging.  The other rider/clients of the ranch are a wonderful group of people.  This is a great place for kids and adults as well (and for adults that sometimes wish they were kids) Users/kellyg

20090427 12:05:32 nbsp look quarabaloosa, you seem to be really critical. is there any place you actually WANT to ride? or do you just go around and bitch and moan about how nothing is perfect enough for you.  ptranch is a really nice stable and leasers can shw. students can slo show once in a while during the schooling shows at the ranch. the trainers typically take it slow jumping so you can perfect your form.  Users/pompomcakes

20090815 21:39:15 nbsp i love this place and the dogs and cats are so careful about not getting in the way, when i did a camp i fell off and they were so nice about asking me if i wanted to keep riding or rest for a little while (the lesson had about 5 minutes left) although some stuff is a little bit jumble its a great facility if you are a beginner Users/Churro615

20090926 15:08:44 nbsp Ummm...seriously? People like this place? Maybe we are talking about two different Pine Trails because the one I viewed had trash on the grounds, a bunch of skinny horses, some scarylooking enclosures.... why would anyone want to keep their horses there? A real letdown from the website for me, which made it look like a nice place with grass in the turnouts and happy horses not horses crowded on bonedry land with poo all over and flies swarming. Yuck no thanks. Its not like I was expecting the Ritz, but some basic good horse care would be nice. Who cares if there are lots of arenas if your horse is thin and miserable? Users/imjustjoking22

20100408 19:10:29 nbsp My horse and I have been at Pine Trails for almost two years now, and I am very happy there. The only skinny horses that are there are super old/hard keepers. All of the horses get plenty of feed and care. (Many are actually on the plump side) This is one of the best lesson/training facilities Ive had my horse at. When I first moved my horse here it looked a little shabby, but a lot has been renovated over the last year. Theres a brand new (beautiful) roundpen, and all of the fencing on the pastures have been redone. There are two big outdoor arenas with nice footing, and a smaller covered arena for the rainy season. The people who work at Pine Trails KNOW horses, and you can tell they care about all of them.  As for the lessons being only good for beginners, I disagree. Ive been riding for most of my life and I find most of my lessons to be both challenging and engaging. I also love that the lessons are catered to what the students want to be learning, and are targeted to work on individual issues that the horse and rider may be having. If youre only looking to jump gigantic fences and you dont care about safety or looking like a sack of crapprobably not the place for you. But if you want to improve your equitation and communication with your horse, I would say give the lessons here a try.  ) Users/chmorr

20110430 19:25:05 nbsp I am updating my review for PTR :)

PTR is one of the few places in Yolo County that offers affordable horse boarding, training, riding lessons and summer camps.

I have boarded by horse at Pine Trails for several years now, and I have also been teaching lessons as Michelles Assistant Instructor for about a year. As a boarder, I have always felt that my horse is in good hands while I am away. PTR feeds good quality hay and the horses are always under constant supervision. There are even affordable schooling options available for those who dont have time to ride everyday, or for when you go out of town.

2011 Summer Horse Camps are MondayFriday 8:30am5pm
$275 per child, discount for siblings and repeat campers.

Instructors teach both Western and English. This barn is perfect for beginner/intermediate riders, or those trying to build up their confidence. The horses are safe, well behaved and well trained. Although PTR has a show team, the environment around the barn is really friendly and relaxed. The perfect atmosphere to ride for fun, learn more about horses, improve your riding (and your horses training) without the snobby looks or comments.

Michelle Haseltine is the owner of PTR, as well as the Head Trainer and Instructor. Dont forget to call 24 hours in advance to cancel a lesson. No shows will be charged the fullamount for missed lessons, (with the exceptions of real emergencies)

Stop by any day of the week to check out the horses, lesson program and meet the instructors! Users/chmorr

20110430 20:51:57 nbsp I road at Pine Trails all through junior high and highschool, the only reason I stopped riding there is because i had to leave for college! I also bought my first horse from Pine Trails, Ive owned her for years now and couldnt be happier with my horse or my experience at Pine Trails. Pine Trails has grown a lot over the years that I have ridden there, there is a group for those more interested in showing, there are all sorts of wonderful clinics so you can learn about all aspects of riding and horse ownership. Prices are reasonable and a lot of time and money has been put into improving the facilities. I loved my time out at Pine Trails, I grew as both a person and a rider with the help of Michelle and all of the other trainers out at the stable. All of the horses are treated beautifully and the lesson horses are trained very well. I cant say enough good things about Pine Trails, I wish I could still be out there every day with my horse and all of the great people out there! I have been riding since I was 6 but no other stable have I ever grown as much or learned as much as I have at Pine Trails. This ranch is a great place for kids and adults, and it was great for me because I became friends with both kids my age and the adults out there. you couldnt find a more friendly group of people. Pine Trails is a wonderful place and it has changed and grown over the years so you need to see it for yourself to understand what a great place it is  Users/DRandall

20110501 09:23:05 nbsp I have rode at this ranch for a long time and have enjoyed learning how 2 ride there... plus the horses r well taken care of and the lesson teachers r very helpful Users/munch

20110501 12:29:11 nbsp I started riding at at PIne Trails Ranch as an adult beginner at 46.  Two horses later, I can say Pine Trails is a safe, supportive environment for all riders, but especially those starting out or wanting to beef up their riding skills. Good advice can also be had for working with and developing your own horse.   The ranch owners and instructors are knowledgeable, friendly and relaxed.  You can learn at your own pace, and there is a wide variety of horses available for lessons and lease at various riding levels, so you can move up when you are ready. If you are interested in owning a horse, but arent sure about the commitment, leasing is a very good way to determine if you are ready to take the plunge,  In my experience, Pine Trails priority is providing a positive horse experience for everyone at every level and they frequently go out of their way for people with a real interest in learning about horses.  Users/MeganLowry

20110710 16:17:51 nbsp My daughters and I have been riding at Pine Trails for over 2 years.  We love it! There are great teachers and horses there and it is affordable. I am an older rider (over 40) and I feel comfortable and safe there. I definitely encourage anyone interested in horseback riding to try it out! They have great summer camps and birthday parties too! Users/franmaguire

20130201 10:44:41 nbsp My daughter (just turning 5 years old) just had her first horse riding lesson here.  It was a wonderful experience for her.  Colette was her trainer, and I was highly impressed with her skills.  My daughter was very comfortable with her, which seemed due to Colette knowing exactly how to talk and interact with kids.  My daughter got to help get the horse ready for the ride, walk the horse over to the riding area, and then learn some riding techniques. Pine Trails Ranch has riding shoes and helmets for kids to borrow during the lesson (my daughter loved the riding shoes and wanted to take them home).  By the end of the lesson, my daughter was able to get the horse to stop, go, and move to the left or right.  Colette new exactly when to offer assistance and when to let my daughter try things by herself.  My daughter has been wanting to go back ever since.  It was a birthday present from her grandparents, but we unexpectedly are now trying to figure out if we can make it a more regular event.  Great job, Pine Trails Ranch! Users/Oaks

20130207 21:35:44 nbsp If you are looking to develop your riding skills and want to learn more about horsemanship, I definitely recommend this place. Ive been at Pine Trails Ranch for five months now and have learned so much during this time. Ive even learned completely new things!!, and thats not bad since I have about 25 years of experience with horses from my home country Sweden. Michelle is very good at matching horse and rider. It took her about five minutes to see what kind of horses I was used to riding, and she was right of course:) I actually cannot think of one bad thing about this place. I just love Pine Trails Ranch and am so impressed by Michelles work with both horses and riders. Users/CharlotteBhatt

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program main
    use global_variables
    use input
    use derived_parameters
    use rhs
    use DVODE_F90_M
    use calculations
    implicit none

    double precision, allocatable :: ATOL(:), Y(:)
    double precision :: RTOL, T, TOUT, RSTATS
    integer :: NEQ, ITASK, ISTATE, ISTATS, IOUT, IERROR, R, kk, STEPS
    DIMENSION :: RSTATS(22), ISTATS(31)
    TYPE(VODE_OPTS) :: OPTIONS
    OPEN(UNIT=6, FILE = 'results_raw.dat')

    call read_inputfile()
    call derived_param()

    NEQ = (N_modes*2)+1

    allocate(Y(NEQ))
    allocate(ATOL(NEQ))

    n10 = p10/(kB*Temp);
    Y(NEQ) = n10;
    do kk = 1, N_modes
        Y(2*kk-1) = N_init(kk);
		Y(2*kk-0) = v_init(kk)*Y(2*kk-1);
    end do

	! Error Tolerances
    do kk = 2, N_modes
        ATOL(2*kk-1) = (1D-20)
        ATOL(2*kk-0) = (1D-30)
    end do
    ATOL(NEQ) = (1D-10)
    RTOL = 1D-6

    T = 0
    TOUT = dt
    STEPS = t_final/dt
    ITASK = 1
    ISTATE = 1
    OPTIONS = SET_OPTS(RELERR=RTOL,ABSERR_VECTOR=ATOL)

	WRITE(6,63)T,Y
	time(1) = T
	nlines = 1
    do IOUT = 1,STEPS
        call DVODE_F90(FEX,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS)
        CALL GET_STATS(RSTATS,ISTATS)
        WRITE(6,63)T,Y
       TOUT = TOUT + dt
       time(IOUT+1) = T
	   nlines = nlines + 1
    end do

63  FORMAT(100E14.6)
call calc()
end program main

max_stars_repo_path: AAQRL-ADM/Modal model/Modal_GDM.f90

max_issues_repo_path: AAQRL-ADM/Modal model/Modal_GDM.f90

max_forks_repo_path: AAQRL-ADM/Modal model/Modal_GDM.f90


text:
      REAL*8 FUNCTION STEPI(F,FDOT,F2DOT,F3DOT,ETA)
*
*
*       Fast time-step expression.
*       --------------------------
*
      IMPLICIT REAL*8  (A-H,O-Z)
      REAL*8  F(3),FDOT(3),F2DOT(3),F3DOT(3)
*
*
*       Set square FDOT & F2DOT and sum of all absolute values.
      FD2 = FDOT(1)**2 + FDOT(2)**2 + FDOT(3)**2
      F2D2 = F2DOT(1)**2 + F2DOT(2)**2 + F2DOT(3)**2
      FI = ABS(F(1)) + ABS(F(2)) + ABS(F(3))
      FD = ABS(FDOT(1)) + ABS(FDOT(2)) + ABS(FDOT(3))
      F2D = ABS(F2DOT(1)) + ABS(F2DOT(2)) + ABS(F2DOT(3))
      F3D = ABS(F3DOT(1)) + ABS(F3DOT(2)) + ABS(F3DOT(3))
*
*       Obtain time-step by simplified relative criterion.
      STEPI = ETA*(2.0*FI*F2D + FD2)/(2.0*FD*F3D + F2D2)
      STEPI = SQRT(STEPI)
*
      RETURN
*
      END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
      subroutine addabd (n, t, y, ml, mu, pa, m0)
c This subroutine computes the matrix A in band form, adds it to pa,
c and returns the sum in pa.   The matrix A is tridiagonal, of order n,
c with nonzero elements (reading across) of  1/6, 4/6, 1/6.
c
      integer i, n, m0, ml, mu, mup1, mup2
      double precision t, y, pa, fact1, fact4, one, four, six
      dimension y(*), pa(m0,*)
      dimension n(*)
      data one/1.0d0/, four/4.0d0/, six/6.0d0/
c
c Set the pointers.
      mup1 = mu + 1
      mup2 = mu + 2
c Compute the elements of A.
      fact1 = one/six
      fact4 = four/six
c Add the matrix A to the matrix pa (banded).
      do 10 i = 1,n(1)
        pa(mu,i) = pa(mu,i) + fact1
        pa(mup1,i) = pa(mup1,i) + fact4
        pa(mup2,i) = pa(mup2,i) + fact1
   10   continue
      return
c end of subroutine addabd for the DLSODI demonstration problem.
      end

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg