!-----------------------------------------------------------------------
!
! RAdS Jan12.
! Copyright (C) 2011-2014 by Rodrigo Alvares de Souza.
! Mail: <rsouza01@gmail.com>. Web: "http://www.astro.iag.usp.br/~rsouza/".
! This program may be copied and/or distributed freely. See the
! _ terms and conditions in the files in the doc/ subdirectory.
!
!
! MODULE: solver_module
!
!> @author
!> Rodrigo Souza
!
! DESCRIPTION:
!> Runge-kutta related routines
!
!-----------------------------------------------------------------------
module solver_module
use global_constants
use config_module
use tov_module

    implicit none

    contains

    !---------------------------------------------------------------------------
    ! DESCRIPTION:
    !> The hard work
    !---------------------------------------------------------------------------

    subroutine process_calculations(parameters, error)
    implicit none
        !Parameters
        type(ConfigParameters) :: parameters
        integer :: error

        integer :: j = 0
        double precision :: t = 0

        !Variables that hold the integrations values for a single call to Runge-Kutta
        double precision :: Y(N_DIFF_EQUATIONS)

        !Variables that hold the physical relevant values
        double precision :: V(N_VARIABLES)

        t = j * parameters%diff_eq_step

        call set_initial_conditions(parameters, Y, V)

        if(.not. parameters%output_summary_only) then
            call print_header_output(parameters, t, V)
        end if

        do while ( j < parameters%max_rk_steps .and. (.not. can_stop(parameters, t, Y, V)))

            !In the first iteration, it must show the initial conditions.
            if(.not. parameters%output_summary_only) then
                call print_formated_output(parameters, t, V)
            end if

            call runge_kutta_4(parameters, t, Y, V)

            call process_specific_calculations(parameters, t, Y, V)

            j = j + 1

            t = j * parameters%diff_eq_step

        end do

        if(.not. parameters%output_summary_only) then
            call print_footer_output(parameters, t, V)
        end if

        call print_summary(parameters, t, V)

    end subroutine process_calculations

    !---------------------------------------------------------------------------
    ! DESCRIPTION:
    !> The hard work
    !---------------------------------------------------------------------------

    subroutine process_calculations_variable_step_size(parameters, error)
    implicit none
        !Parameters
        type(ConfigParameters) :: parameters
        integer :: error

        integer :: j = 0
        double precision :: t = 0

        !Variables that hold the integrations values for a single call to Runge Kutta
        double precision :: Y(N_DIFF_EQUATIONS)

        !Variables that hold the physical relevant values
        double precision :: V(N_VARIABLES)

        logical :: flag = .true.

        double precision :: calculated_h = 0.

        calculated_h = parameters%max_diff_eq_step

        t = j * calculated_h

        call set_initial_conditions(parameters, Y, V)

        if(.not. parameters%output_summary_only) then
            call print_header_output(parameters, t, V)
        end if

        do while ( flag .eqv. .true. .and. j < parameters%max_rk_steps .and. (.not. can_stop(parameters, t, Y, V)))

            !In the first iteration, it must show the initial conditions.
            if(.not. parameters%output_summary_only) then
                call print_formated_output(parameters, t, V)
            end if

            call runge_kutta_fehlberg(parameters, t, Y, V, calculated_h, flag)

            call process_specific_calculations(parameters, t, Y, V)

            j = j + 1

            t = j * calculated_h

        end do

        if(.not. parameters%output_summary_only) then
            call print_footer_output(parameters, t, V)
        end if

        call print_summary(parameters, t, V)

    end subroutine process_calculations_variable_step_size

    !---------------------------------------------------------------------------
    ! DESCRIPTION:
    !> runge_kutta_fehlberg
    !---------------------------------------------------------------------------
    subroutine runge_kutta_fehlberg(parameters, t, Y, V, calculated_h, flag)
    implicit none

        type(ConfigParameters) :: parameters
        double precision :: t, Y(N_DIFF_EQUATIONS), V(N_VARIABLES)
        logical :: flag

        double precision :: k1(N_DIFF_EQUATIONS), k2(N_DIFF_EQUATIONS), &
                            k3(N_DIFF_EQUATIONS), k4(N_DIFF_EQUATIONS), &
                            k5(N_DIFF_EQUATIONS), k6(N_DIFF_EQUATIONS)

        integer :: equationNumber
        double precision :: t_linha, calculated_h
        logical :: Rs_leq_TOL = .true. !R's < TOL
        double precision :: R(N_DIFF_EQUATIONS), delta, average_R

        t_linha = t

        !Small adjustment for t=0 in order to avoid the singularity
        if (t == 0.) then
            t_linha = t + calculated_h * 0.001
        end if

        !write (*,'(A)') '----------------------------------------'
        !1st step
        !write (*,'(A)'),'Eq.             K1'
        !write (*,'(A)') '--- --------------'
        do equationNumber = 1, N_DIFF_EQUATIONS, 1
            k1(equationNumber) = calculated_h * derivative(parameters, t_linha, Y, V, equationNumber)
        !    write (*,'(i3, f15.7)') , equationNumber, k1(equationNumber)
        end do

        !2nd step
        !write (*,'(A)'),'Eq.             K2'
        !write (*,'(A)') '--- --------------'
        do equationNumber = 1, N_DIFF_EQUATIONS, 1
            k2(equationNumber) = calculated_h * derivative(parameters, t_linha + calculated_h/4., Y + k1/4., V, equationNumber)
        !    write (*,'(i3, f15.7)') , equationNumber, k2(equationNumber)
        end do

        !3rd step
        !write (*,'(A)'),'Eq.             K3'
        !write (*,'(A)') '--- --------------'
        do equationNumber = 1, N_DIFF_EQUATIONS, 1
            k3(equationNumber) = calculated_h * derivative(parameters, &
                t_linha + 3.*calculated_h/8.,  &
                Y + 3.*k1/32.  + 9.*k2/32., &
                V, equationNumber)
        !    write (*,'(i3, f15.7)') , equationNumber, k3(equationNumber)
        end do

        !4th step
        !write (*,'(A)'),'Eq.             K4'
        !write (*,'(A)') '--- --------------'
        do equationNumber = 1, N_DIFF_EQUATIONS, 1
            k4(equationNumber) = calculated_h * derivative(parameters, &
                t_linha + (12./13.)*calculated_h,  &
                Y + (1932./2197)*k1 - (7200./2197)*k2 + (7296./2197.)*k3, &
                V, equationNumber)
        !    write (*,'(i3, f15.7)') , equationNumber, k4(equationNumber)
        end do

        !5th step
        !write (*,'(A)'),'Eq.             K4'
        !write (*,'(A)') '--- --------------'
        do equationNumber = 1, N_DIFF_EQUATIONS, 1
            k5(equationNumber) = calculated_h * derivative(parameters, &
                t_linha + calculated_h, &
                Y + (439./216.)*k1 - 8.*k2 + (3680./513)*k3 - (845./4104.)*k4, &
                V, equationNumber)
        !    write (*,'(i3, f15.7)') , equationNumber, k5(equationNumber)
        end do

        !6th step
        !write (*,'(A)'),'Eq.             K4'
        !write (*,'(A)') '--- --------------'
        do equationNumber = 1, N_DIFF_EQUATIONS, 1
            k6(equationNumber) = calculated_h * derivative(parameters, &
                t_linha + 0.5*calculated_h,  &
                Y - (8./27.)*k1 - 2.*k2 - (3544./2565.)*k3 - (1859./4104.)*k4 - (11./40.)*k5, &
                V, equationNumber)
        !    write (*,'(i3, f15.7)') , equationNumber, k6(equationNumber)
        end do

        !Calculations
        do equationNumber = 1, N_DIFF_EQUATIONS, 1
            R(equationNumber) = (1./calculated_h) * abs( &
                (1./360.)*k1(equationNumber) - &
                (128./4275.) * k2(equationNumber) - &
                (2197./75240.) * k4(equationNumber) + &
                (1./50.) * k5(equationNumber) + &
                (2./55.)*k6(equationNumber));

            !write (*,*) 'R(', equationNumber, ') = ', R(equationNumber)
            if(R(equationNumber) > parameters%interpolation_tolerance) then
                Rs_leq_TOL = .false.
            end if

            average_R = R(equationNumber) / N_DIFF_EQUATIONS

        end do

        if(Rs_leq_TOL .eqv. .true.) then
            do equationNumber = 1, N_DIFF_EQUATIONS, 1
                Y(equationNumber) = Y(equationNumber) + &
                    (25./216.)*k1(equationNumber) + &
                    (1408./2565.)*k3(equationNumber) + &
                    (2197./4104.)*k4(equationNumber) + &
                    (1./5.)*k5(equationNumber)
            end do

            calculated_h = parameters%diff_eq_step

        end if

        delta = 0.84 * (parameters%interpolation_tolerance / average_R)** (1./4.)

        if(delta < 0.1) then
            calculated_h = 0.1 * calculated_h
        elseif(delta >= 4.) then
            calculated_h = 4. * calculated_h
        else
            calculated_h = delta * calculated_h
        end if

        if(calculated_h > parameters%max_diff_eq_step) then
            calculated_h = parameters%max_diff_eq_step
        end if

        if (calculated_h < parameters%min_diff_eq_step) then
            flag = .false.
        end if

    end subroutine runge_kutta_fehlberg

    !---------------------------------------------------------------------------
    ! DESCRIPTION:
    !> runge_kutta_4
    !---------------------------------------------------------------------------
    subroutine runge_kutta_4(parameters, t, Y, V)
    implicit none

        type(ConfigParameters) :: parameters
        double precision :: t, Y(N_DIFF_EQUATIONS), V(N_VARIABLES)

        double precision :: k1(N_DIFF_EQUATIONS), k2(N_DIFF_EQUATIONS), &
                            k3(N_DIFF_EQUATIONS), k4(N_DIFF_EQUATIONS)

        integer :: equationNumber
        double precision :: t_linha, h

        h = parameters%diff_eq_step
        t_linha = t

        !Small adjustment for t=0 in order to avoid the singularity
        if (t == 0.) then
            t_linha = t + parameters%diff_eq_step * 0.001
        end if

        !write (*,'(A)') '----------------------------------------'
        !1st step
        !write (*,'(A)'),'Eq.             K1'
        !write (*,'(A)') '--- --------------'
        do equationNumber = 1, N_DIFF_EQUATIONS, 1
            k1(equationNumber) = h * derivative(parameters, t_linha, Y, V, equationNumber)
        !    write (*,'(i3, f15.7)') , equationNumber, k1(equationNumber)
        end do

        !2nd step
        !write (*,'(A)'),'Eq.             K2'
        !write (*,'(A)') '--- --------------'
        do equationNumber = 1, N_DIFF_EQUATIONS, 1
            k2(equationNumber) = h * derivative(parameters, t_linha + h/2., Y + k1/2., V, equationNumber)
        !    write (*,'(i3, f15.7)') , equationNumber, k2(equationNumber)
        end do

        !3rd step
        !write (*,'(A)'),'Eq.             K3'
        !write (*,'(A)') '--- --------------'
        do equationNumber = 1, N_DIFF_EQUATIONS, 1
            k3(equationNumber) = h * derivative(parameters, t_linha + h/2.,  Y + k2/2., V, equationNumber)
        !    write (*,'(i3, f15.7)') , equationNumber, k3(equationNumber)
        end do

        !4th step
        !write (*,'(A)'),'Eq.             K4'
        !write (*,'(A)') '--- --------------'
        do equationNumber = 1, N_DIFF_EQUATIONS, 1
            k4(equationNumber) = h * derivative(parameters, t_linha + h,  Y + k3, V, equationNumber)
        !    write (*,'(i3, f15.7)') , equationNumber, k4(equationNumber)
        end do

        !Calculations
        do equationNumber = 1, N_DIFF_EQUATIONS, 1
            Y(equationNumber) = Y(equationNumber) + ( k1(equationNumber) + 2. * k2(equationNumber) + &
                2. * k3(equationNumber) + k4(equationNumber)) / 6.;
        end do

        !write (*,'(A)') '----------------------------------------'
        !write (*,'(A)'),'         t          Massa        Pressao'
        !write (*,'(A)') '---------- -------------- --------------'

        !write (*,'(f10.5, f15.7 f15.7)') , t, Y(1), Y(2)
        !write (*,'(A)') '----------------------------------------'

    end subroutine runge_kutta_4

end module solver_module


text:
// Copyright 2019 © Lockheed Martin Corporation
// TLRB AIB PHY top-level file list
./tlrb_aib_phy_mc.sv
-F ../ip/itrx_aib_phy/rtl/itrx_aib_phy_mc.f

ext: f


text:
!#########################################################################
!# ***********************************************************************
!# <name> transform_basicRoutines </name>
!# ***********************************************************************
!#
!# <purpose>
!# This module contains basic transformation routines for lines and
!# line segments.
!#
!# </purpose>
!#########################################################################

module transform_basicRoutines

  use fsystem

  implicit none

  contains

!<function>
  !REAL(DP) FUNCTION DPO2SG(dpx,dpy,Pstart,Pend)
  real(DP) function pdist_point_StrghtLineSeg(dpx,dpy,Pstart,Pend)

  !<description>
  ! Determine minimum (Euclidean) distance from a point (dpx,dpy) to
  ! a straight line segment.
  !
  ! Return value = minimum distance
  !</description>

  !<input>

  ! dpx - x-coodinate of point
  real(DP), intent(IN) :: dpx

  ! dpy - y-coodinate of point
  real(DP), intent(IN) :: dpy

  ! Pstart - starting point of the line segment.
  !          Pstart(1) = x-coordinate
  !          Pstart(2) = y-coordinate
  real(DP), dimension(:), intent(IN) :: Pstart

  ! Pend - endpoint of the line segment.
  !        Pend(1) = x-coordinate
  !        Pend(2) = y-coordinate
  real(DP), dimension(:), intent(IN) :: Pend

  !</input>

!</function>

  ! local variables
  real(DP) :: dsx1,dsy1,dsx2,dsy2,dvx,dvy,dwx,dwy
  real(DP) :: dbx,dby,db
  real(DP) :: dvl,dpl

  dsx1 = Pstart(1)
  dsy1 = Pstart(2)
  dsx2 = Pend(1)
  dsy2 = Pend(2)

  ! Vector Pstart -> Pend and its length
  dvx = dsx2-dsx1
  dvy = dsy2-dsy1
  dvl = sqrt(dvx*dvx+dvy*dvy)

  ! Vector Pstart -> (dpx,dpy)
  dwx = dpx-dsx1
  dwy = dpy-dsy1

  ! Scalar product to calculate the length of the
  ! projected vector
  dpl = dwx*dvx+dwy*dvy

  ! Relative length of the projection:
  db  = dpl/(dvl*dvl)

  ! Relative Length <= 0
  ! => Connection between Pstart and (dpx,dpy) is shortest distance
  !
  ! Relative Length >= 1
  ! => Connection between Pend and (dpx,dpy) is shortest distance
  if (db.le.0D0) then
    pdist_point_StrghtLineSeg = sqrt((dpx-dsx1)**2+(dpy-dsy1)**2)
  else if (db.ge.1D0) then
    pdist_point_StrghtLineSeg = sqrt((dpx-dsx2)**2+(dpy-dsy2)**2)
  else
    ! Calculate the projection and the distance to that point.
    ! Remember to take the square root of the distances as they
    ! are squared...
    dbx = dsx1+dvx*db
    dby = dsy1+dvy*db
    pdist_point_StrghtLineSeg = sqrt((dpx-dbx)**2+(dpy-dby)**2)
  endif

  end function pdist_point_StrghtLineSeg

!<subroutine>
  !SUBROUTINE PPO2SG(dpx,dpy,Pstart,Pend,DPARM,DPROJ)
  subroutine pproj_point_StrghtLineSeg(dpx,dpy,Pstart,Pend,dparm,Dproj)

  !<description>
  ! Project a point (dpx,dpy) onto a line segment.
  !</description>

  !<input>

  ! dpx - x-coodinate of point
  real(DP), intent(IN) :: dpx

  ! dpy - y-coodinate of point
  real(DP), intent(IN) :: dpy

  ! Pstart - starting point of the line segment.
  !          Pstart(1) = x-coordinate
  !          Pstart(2) = y-coordinate
  real(DP), dimension(:), intent(IN) :: Pstart

  ! Pend - endpoint of the line segment.
  !        Pend(1) = x-coordinate
  !        Pend(2) = y-coordinate
  real(DP), dimension(:), intent(IN) :: Pend

  !</input>

  !<output>

  ! Dproj - Projection of (dpx,dpy) onto the line segment.
  !         Dproj(1) = x-coordinate
  !         Dproj(2) = y-coordinate
  real(DP), dimension(:), intent(OUT) :: Dproj

  ! dparm  - Parameter value of the projected point along the line
  !          segment. Values in [0,1].
  !          Dproj = Pstart + dparm*(Pend-Pstart)
  real(DP), intent(OUT) :: dparm

  !</output>

! </subroutine>

  ! local variables
  real(DP) :: dsx1,dsy1,dsx2,dsy2,dvx,dvy,dwx,dwy
  real(DP) :: dvl, dpl

  dsx1 = Pstart(1)
  dsy1 = Pstart(2)
  dsx2 = Pend(1)
  dsy2 = Pend(2)

  ! Vector Pstart -> Pend and its length
  dvx = dsx2-dsx1
  dvy = dsy2-dsy1
  dvl = sqrt(dvx*dvx+dvy*dvy)

  ! Vector Pstart -> (dpx,dpy)
  dwx = dpx-dsx1
  dwy = dpy-dsy1

  ! Scalar product to calculate the length of the
  ! projected vector
  dpl = dwx*dvx+dwy*dvy

  ! Relative length of our projected vector:
  dparm  = dpl/(dvl*dvl)

  ! Length <= 0
  ! => Point is projected onto the starting point Pstart
  !
  ! Length >= 1
  ! => Point is projected onto the endpoint Pend
  if (dparm.le.0D0) then
    dparm = 0D0
    Dproj(1)=dsx1
    Dproj(2)=dsy1
  else if (dparm.ge.1D0) then
    dparm = 1D0
    Dproj(1)=dsx2
    Dproj(2)=dsy2
  else
    ! Calculate the projection and the distance to that point.
    ! Remember to take the square root of the distances as they
    ! are squared...
    Dproj(1) = dsx1+dvx*dparm
    Dproj(2) = dsy1+dvy*dparm
  endif

  end subroutine pproj_point_StrghtLineSeg

!<function>
  !REAL(DP) FUNCTION DPO2SL(dpx,dpy,PPT1,PPT2)
  real(DP) function pdist_point_StrghtLine(dpx,dpy,Ppt1,Ppt2)

  !<description>
  ! Determine minimum (Euclidean) distance from a point (dpx,dpy) to
  ! a straight line, given by two points on the line.
  !
  ! In contrast to pdist_point_StrghtLineSeg the line has no starting/ending point!
  !
  ! Return value = minimum distance, >= 0.
  ! < 0 indicates an error: Ppt1=Ppt2
  !</description>

  !<input>

  ! dpx - x-coodinate of point
  real(DP), intent(IN) :: dpx

  ! dpy - y-coodinate of point
  real(DP), intent(IN) :: dpy

  ! Ppt1 - one point of the line segment.
  !          Ppt1(1) = x-coordinate
  !          Ppt1(2) = y-coordinate
  real(DP), dimension(:), intent(IN) :: Ppt1

  ! Ppt2   - another point of the line segment.
  !          Ppt2(1) = x-coordinate
  !          Ppt2(2) = y-coordinate
  real(DP), dimension(:), intent(IN) :: Ppt2

  !</input>

!</function>

  ! local variables
  real(DP) :: dsx1,dsy1,dsx2,dsy2,dvx,dvy,dwx,dwy
  real(DP) :: db,dbx,dby
  real(DP) :: dvl,dpl

  dsx1 = Ppt1(1)
  dsy1 = Ppt1(2)
  dsx2 = Ppt2(1)
  dsy2 = Ppt2(2)

  ! Vector Pstart -> Pend and its length DC2
  dvx = dsx2-dsx1
  dvy = dsy2-dsy1
  dvl = sqrt(dvx*dvx+dvy*dvy)

  ! Stop calculation if both points are the same.
  if (dvl.eq.0D0) then
    pdist_point_StrghtLine = -1
    return
  end if

  ! Vector Pstart -> (dpx,dpy)
   dwx = dpx-dsx1
   dwy = dpy-dsy1

   !Calculate length of projected vector
   dpl = dwx*dvx+dwy*dvy

   ! Calculate the projection
   db  = dpl/dvl
   dbx = dsx1+dvx*db
   dby = dsy1+dvy*db
   pdist_point_StrghtLine = sqrt((dpx-dby)**2+(dpy-dby)**2)

   end function pdist_point_StrghtLine

!<subroutine>
  !SUBROUTINE PPO2SL(dpx,dpy,Ppt1,Ppt2,dparm,Dproj)
  subroutine pproj_point_StrghtLine(dpx,dpy,Ppt1,Ppt2,dparm,Dproj)

  !<description>
  ! In contrast to pproj_point_StrghtLineSeg the line has no starting/ending point!
  !</description>

  !<input>

  ! dpx - x-coodinate of point
  real(DP), intent(IN) :: dpx

  ! dpy - y-coodinate of point
  real(DP), intent(IN) :: dpy

  ! Ppt1 - one point of the line segment.
  !          Ppt1(1) = x-coordinate
  !          Ppt1(2) = y-coordinate
  real(DP), dimension(:), intent(IN) :: Ppt1

  ! Ppt2   - another point of the line segment.
  !          Ppt2(1) = x-coordinate
  !          Ppt2(2) = y-coordinate
  real(DP), dimension(:), intent(IN) :: Ppt2

  !</input>

  !<output>

  ! Dproj - Projection of (dpx,dpy) onto the line segment.
  !         Dproj(1) = x-coordinate
  !         Dproj(2) = y-coordinate
  real(DP), dimension(:), intent(OUT) :: Dproj

  ! dparm  - Parameter value of the projected point along the line
  !          segment. Values in [0,1].
  !          Dproj = Pstart + dparm*(Pend-Pstart)
  real(DP), intent(OUT) :: dparm

  !</output>

! </subroutine>

  ! local variables
  real(DP) :: dsx1,dsy1,dsx2,dsy2,dvx,dvy,dwx,dwy
  real(DP) :: dvl,dpl

  dsx1 = Ppt1(1)
  dsy1 = Ppt1(2)
  dsx2 = Ppt2(1)
  dsy2 = Ppt2(2)

  ! Vector Pstart -> Pend and its length DC2
  dvx = dsx2-dsx1
  dvy = dsy2-dsy1
  dvl = sqrt(dvx*dvx+dvy*dvy)

  ! Stop calculation if both points are the same.
  if (dvl.eq.0D0) then
    dparm = 1D99
    return
  end if

  ! Vector Pstart -> (dpx,dpy)
  dwx = dpx-dsx1
  dwy = dpy-dsy1

  ! Calculate length of projected vector
  dpl = dwx*dvx+dwy*dvy

  ! Calculate the projection and the distance to that point
  dparm  = dpl/(dvl*dvl)
  Dproj(1) = dsx1+dvx*dparm
  Dproj(2) = dsy1+dvy*dparm

  end subroutine pproj_point_StrghtLine

!<function>
  !INTEGER FUNCTION PTROLN (dpx,dpy,Ppt1,Ppt2)
  integer function ptest_pointRightLeftOnLine(dpx,dpy,Ppt1,Ppt2)

  !<description>
  ! Test if a point is right-side of a line
  !
  ! This routine tests, if a given point (dpx,dpy) is right-side or
  ! left-side of a line, given by two points.
  !
  ! Return value >  0, if the point is right-side of the line Ppt1->Ppt2
  !              <  0, if the point if left-side of the line
  !              =  0, if the point is on the line or if an error
  !              occurred (Ppt1=Ppt2)
  !
  ! Remark: More specifically the return value is the squared scalar
  !         product of the normal vector of the line with the vector
  !         Ppt1->point...
  !</description>

  !<input>

  ! dpx - x-coodinate of point
  real(DP), intent(IN) :: dpx

  ! dpy - y-coodinate of point
  real(DP), intent(IN) :: dpy

  ! Ppt1 - one point of the line segment.
  !          Ppt1(1) = x-coordinate
  !          Ppt1(2) = y-coordinate
  real(DP), dimension(:), intent(IN) :: Ppt1

  ! Ppt2   - another point of the line segment.
  !          Ppt2(1) = x-coordinate
  !          Ppt2(2) = y-coordinate
  real(DP), dimension(:), intent(IN) :: Ppt2

  !</input>

  ! local variables
  real(DP) :: dnx,dny

  ! To distinguish where the point is, relative to the line, we check
  ! the vector Ppt1->point in relation to the normal vector of the line.
  !
  ! At first build the normal vector by rotating the tangential vector
  ! by -90 degrees; this points to the right side of the line.
  dnx = (Ppt2(2)-Ppt1(2))
  dny = -(Ppt2(1)-Ppt1(1))

  ! If this normal vector points into the "same" direction like the
  ! vector Ppt1->point, the point is on the right.
  ! So we build the scalar product. This is:
  !  > 0, if the point is on the right
  !  < 0, if the point is on the left
  !  = 0, if the point is on the line
  ! and thus a candidate to serve as a return value of this function!
  ptest_pointRightLeftOnLine = dnx*(dpx-Ppt1(1)) + dny*(dpy-Ppt1(2))

  end function ptest_pointRightLeftOnLine

end module transform_basicRoutines

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
!! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=45786
!!
!! Compiler (version ?) incorrectly complains (essentially) that it has
!! no == operator for the operands when in fact it should -- it appears
!! that the defined .EQ. operator is not being treated as the same
!! as == in the module.  Error occurs for other equivalent pairs of
!! operators.  Here's the compiler error:
!!
!!  print *, a == b
!!           1
!! Error: Operands of comparison operator '==' at (1) are TYPE(foo)/TYPE(foo)

module foo_type
  private
  public :: foo, operator(==)
  type :: foo
    integer :: bar
  end type
  interface operator(.eq.)
    module procedure eq_foo
  end interface
contains
  logical function eq_foo (a, b)
    type(foo), intent(in) :: a, b
    eq_foo = (a%bar == b%bar)
  end function
end module

subroutine use_it (a, b)
  use foo_type
  type(foo) :: a, b
  print *, a == b
end subroutine

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


text:
Davis Municipal Code

 Chapter 37 Tree Planting, Preservation and Protection

http://www.city.davis.ca.us/cmo/citycode/chapter.cfm?chapter37 Chapter 37 of the Davis Municipal Code.

Wiki Pages linked from here are created to discuss that section of the Code. See Davis Municipal Code for more information on link format if you want to create a page for discussion.

  37.01.0  ADMINISTRATIVE PROVISIONS
  37.01.010 Introduction and Purpose of Chapter.
  37.01.020 Definitions.
  37.01.030 Repealed by 2269.
  37.01.040 Responsibility of Parks and Community Services Director.
  37.01.050 Street Trees Planted by Subdivider.
  37.01.060 Street tree fees required when buildings constructed or remodeled.
  37.01.070 Care and planting generally.
  37.01.080 Permits.
  37.01.090 Prohibited acts generally; planting by property owners.
  37.01.100 Protection while constructing, etc., buildings.
  37.01.110 Right of inspection.
  37.01.120 Public nuisanceCertain conditions designated.
  37.01.130 SameAbatement; liability of city, etc.
  37.01.140 Appeals from decision of park and grounds superintendent; appeals from street tree commission to city council.
  37.01.150 Interference with enforcement officers prohibited.

  37.02.0  STREET TREES AND CITY TREES
  37.02.010 Street Trees Planted by Subdivider
  37.02.020 Street Trees and Street Tree Fees  Associated with Building Permits or Discretionary Projects
  37.02.030 Right and Oblication of Private Property Owners for Street Trees and City Trees
  Davis Municipal Code/37.02.040 37.02.040 Prohibited Acts
  37.02.050 Tree Modification Permit
  37.02.060 Tree Permits or Public Utilities
  37.02.070 Tree Removal Request

  37.03.0  LANDMARK TREES AND TREES OF SIGNIFICANCE LOCATED ON CERTAIN PRIVATE PROPERTY
  37.03.010 Applicability
  37.03.020 Landmark Tree  Designation Criteria
  37.03.030 Landmark Tree Designation  Benefits to Property Owners
  37.03.040 Landmark Tree Designation  Removal of Designation
  37.03.050 Trees of Significance  Identification and Classification
  Davis Municipal Code/37.03.060 37.03.060 Landmark Trees and Trees of Significance  Removal or Modification Not Associated with Building Permits or Discretionary Projects
  37.03.070 Landmark Trees and Trees of Significance  Removal or Modification Associated with Building Permits or Discretionary Projects

  37.04.0  PARKING LOT TREES
  37.04.010 Applicability.
  37.04.020 Parking Lot Shading Tree Selection and Shading Plan
  37.04.030 Parking Lot Planting/Maintenance Standards

  37.05.0  PROTECTION OF TREES DURING NEW CONSTRUCTION
  37.05.010 Protection While Trenching, Grading, Performing Construction, or Other Related Site Work

  37.06.0  VIOLATION OF CHAPTER AS PUBLIC NUISANCE
  37.06.010 Right of Inspection
  37.06.020 Public Nuisance  Certain Conditions Designated
  37.06.030 Same  Abatement Notice; Liability of City, etc.
  37.06.040 Violation  Enforcement

  37.07.0  APPENDICES
  37.07.010 Documents Available on Request

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c> \file hpsolb.f

c> This subroutine sorts out the least element of t, and puts the
c>   remaining elements of t in a heap.
c>
c> @param n On entry n is the dimension of the arrays t and iorder.

c>          On exit n is unchanged.
c>
c> @param t On entry t stores the elements to be sorted.

c>          On exit t(n) stores the least elements of t, and t(1) to t(n-1)
c>             stores the remaining elements in the form of a heap.
c>
c> @param iorder On entry iorder(i) is the index of t(i).

c>               On exit iorder(i) is still the index of t(i), but iorder may be
c>                  permuted in accordance with t.
c>
c> @param iheap On entry iheap should be set as follows:

c>
iheap .eq. 0 if t(1) to t(n) is not in the form of a heap,

c>
iheap .ne. 0 if otherwise.

c>              On exit iheap is unchanged.
      subroutine hpsolb(n, t, iorder, iheap)
      integer          iheap, n, iorder(n)
      double precision t(n)
c
c     References:
c       Algorithm 232 of CACM (J. W. J. Williams): HEAPSORT.
c
c                           *  *  *
c
c     NEOS, November 1994. (Latest revision June 1996.)
c     Optimization Technology Center.
c     Argonne National Laboratory and Northwestern University.
c     Written by
c                        Ciyou Zhu
c     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
c
c     ************

      integer          i,j,k,indxin,indxou
      double precision ddum,out

      if (iheap .eq. 0) then

c        Rearrange the elements t(1) to t(n) to form a heap.

         do 20 k = 2, n
            ddum  = t(k)
            indxin = iorder(k)

c           Add ddum to the heap.
            i = k
   10       continue
            if (i.gt.1) then
               j = i/2
               if (ddum .lt. t(j)) then
                  t(i) = t(j)
                  iorder(i) = iorder(j)
                  i = j
                  goto 10
               endif
            endif
            t(i) = ddum
            iorder(i) = indxin
   20    continue
      endif

c     Assign to 'out' the value of t(1), the least member of the heap,
c        and rearrange the remaining members to form a heap as
c        elements 1 to n-1 of t.

      if (n .gt. 1) then
         i = 1
         out = t(1)
         indxou = iorder(1)
         ddum  = t(n)
         indxin  = iorder(n)

c        Restore the heap
   30    continue
         j = i+i
         if (j .le. n-1) then
            if (t(j+1) .lt. t(j)) j = j+1
            if (t(j) .lt. ddum ) then
               t(i) = t(j)
               iorder(i) = iorder(j)
               i = j
               goto 30
            endif
         endif
         t(i) = ddum
         iorder(i) = indxin

c     Put the least member in t(n).

         t(n) = out
         iorder(n) = indxou
      endif

      return

      end

max_stars_repo_name: jonathanschilling/L-BFGS-B

max_issues_repo_name: jonathanschilling/L-BFGS-B

max_forks_repo_name: jonathanschilling/L-BFGS-B


text:
program main

  use init_data_module, only : multifab, init_data

  implicit none

  integer, parameter :: dm = 3
  type(multifab) :: data
  double precision, allocatable :: dx(:)

  data%dim = dm

  allocate(dx(data%dim))
  dx = 1.d-1

  call init_data(data,dx)

end program main

max_stars_repo_name: memmett/BoxLib

max_issues_repo_name: memmett/BoxLib

max_forks_repo_name: memmett/BoxLib


submodule (reconstruction) weno3
implicit none
! third order weno coefficients
double precision,parameter :: a200 =-0.5d0
double precision,parameter :: a201 = 1.5d0
double precision,parameter :: a210 = 0.5d0
double precision,parameter :: a211 = 0.5d0
double precision,parameter :: c20  = 1.d0/3.d0
double precision,parameter :: c21  = 2.d0/3.d0

double precision,parameter :: ss   = 1.0e-20

contains

module subroutine reconst_component_scalar_wenojs3(h,f,lo,hi,padding,bias,axis)
    implicit none
    type(rscalarfield_t),intent(inout) :: h
    type(rscalarfield_t),target,intent(in) :: f
    type(fscalarshape_t),intent(in) :: lo,hi,padding
    integer,intent(in) :: bias,axis

    double precision,dimension(:,:),pointer :: f0
    double precision,dimension(:,:),pointer :: f1
    double precision,dimension(:,:),pointer :: f2
    integer :: filo,fihi, &
               fjlo,fjhi

    double precision :: is0,is1
    double precision :: aa0,aa1,aa_sum
    double precision :: w0,w1
    double precision :: q20,q21
    integer,dimension(3) :: shift
    integer :: i,j
    if(bias.eq.bias_upwind) then
      shift = (/-1,0,1/)
    else
      shift = (/2,1,0/)
    endif

    filo = lo%d(1) + padding%d(1)
    fihi = hi%d(1) + padding%d(1)
    fjlo = lo%d(2) + padding%d(2)
    fjhi = hi%d(2) + padding%d(2)
    select case(axis)
    case(axis_xi)
      f0 => f%p(filo+shift(1):fihi+shift(1),fjlo:fjhi)
      f1 => f%p(filo+shift(2):fihi+shift(2),fjlo:fjhi)
      f2 => f%p(filo+shift(3):fihi+shift(3),fjlo:fjhi)
    case(axis_eta)
      f0 => f%p(filo:fihi,fjlo+shift(1):fjhi+shift(1))
      f1 => f%p(filo:fihi,fjlo+shift(2):fjhi+shift(2))
      f2 => f%p(filo:fihi,fjlo+shift(3):fjhi+shift(3))
    case default
      print*,"Error: reconst_component_scalar_wenojs3()"
      print*,"No such ireconst: ",axis
      stop
    end select

    do j=lo%d(2), hi%d(2)
      do i=lo%d(1), hi%d(1)
        is0 = (f0(i,j)-f1(i,j))**2
        is1 = (f1(i,j)-f2(i,j))**2

        aa0    = c20/(ss+is0)**2
        aa1    = c21/(ss+is1)**2
        aa_sum = 1.d0/(aa0+aa1)

        w0 = aa0*aa_sum
        w1 = aa1*aa_sum

        q20 = a200*f0(i,j) + a201*f1(i,j)
        q21 = a210*f1(i,j) + a211*f2(i,j)

        h%p(i,j) = w0*q20 + w1*q21
      enddo
    enddo
  end subroutine reconst_component_scalar_wenojs3

  module subroutine reconst_component_vector_wenojs3(h,f,lo,hi,padding,bias,axis)
    implicit none
    type(rvectorfield_t),intent(inout) :: h
    type(rvectorfield_t),target,intent(in) :: f
    type(fvectorshape_t),intent(in) :: lo,hi,padding
    integer,intent(in) :: bias,axis

    double precision,dimension(:,:,:),pointer :: f0
    double precision,dimension(:,:,:),pointer :: f1
    double precision,dimension(:,:,:),pointer :: f2
    integer :: fnlo,fnhi, &
               filo,fihi, &
               fjlo,fjhi

    double precision,dimension(lo%d(1):hi%d(1)) :: is0,is1
    double precision,dimension(lo%d(1):hi%d(1)) :: aa0,aa1,aa_sum
    double precision,dimension(lo%d(1):hi%d(1)) :: w0,w1
    double precision,dimension(lo%d(1):hi%d(1)) :: q20,q21
    integer :: i,j
    integer :: nlo,nhi
    integer,dimension(3) :: shift
    if(bias.eq.bias_upwind) then
      shift = (/-1,0,1/)
    else
      shift = (/2,1,0/)
    endif
    fnlo = lo%d(1) + padding%d(1)
    fnhi = hi%d(1) + padding%d(1)
    filo = lo%d(2) + padding%d(2)
    fihi = hi%d(2) + padding%d(2)
    fjlo = lo%d(3) + padding%d(3)
    fjhi = hi%d(3) + padding%d(3)
    select case(axis)
    case(axis_xi)
      f0 => f%p(fnlo:fnhi,filo+shift(1):fihi+shift(1),fjlo:fjhi)
      f1 => f%p(fnlo:fnhi,filo+shift(2):fihi+shift(2),fjlo:fjhi)
      f2 => f%p(fnlo:fnhi,filo+shift(3):fihi+shift(3),fjlo:fjhi)
    case(axis_eta)
      f0 => f%p(fnlo:fnhi,filo:fihi,fjlo+shift(1):fjhi+shift(1))
      f1 => f%p(fnlo:fnhi,filo:fihi,fjlo+shift(2):fjhi+shift(2))
      f2 => f%p(fnlo:fnhi,filo:fihi,fjlo+shift(3):fjhi+shift(3))
    case default
      print*,"Error: reconst_component_vector_wenojs3()"
      print*,"No such ireconst: ",axis
      stop
    end select

    do j=lo%d(3), hi%d(3)
      do i=lo%d(2), hi%d(2)
        is0(:) = (f0(:,i,j)-f1(:,i,j))**2
        is1(:) = (f1(:,i,j)-f2(:,i,j))**2

        aa0(:)    = c20/(ss+is0(:))**2
        aa1(:)    = c21/(ss+is1(:))**2
        aa_sum(:) = 1.d0/(aa0(:)+aa1(:))

        w0(:) = aa0(:)*aa_sum(:)
        w1(:) = aa1(:)*aa_sum(:)

        q20(:) = a200*f0(:,i,j) + a201*f1(:,i,j)
        q21(:) = a210*f1(:,i,j) + a211*f2(:,i,j)

        h%p(lo%d(1):hi%d(1),i,j) = w0*q20 + w1*q21
      enddo
    enddo
  end subroutine reconst_component_vector_wenojs3

  module subroutine reconst_component_matrix_wenojs3(h,f,lo,hi,padding,bias,axis)
    implicit none
    type(rmatrixfield_t),intent(inout) :: h
    type(rmatrixfield_t),target,intent(in) :: f
    type(fmatrixshape_t),intent(in) :: lo,hi,padding
    integer,intent(in) :: bias,axis

    double precision,dimension(:,:,:,:),pointer :: f0
    double precision,dimension(:,:,:,:),pointer :: f1
    double precision,dimension(:,:,:,:),pointer :: f2
    integer :: fmlo,fmhi, &
               fnlo,fnhi, &
               filo,fihi, &
               fjlo,fjhi

    double precision,dimension(lo%d(1):hi%d(1),lo%d(2):hi%d(2)) :: is0,is1
    double precision,dimension(lo%d(1):hi%d(1),lo%d(2):hi%d(2)) :: aa0,aa1,aa_sum
    double precision,dimension(lo%d(1):hi%d(1),lo%d(2):hi%d(2)) :: w0,w1
    double precision,dimension(lo%d(1):hi%d(1),lo%d(2):hi%d(2)) :: q20,q21
    integer :: i,j
    integer,dimension(3) :: shift
    if(bias.eq.bias_upwind) then
      shift = (/-1,0,1/)
    else
      shift = (/2,1,0/)
    endif

    fmlo = lo%d(1) + padding%d(1)
    fmhi = hi%d(1) + padding%d(1)
    fnlo = lo%d(2) + padding%d(2)
    fnhi = hi%d(2) + padding%d(2)
    filo = lo%d(3) + padding%d(3)
    fihi = hi%d(3) + padding%d(3)
    fjlo = lo%d(4) + padding%d(4)
    fjhi = hi%d(4) + padding%d(4)
    select case(axis)
    case(axis_xi)
      f0 => f%p(fmlo:fmhi,fnlo:fnhi,filo+shift(1):fihi+shift(1),fjlo:fjhi)
      f1 => f%p(fmlo:fmhi,fnlo:fnhi,filo+shift(2):fihi+shift(2),fjlo:fjhi)
      f2 => f%p(fmlo:fmhi,fnlo:fnhi,filo+shift(3):fihi+shift(3),fjlo:fjhi)
    case(axis_eta)
      f0 => f%p(fmlo:fmhi,fnlo:fnhi,filo:fihi,fjlo+shift(1):fjhi+shift(1))
      f1 => f%p(fmlo:fmhi,fnlo:fnhi,filo:fihi,fjlo+shift(2):fjhi+shift(2))
      f2 => f%p(fmlo:fmhi,fnlo:fnhi,filo:fihi,fjlo+shift(3):fjhi+shift(3))
    case default
      print*,"Error: reconst_component_matrix_wenojs3()"
      print*,"No such ireconst: ",axis
      stop
    end select

    do j = lo%d(4), hi%d(4)
      do i = lo%d(3), hi%d(3)
        is0(:,:) = (f0(:,:,i,j)-f1(:,:,i,j))**2
        is1(:,:) = (f1(:,:,i,j)-f2(:,:,i,j))**2

        aa0(:,:)    = c20/(ss+is0(:,:))**2
        aa1(:,:)    = c21/(ss+is1(:,:))**2
        aa_sum(:,:) = 1.d0/(aa0(:,:)+aa1(:,:))

        w0(:,:) = aa0(:,:)*aa_sum(:,:)
        w1(:,:) = aa1(:,:)*aa_sum(:,:)

        q20(:,:) = a200*f0(:,:,i,j) + a201*f1(:,:,i,j)
        q21(:,:) = a210*f1(:,:,i,j) + a211*f2(:,:,i,j)

        h%p(lo%d(1):hi%d(1),lo%d(2):hi%d(2),i,j) = w0*q20 + w1*q21
      enddo
    enddo
  end subroutine reconst_component_matrix_wenojs3

  module subroutine reconst_component_tensor_wenojs3(h,f,lo,hi,padding,bias,axis)
    implicit none
    type(rtensorfield_t),intent(inout) :: h
    type(rtensorfield_t),target,intent(in) :: f
    type(ftensorshape_t),intent(in) :: lo,hi,padding
    integer,intent(in) :: bias,axis

    double precision,dimension(:,:,:,:,:),pointer :: f0
    double precision,dimension(:,:,:,:,:),pointer :: f1
    double precision,dimension(:,:,:,:,:),pointer :: f2
    integer :: fllo,flhi, &
               fmlo,fmhi, &
               fnlo,fnhi, &
               filo,fihi, &
               fjlo,fjhi

    double precision,dimension(lo%d(1):hi%d(1),lo%d(2):hi%d(2),lo%d(3):hi%d(3)) :: is0,is1
    double precision,dimension(lo%d(1):hi%d(1),lo%d(2):hi%d(2),lo%d(3):hi%d(3)) :: aa0,aa1,aa_sum
    double precision,dimension(lo%d(1):hi%d(1),lo%d(2):hi%d(2),lo%d(3):hi%d(3)) :: w0,w1
    double precision,dimension(lo%d(1):hi%d(1),lo%d(2):hi%d(2),lo%d(3):hi%d(3)) :: q20,q21
    integer :: i,j
    integer,dimension(3) :: shift
    if(bias.eq.bias_upwind) then
      shift = (/-1,0,1/)
    else
      shift = (/2,1,0/)
    endif
    fllo = lo%d(1) + padding%d(1)
    flhi = hi%d(1) + padding%d(1)
    fmlo = lo%d(2) + padding%d(2)
    fmhi = hi%d(2) + padding%d(2)
    fnlo = lo%d(3) + padding%d(3)
    fnhi = hi%d(3) + padding%d(3)
    filo = lo%d(4) + padding%d(4)
    fihi = hi%d(4) + padding%d(4)
    fjlo = lo%d(5) + padding%d(5)
    fjhi = hi%d(5) + padding%d(5)
    select case(axis)
    case(axis_xi)
      f0 => f%p(fllo:flhi,fmlo:fmhi,fnlo:fnhi,filo+shift(1):fihi+shift(1),fjlo:fjhi)
      f1 => f%p(fllo:flhi,fmlo:fmhi,fnlo:fnhi,filo+shift(2):fihi+shift(2),fjlo:fjhi)
      f2 => f%p(fllo:flhi,fmlo:fmhi,fnlo:fnhi,filo+shift(3):fihi+shift(3),fjlo:fjhi)
    case(axis_eta)
      f0 => f%p(fllo:flhi,fmlo:fmhi,fnlo:fnhi,filo:fihi,fjlo+shift(1):fjhi+shift(1))
      f1 => f%p(fllo:flhi,fmlo:fmhi,fnlo:fnhi,filo:fihi,fjlo+shift(2):fjhi+shift(2))
      f2 => f%p(fllo:flhi,fmlo:fmhi,fnlo:fnhi,filo:fihi,fjlo+shift(3):fjhi+shift(3))
    case default
      print*,"Error: reconst_component_tensor_wenojs3()"
      print*,"No such ireconst: ",axis
      stop
    end select
    do j = lo%d(5), hi%d(5)
      do i = lo%d(4), hi%d(4)
        is0(:,:,:) = (f0(:,:,:,i,j)-f1(:,:,:,i,j))**2
        is1(:,:,:) = (f1(:,:,:,i,j)-f2(:,:,:,i,j))**2

        aa0(:,:,:)    = c20/(ss+is0(:,:,:))**2
        aa1(:,:,:)    = c21/(ss+is1(:,:,:))**2
        aa_sum(:,:,:) = 1.d0/(aa0(:,:,:)+aa1(:,:,:))

        w0(:,:,:) = aa0(:,:,:)*aa_sum(:,:,:)
        w1(:,:,:) = aa1(:,:,:)*aa_sum(:,:,:)

        q20(:,:,:) = a200*f0(:,:,:,i,j) + a201*f1(:,:,:,i,j)
        q21(:,:,:) = a210*f1(:,:,:,i,j) + a211*f2(:,:,:,i,j)

        h%p(lo%d(1):hi%d(1),lo%d(2):hi%d(2),lo%d(3):hi%d(3),i,j) = w0*q20 + w1*q21
      enddo
    enddo
  end subroutine reconst_component_tensor_wenojs3
end submodule


text:
c# outarea.f
      SUBROUTINE OUTciraf(ix,iy)
C--------------------------------
C
C     THIS ROUTINE OUTPUTS all variables for area coverage contour plotting
C
      use cversn
      common /chours/ nhours,ihours(24)       !  which hours are active
      common /cdaily/ idaily(12)       !  day of the month modification
      common /cantenna/ numants,iats(20),anttype(20),antname(20),
     +                  xfqs(20),xfqe(20),designfreq(20),antfile(20),
     +                  beammain(20),offazim(20),cond(20),diel(20),
     +                  array(30,91,20),aeff(30,20)
      character anttype*10,antname*70,antfile*24
      common /pantenna/ pwrkw(20),pwrdba(20)
      common /sncom/ snxx(13)
      common /cgains/ gaint(13),gainr(13)
      COMMON/TON/ADJ,ADS,GNOS,GOT,REL,SL,SLS,SPR,SU,SUS
     A ,XNOISE,ZNOISE,NF
      COMMON /DUDL_NOIS/ DU_NOIS(13),DL_NOIS(13)
      COMMON / ION / IHRE, IHRO, IHRS, LUFP, METHOD, NOISE, NPSL
      COMMON / ALPHA / IMON(12), ITRAN(2), IRCVR(2), LBMAP(2), MODE(13),
     A MODER(13), ITLAT, ITLONG, IRLAT, IRLONG, IRLATU, IRLONGU, NYEAR
      CHARACTER IMON*3, NYEAR*5, ITRAN*10, IRCVR*10, LBMAP*10, ITLAT*1,
     A ITLONG*1, IRLAT*1, IRLONG*1, IRLATU*1, IRLONGU*1, MODE*2, MODER*2
      COMMON / CONTRL / IELECT(3), KTOUT(12), MONTHS(12), SUNSP(12),
     A IANTOU, ICARD, INTEG, IRED, ISOUT, ISPROC, ISSN, ITYPE, JDASH,
     B JFREQ, JLONG, KCARD, KRUN, MAPIN, MAXNAM, MONOLD, MOREM, MORES,
     C NUMNAM, NUPROC, MAXMET, MSPEC, M100
      COMMON / FILES / LUI, LUO, LU2, LU5, LU6, LU15, LU16, LU20, LU25,
     A LU26, LU35
      COMMON /DON /ALATD, AMIN, AMIND, BTR, BTRD, DLONG, DMP, ERTR, GCD,
     1 GCDKM, PMP, PWR, TLAT, TLATD, TLONG, TLONGD, RSN, SIGTR, RLAT,
     2 RLATD,RLONG,RLONGD,BRTD,FLUX,ULAT,ULATD,ULONG,ULONGD,SSN,D90R,
     3 D50R,D10R,D90S,D50S,D10S
      COMMON/LPATH/ GCDLNG,TXRGML(45,2),DELOPT,GMIN,YMIN,LTXRGM(2)
      COMMON/FRQ/FREA(13),FREL(29),FREQ,JMODE,ITXRCP(2)
      COMMON/MUFS/EMUF(24),F1MUF(24),F2MUF(24),ESMUF(24),ALLMUF(24)
     1,FOT(24),XLUF(24),HPF(24),ANGMUF(24),MODMUF,SIGL(4),SIGU(4),DELMU
     2F(4),HPMUF(4),HTMUF(4),FVMUF(4),AFMUF(4),NHOPMF(4),YFOT(4),YHPF(4)
     3 ,YMUF(4)
      COMMON / OUTFMT / KLINE(26)
      COMMON /OUTPRT/ LINBOT(26), LINBD(14), LINTOP(15), LINTP(14),
     A GRPTYP, JOUT, LINBYP, LINES, LINMAX, LINTYP, LPAGES, NLINE
      COMMON / SON / ANGLE(13), ANGLER(13), CPROB(13), DBLOS(13),
     A DBLOSL(13), DBLOSU(13), DBU(13), DELAY(13), DBW(13), NHP(13),
     B XNYNOIS(13), PROBMP(13), RELIAB(13), SNDB(13), SNPR(13),
     C SNRLW(13), SNRUP(13), SPROB(13), VHIGH(13), RNEFF(13),MDL(13)
       CHARACTER MDL*1
      COMMON / TIME / IT, GMT, UTIME(24), GMTR, XLMT(24), ITIM, JTX
      CHARACTER KLINE*6
      COMMON / METSET / ITRUN, ITOUT, JTRUN(40), JTOUT(40)
c jw      common /CVERSN/ VERSN
c jw         character VERSN*8

      character antbear*12,slew*4
      character alfs(23)*6,alfPwr*10,label*90,path*5,daily*3
      equivalence (alfs,label)
      data SNR88,SNR91/88.,91./

      if(ix.lt.0) go to 100                  !  output header

      write(alfs(1),'(f6.2)') frel(12)            !  MUF
      if(JLONG.lt.0) then
        write(alfs(2),'(2x,i2,a2)') NHP(1),MODE(1)  !  MODE
      else
        write(alfs(2),'(2x,a2,a2)') MODE(1),MODER(1)  !  T_MODE,R_MODE
      end if
      do 50 i=2,nf
      if(dbu(i).gt.dbu(1)) dbu(1)=dbu(i)
      if(dbw(i).gt.dbw(1)) dbw(1)=dbw(i)
      if(sndb(i).gt.sndb(1)) sndb(1)=sndb(i)
      if(reliab(i).gt.reliab(1)) reliab(1)=reliab(i)
      if(sprob(i).gt.sprob(1)) sprob(1)=sprob(i)
      if(snxx(i).gt.snxx(1)) snxx(1)=snxx(i)
50    continue
c          calculate the George Lane % Power Cut
      SNR_LW=SNRLW(1)
      SNR_UP=SNRUP(1)
      SNR50=sndb(1)
      call pwrcut(SNR50,SNR_LW,SNR_UP,SNR88,SNR91,power_cut)
      write(alfs(3),'(f6.2)') xlimit6(angle(1),2)            !  ANGLE
      write(alfs(4),'(f6.2)') xlimit6(delay(1),2)            !  DELAY
      write(alfs(5),'(f6.1)') xlimit6(vhigh(1),1)            !  V HITE
      write(alfs(6),'(f6.3)') xlimit6(cprob(1),3)            !  MUF days
      write(alfs(7),'(f6.1)') xlimit6(dblos(1),1)            !  LOSS
      write(alfs(8),'(f6.1)') xlimit6(dbu(1)  ,1)            !  DBU
      write(alfs(9),'(f6.1)') xlimit6(dbw(1),1)              !  S DBW
      write(alfs(10),'(f6.1)') xlimit6(XNYNOIS(1)+RNEFF(1),1)      !  N DBW
      write(alfs(11),'(f6.1)') xlimit6(sndb(1),1)            !  SNR median
      write(alfs(12),'(f6.1)') xlimit6(snpr(1),1)            !  RPWRG
      write(alfs(13),'(f6.3)') xlimit6(reliab(1),3)          !  REL
      write(alfs(14),'(f6.3)') xlimit6(probmp(1),3)          !  MPROB
      write(alfs(15),'(f6.3)') xlimit6(sprob(1),3)           !  S PRB
      write(alfs(16),'(f6.2)') xlimit6(gaint(1),2)           !  TGAIN
      write(alfs(17),'(f6.2)') xlimit6(gainr(1),2)           !  RGAIN
      write(alfs(18),'(f6.1)') xlimit6(snxx(1),1)            !  SNRxx
ccc      write(alfs(19),'(f6.1)') xlimit6(dbw(1)+30.,1)         !  DBM
      write(alfs(19),'(f6.2)') xlimit6(DU_NOIS(1),2)         !  DU NOISE
      write(alfs(20),'(f6.2)') xlimit6(DL_NOIS(1),2)         !  DL NOISE
      write(alfs(21),'(f6.2)') xlimit6(DBLOSL(1),2)          !  SIG LW
      write(alfs(22),'(f6.2)') xlimit6(DBLOSU(1),2)          !  SIG UP
      write(alfs(23),'(f6.3)') xlimit6(power_cut,3)          !  PWRCT
      if(nf.eq.1) then
      write(LUO,'(2i3,2f10.4,23a6)') ix,iy,rlatd,rlongd,alfs
      else
      write(LUO,'(2i3,2f10.4,23a6)') ix,iy,rlatd,rlongd,alfs(1),alfs(8),
     +            alfs(9),alfs(11),alfs(13),alfs(15),alfs(18)
      end if
      return

100   write(LUO,101) VERSN
c jw 101   format('VOACAP Version',1x,a8)
101   format('VOACAPL Version',1x,a8)
      path='     '
      if(NPSL.ne.0) path='/Long'     !  using long path
      itr=1
      call TxPwr(pwrkw(itr),alfPwr,nchp)
      daily='   '
      if(idaily(1).ne.0) write(daily,'(1h.,i2.2)') idaily(1)
c**********************************************************
c         get slew angle for HFCC antennas
      islew=0
      ncha=lcount(antfile(itr),24)
      slew=antfile(itr)(ncha-3:ncha)
      if(slew(1:2).eq.'.m') then    !  -slew angle
         read(slew,'(2x,i2)',err=110) islew
         islew=-islew
      else if(slew(1:2).eq.'.p') then    !  +slew angle
         read(slew,'(2x,i2)',err=110) islew
      end if
110   if(islew.eq.0) then
         write(antbear,'(i3,3hdeg)') iabs(nint(beammain(itr)))
      else
         mainbeam=nint(beammain(itr))+islew       !  add slew back in
         if(mainbeam.lt.0.) mainbeam=mainbeam+360
         write(antbear,'(1h(,i3,1h,,i3,4h)deg)') mainbeam,islew
      end if
c**********************************************************
      if(nf.eq.1) then
      write(label,102) ITRAN,antname(itr)(1:10),alfPwr(1:nchp),
     +                 antbear,ihours(1),Frel(1),
     +                 IMON(months(1)),daily,nint(sunsp(1)),path
102   format(2a10,1h[,a10,2h] ,a,1x,a,i3.2,'ut',f7.3,'MHz ',2a3,
     +       i4,'ssn',a)
      else
      write(label,202) ITRAN,antname(itr)(1:10),alfPwr(1:nchp),
     +                 antbear,ihours(1),nf,
     +                 IMON(months(1)),daily,nint(sunsp(1)),path
202   format(2a10,1h[,a10,2h] ,a,1x,a,i3.2,'ut',i3,'Freqs ',2a3,
     +       i4,'ssn',a)
      end if
      call squeez(label,90)      !  squeeze blanks out
      nch=lcount(label,90)
ccc      write(*,'(1x,a)') label(1:nch)
      call progress_label(label)
      write(LUO,'(a)') label(1:nch)
      if(nf.eq.1) then
      write(LUO,103) -ix,-iy
103   format(2i3,'  Latitude Longitude',
     +               '   MUF  MODE ANGLE DELAY VHITE MUFda  LOSS',
     +               '   DBU  SDBW  NDBW   SNR RPWRG   REL MPROB',
     +         ' SPROB TGAIN RGAIN SNRxx    DU    DL SIGLW SIGUP PWRCT')
      else
      write(LUO,104) -ix,-iy
104   format(2i3,'  Latitude Longitude',
     +               '   MUF   DBU  SDBW   SNR   REL SPROB SNRxx')
      end if
      RETURN
      END
C---------------------------------

ext: for

max_stars_repo_name: jawatson/voacapl

max_issues_repo_name: jawatson/voacapl

max_forks_repo_name: jawatson/voacapl


text:
      SUBROUTINE DEGEN(I1,I2,ICASE)
*
*
*       Binary output for degenerate stars.
*       -----------------------------------
*
      INCLUDE 'common6.h'
      REAL*8 EB(KMAX),SEMI(KMAX),ECC(KMAX),RCM(KMAX),ECM(KMAX)
      CHARACTER*8  WHICH1
      LOGICAL FIRST,FIRST2,FIRST3
      SAVE FIRST,FIRST2,FIRST3
      DATA FIRST,FIRST2,FIRST3 /.TRUE.,.TRUE.,.TRUE./
*
*
*       Skip output for start and end of merger.
      IF (IPHASE.EQ.6.OR.IPHASE.EQ.7) GO TO 50
      IF (ICASE.EQ.7) GO TO 40
*
*       See whether KS binaries contain any degenerate stars.
      NB = 0
      DO 1 IPAIR = I1,I2
          J2 = 2*IPAIR
          J1 = J2 - 1
          IF (KSTAR(J1).GT.9.OR.KSTAR(J2).GT.9) THEN
              NB = NB + 1
          END IF
    1 CONTINUE
*
*       Open unit #4 the first time.
      IF (NB.GT.0.AND.FIRST) THEN
          OPEN (UNIT=4,STATUS='UNKNOWN',FORM='FORMATTED',FILE='DEGEN')
          FIRST = .FALSE.
*
*       Print cluster scaling parameters at start of the run.
          if(rank.eq.0)then
          WRITE (4,2)  RBAR, BODYM*ZMBAR, BODY1*ZMBAR, TSCALE,
     &                 NBIN0, NZERO
    2     FORMAT (/,6X,'MODEL:    RBAR =',F5.1,'  <M> =',F6.2,
     &                 '  M1 =',F6.1,'  TSCALE =',F6.2,
     &                 '  NB =',I4,'  N0 =',I6,//)
*
          WRITE (4,3)
    3     FORMAT (' #   TPHYS    A     E     Rp/R*      P    r',
     &            '     M1   M2    K*         NAME',/)
          end if
      END IF
*
*       Form binding energy and central distance for degenerate stars.
      JPAIR = 0
      IPRINT = 0
      DO 20 IPAIR = I1,I2
          J2 = 2*IPAIR
          J1 = J2 - 1
          IF (KSTAR(J1).LT.10.AND.KSTAR(J2).LT.10) GO TO 20
          JPAIR = JPAIR + 1
          ICM = N + IPAIR
*       Avoid division by zero for merged or synchronous ghost binary.
          IF (BODY(J1).GT.0.0) THEN
              EB(JPAIR) = BODY(J1)*BODY(J2)*H(IPAIR)/
     &                                             (BODY(J1) + BODY(J2))
              SEMI(JPAIR) = -0.5d0*BODY(ICM)/H(IPAIR)
              ECC2 = (1.d0 - R(IPAIR)/SEMI(JPAIR))**2 +
     &                       TDOT2(IPAIR)**2/(BODY(ICM)*SEMI(JPAIR))
              ECC(JPAIR) = SQRT(ECC2)
*       Set zero eccentricity after common envelope stage (still large R).
              IF (R(IPAIR).GT.2.0*SEMI(JPAIR)) THEN
                  ECC(JPAIR) = 0.d0
              END IF
              EB(JPAIR) = MAX(EB(JPAIR),-9.99999d0)
          ELSE
              EB(JPAIR) = 0.d0
              SEMI(JPAIR) = R(IPAIR)
              ECC(JPAIR) = 0.d0
          END IF
          RCM(JPAIR) = SQRT((X(1,ICM) - RDENS(1))**2 +
     &                      (X(2,ICM) - RDENS(2))**2 +
     &                      (X(3,ICM) - RDENS(3))**2)
*       Obtain binding energy (per unit mass) of c.m. motion.
          VJ2 = XDOT(1,ICM)**2 + XDOT(2,ICM)**2 +
     &                           XDOT(3,ICM)**2
*         POTJ = 0.0
*         DO 5 J = IFIRST,NTOT
*             IF (J.EQ.ICM) GO TO 5
*             RIJ2 = (X(1,ICM) - X(1,J))**2 + (X(2,ICM) - X(2,J))**2 +
*    &                                        (X(3,ICM) - X(3,J))**2
*             POTJ = POTJ + BODY(J)/SQRT(RIJ2)
*   5     CONTINUE
*         IF (TIME.EQ.TADJ) THEN
*             POTJ = -PHI(ICM)
*         ELSE
              POTJ = 0.d0
*         END IF
          ECM(JPAIR) = 0.5d0*VJ2 - POTJ
*       Check for external tidal field (NB! already included on GRAPE).
*         IF (KZ(14).GT.0) THEN
*             CALL XTRNLV(ICM,ICM)
*             ECM(JPAIR) = ECM(JPAIR) + HT/(BODY(ICM) + 1.0D-20)
*         END IF
          TPHYS = TTOT*TSTAR
          TK = SEMI(JPAIR)*SQRT(ABS(SEMI(JPAIR))/(BODY(ICM) + 1.0d-20))
          TK = DAYS*TK
          TK = MIN(TK,999999.9d0)
          RBIG = MAX(RADIUS(J1),RADIUS(J2))
          RATIO = SEMI(JPAIR)*(1.d0 - ECC(JPAIR))/RBIG
          RATIO = MIN(RATIO,99.9d0)
          IF (SEMI(JPAIR).LT.0.0.AND.RATIO.GT.5.0) GO TO 20
          SEMI(JPAIR) = SEMI(JPAIR)*RBAR*AU
          IF (IPRINT.EQ.0.AND.ICASE.EQ.0) THEN
              if(rank.eq.0) WRITE (4,*)
          END IF
          if(rank.eq.0)then
          WRITE (4,10)  ICASE, TPHYS, SEMI(JPAIR), ECC(JPAIR), RATIO,
     &                  TK, RCM(JPAIR), BODY(J1)*ZMBAR, BODY(J2)*ZMBAR,
     &                  KSTAR(J1), KSTAR(J2), KSTAR(ICM),
     &                  NAME(J1), NAME(J2)
   10     FORMAT (I2,F8.1,F8.2,F7.3,F6.1,F9.1,F6.2,2F5.1,3I4,2I6)
          end if
          IPRINT = IPRINT + 1
   20 CONTINUE
*
*       Close file at end of main output.
      IF (IPRINT.GT.0.AND.ICASE.EQ.0) THEN
          CALL FLUSH(4)
      END IF
*
*       Search for neutron stars at main output.
      IF (ICASE.EQ.0) THEN
          DO 30 J = 1,N
              IF (KSTAR(J).EQ.13) THEN
                  IF (FIRST2) THEN
                      OPEN (UNIT=33,STATUS='UNKNOWN',FORM='FORMATTED',
     &                                                 FILE='NS')
                      FIRST2 = .FALSE.
                  END IF
                  IF (J.LT.IFIRST) THEN
                      WHICH1 = ' BINARY '
                      I = KVEC(J) + N
                      VI2 = XDOT(1,I)**2 + XDOT(2,I)**2 +
     &                                     XDOT(3,I)**2
                  ELSE
                      WHICH1 = ' SINGLE '
                      VI2 = XDOT(1,J)**2 + XDOT(2,J)**2 +
     &                                     XDOT(3,J)**2
                  END IF
                  if(rank.eq.0)then
                  WRITE (33,25)  WHICH1, J, NAME(J), IFIRST, KSTAR(J),
     &                           TPHYS, SQRT(VI2)*VSTAR
   25             FORMAT (1X,A8,'NS','  J NAM I* K* TPH V ',
     &                                  2I6,I5,I4,2F7.1)
                  CALL FLUSH(33)
                  end if
              ELSE IF (KSTAR(J).GT.13) THEN
                  IF (J.LT.IFIRST) THEN
                      JCM = N + KVEC(J)
                      IF (NAME(JCM).LT.0) GO TO 30
                  END IF
                  IF (FIRST3) THEN
                      OPEN (UNIT=34,STATUS='UNKNOWN',FORM='FORMATTED',
     &                                                 FILE='BH')
                      FIRST3 = .FALSE.
                  END IF
                  VI2 = XDOT(1,J)**2 + XDOT(2,J)**2 +
     &                                 XDOT(3,J)**2
                  VI = SQRT(VI2)
                  WHICH1 = ' AIC    '
                  IF (KSTAR(J).EQ.14) WHICH1 = ' BH     '
                  IF (KSTAR(J).EQ.15) THEN
*       Ensure that errant massless remnant will escape.
                      RI = SQRT(X(1,J)**2 + X(2,J)**2 + X(3,J)**2)
                      DO 26 L = 1,3
                          X0(L,J) = 1000.0*RSCALE*X(L,J)/RI
                          X(L,J) = X0(L,J)
                          X0DOT(L,J) = SQRT(0.004*ZMASS/RSCALE)*
     &                                 XDOT(L,J)/VI
                          XDOT(L,J) = X0DOT(L,J)
   26                 CONTINUE
                  ELSE
                      VI = VI*VSTAR
                      if(rank.eq.0)then
                      WRITE (34,28) WHICH1,J,NAME(J),KSTAR(J),TPHYS,VI
   28                 FORMAT (1X,A8,'J NAM K* TPH V ',2I6,I4,2F7.1)
                      CALL FLUSH(34)
                      end if
                  END IF
              END IF
   30     CONTINUE
      END IF
*
*       Print single neutron stars at escape time.
   40 IF (ICASE.EQ.7) THEN
          J = I1
          VI2 = XDOT(1,J)**2 + XDOT(2,J)**2 + XDOT(3,J)**2
          VI = SQRT(VI2)*VSTAR
          WHICH1 = ' ESCAPE '
          if(rank.eq.0)then
          WRITE (33,25)  WHICH1, J, NAME(J), IFIRST, KSTAR(J), TPHYS, VI
          end if
      END IF
*
*       Include counter for doubly degenerate binary.
      IF (ICASE.EQ.3.OR.ICASE.EQ.4) THEN
          IPAIR = I1
          J1 = 2*IPAIR - 1
          J2 = J1 + 1
          IF (KSTAR(J1).GE.10.AND.KSTAR(J2).GE.10) THEN
              NDD = NDD + 1
              A = -0.5d0*SU*BODY(N+IPAIR)/H(IPAIR)
              if(rank.eq.0)then
              WRITE (6,48)  IPAIR, NAME(J1), NAME(J2), KSTAR(J1),
     &                      KSTAR(J2), KSTAR(N+IPAIR), R(IPAIR), A, TK
   48         FORMAT (' NEW DD    KS NM K* R A P ',I4,2I6,3I4,1P,3E10.2)
              end if
          END IF
      END IF
*
   50 RETURN
*
      END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
Describe Users/kyle2008 here.

20080706 21:59:58 nbsp  Welcome to the Wiki. Rather than just commenting on Bizarro World, how about creating a page for the other comic book store? And... just in case please read Welcome to the Wiki/Business Owner; it will help explain how you can make the wiki a positive experience for you without clashing with established wiki social norms. Users/JasonAller

20080821 03:29:00 nbsp who is this guy anyway? is he even a real person? Users/DonovanKelley

20080827 05:44:04 nbsp You are a liar, a coward and a scumbag. There will be dire consequences for your trespasses. Users/JoshLawson


text:
      SUBROUTINE ML5_0_MP_HELAS_CALLS_AMPB_1(P,NHEL,H,IC)
C
      USE ML5_0_POLYNOMIAL_CONSTANTS
      IMPLICIT NONE
C
C     CONSTANTS
C
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      INTEGER    NCOMB
      PARAMETER (NCOMB=32)

      INTEGER NBORNAMPS
      PARAMETER (NBORNAMPS=2)
      INTEGER    NLOOPS, NLOOPGROUPS, NCTAMPS
      PARAMETER (NLOOPS=11, NLOOPGROUPS=9, NCTAMPS=28)
      INTEGER    NLOOPAMPS
      PARAMETER (NLOOPAMPS=39)
      INTEGER    NWAVEFUNCS,NLOOPWAVEFUNCS
      PARAMETER (NWAVEFUNCS=10,NLOOPWAVEFUNCS=27)
      REAL*16     ZERO
      PARAMETER (ZERO=0.0E0_16)
      COMPLEX*32     IZERO
      PARAMETER (IZERO=CMPLX(0.0E0_16,0.0E0_16,KIND=16))
C     These are constants related to the split orders
      INTEGER    NSO, NSQUAREDSO, NAMPSO
      PARAMETER (NSO=0, NSQUAREDSO=0, NAMPSO=0)
C
C     ARGUMENTS
C
      REAL*16 P(0:3,NEXTERNAL)
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
      INTEGER H
C
C     LOCAL VARIABLES
C
      INTEGER I,J,K
      COMPLEX*32 COEFS(MAXLWFSIZE,0:VERTEXMAXCOEFS-1,MAXLWFSIZE)
C
C     GLOBAL VARIABLES
C
      INCLUDE 'mp_coupl_same_name.inc'

      INTEGER GOODHEL(NCOMB)
      LOGICAL GOODAMP(NSQUAREDSO,NLOOPGROUPS)
      COMMON/ML5_0_FILTERS/GOODAMP,GOODHEL

      INTEGER SQSO_TARGET
      COMMON/ML5_0_SOCHOICE/SQSO_TARGET

      LOGICAL UVCT_REQ_SO_DONE,MP_UVCT_REQ_SO_DONE,CT_REQ_SO_DONE
     $ ,MP_CT_REQ_SO_DONE,LOOP_REQ_SO_DONE,MP_LOOP_REQ_SO_DONE
     $ ,CTCALL_REQ_SO_DONE,FILTER_SO
      COMMON/ML5_0_SO_REQS/UVCT_REQ_SO_DONE,MP_UVCT_REQ_SO_DONE
     $ ,CT_REQ_SO_DONE,MP_CT_REQ_SO_DONE,LOOP_REQ_SO_DONE
     $ ,MP_LOOP_REQ_SO_DONE,CTCALL_REQ_SO_DONE,FILTER_SO

      COMPLEX*32 AMP(NBORNAMPS)
      COMMON/ML5_0_MP_AMPS/AMP
      COMPLEX*32 W(20,NWAVEFUNCS)
      COMMON/ML5_0_MP_W/W

      COMPLEX*32 WL(MAXLWFSIZE,0:LOOPMAXCOEFS-1,MAXLWFSIZE
     $ ,0:NLOOPWAVEFUNCS)
      COMPLEX*32 PL(0:3,0:NLOOPWAVEFUNCS)
      COMMON/ML5_0_MP_WL/WL,PL

      COMPLEX*32 AMPL(3,NCTAMPS)
      COMMON/ML5_0_MP_AMPL/AMPL

C
C     ----------
C     BEGIN CODE
C     ----------

C     The target squared split order contribution is already reached
C      if true.
      IF (FILTER_SO.AND.MP_CT_REQ_SO_DONE) THEN
        GOTO 1001
      ENDIF

      CALL MP_IXXXXX(P(0,1),ZERO,NHEL(1),+1*IC(1),W(1,1))
      CALL MP_OXXXXX(P(0,2),ZERO,NHEL(2),-1*IC(2),W(1,2))
      CALL MP_OXXXXX(P(0,3),ZERO,NHEL(3),+1*IC(3),W(1,3))
      CALL MP_IXXXXX(P(0,4),ZERO,NHEL(4),-1*IC(4),W(1,4))
      CALL MP_VXXXXX(P(0,5),ZERO,NHEL(5),+1*IC(5),W(1,5))
      CALL MP_FFV1_2(W(1,1),W(1,5),GC_5,ZERO,ZERO,W(1,6))
      CALL MP_FFV2_3(W(1,4),W(1,3),GC_47,MDL_MW,MDL_WW,W(1,7))
C     Amplitude(s) for born diagram with ID 1
      CALL MP_FFV2_0(W(1,6),W(1,2),W(1,7),GC_47,AMP(1))
      CALL MP_FFV1_1(W(1,2),W(1,5),GC_5,ZERO,ZERO,W(1,8))
C     Amplitude(s) for born diagram with ID 2
      CALL MP_FFV2_0(W(1,1),W(1,8),W(1,7),GC_47,AMP(2))
      CALL MP_FFV2_1(W(1,2),W(1,7),GC_47,ZERO,ZERO,W(1,9))
C     Counter-term amplitude(s) for loop diagram number 3
      CALL MP_R2_QQ_1_0(W(1,6),W(1,9),R2_QQQ,AMPL(1,1))
C     Counter-term amplitude(s) for loop diagram number 4
      CALL MP_FFV2_0(W(1,6),W(1,2),W(1,7),R2_SXCW,AMPL(1,2))
C     Counter-term amplitude(s) for loop diagram number 5
      CALL MP_FFV2_0(W(1,1),W(1,8),W(1,7),R2_SXCW,AMPL(1,3))
C     Counter-term amplitude(s) for loop diagram number 7
      CALL MP_FFV1_0(W(1,1),W(1,9),W(1,5),UV_GQQQ_1EPS,AMPL(2,4))
      CALL MP_FFV1_0(W(1,1),W(1,9),W(1,5),UV_GQQQ_1EPS,AMPL(2,5))
      CALL MP_FFV1_0(W(1,1),W(1,9),W(1,5),UV_GQQQ_1EPS,AMPL(2,6))
      CALL MP_FFV1_0(W(1,1),W(1,9),W(1,5),UV_GQQQ_1EPS,AMPL(2,7))
      CALL MP_FFV1_0(W(1,1),W(1,9),W(1,5),UV_GQQB,AMPL(1,8))
      CALL MP_FFV1_0(W(1,1),W(1,9),W(1,5),UV_GQQQ_1EPS,AMPL(2,9))
      CALL MP_FFV1_0(W(1,1),W(1,9),W(1,5),UV_GQQT,AMPL(1,10))
      CALL MP_FFV1_0(W(1,1),W(1,9),W(1,5),UV_GQQQ_1EPS,AMPL(2,11))
      CALL MP_FFV1_0(W(1,1),W(1,9),W(1,5),UV_GQQG_1EPS,AMPL(2,12))
      CALL MP_FFV1_0(W(1,1),W(1,9),W(1,5),R2_GQQ,AMPL(1,13))
      CALL MP_FFV2_2(W(1,1),W(1,7),GC_47,ZERO,ZERO,W(1,10))
C     Counter-term amplitude(s) for loop diagram number 11
      CALL MP_R2_QQ_1_0(W(1,10),W(1,8),R2_QQQ,AMPL(1,14))
C     Counter-term amplitude(s) for loop diagram number 12
      CALL MP_FFV1_0(W(1,10),W(1,2),W(1,5),UV_GQQQ_1EPS,AMPL(2,15))
      CALL MP_FFV1_0(W(1,10),W(1,2),W(1,5),UV_GQQQ_1EPS,AMPL(2,16))
      CALL MP_FFV1_0(W(1,10),W(1,2),W(1,5),UV_GQQQ_1EPS,AMPL(2,17))
      CALL MP_FFV1_0(W(1,10),W(1,2),W(1,5),UV_GQQQ_1EPS,AMPL(2,18))
      CALL MP_FFV1_0(W(1,10),W(1,2),W(1,5),UV_GQQB,AMPL(1,19))
      CALL MP_FFV1_0(W(1,10),W(1,2),W(1,5),UV_GQQQ_1EPS,AMPL(2,20))
      CALL MP_FFV1_0(W(1,10),W(1,2),W(1,5),UV_GQQT,AMPL(1,21))
      CALL MP_FFV1_0(W(1,10),W(1,2),W(1,5),UV_GQQQ_1EPS,AMPL(2,22))
      CALL MP_FFV1_0(W(1,10),W(1,2),W(1,5),UV_GQQG_1EPS,AMPL(2,23))
      CALL MP_FFV1_0(W(1,10),W(1,2),W(1,5),R2_GQQ,AMPL(1,24))

      GOTO 1001
 2000 CONTINUE
      MP_CT_REQ_SO_DONE=.TRUE.
 1001 CONTINUE
      END

max_stars_repo_name: madnklo/madnklo

max_issues_repo_name: madnklo/madnklo

max_forks_repo_name: madnklo/madnklo


text:
c=======================================================================
c
c     subroutine EIGENVD
c
c     Singular Value Decomposition   A = U * Diag(S) * VT
c
c-----------------------------------------------------------------------
      SUBROUTINE eigenvd (n, Q, dwork, D, U, info)
c-----------------------------------------------------------------------
c
c     INPUT
c       n : integer, size of matrix Q.
c       Q : (n,m)-double matrix.
c
c     WORKSPACE
c       dwork: double workspace of size n*n + 6*n + 4*n
c
c     OUTPUT
c       D: double (n)-vector of eigenvalues.
c       U: double (n,n)-matrix U.
c       info: integer,
c
c     CALL
c       YM, DGEEV (LAPACK)
c------------------------------------------------------------------------
c
      implicit none
c
c     i/o arguments
      INTEGER n, info
      DOUBLE PRECISION  Q(n,*), dwork(*), U(*), D(*)
c
c     local variables
      INTEGER lwork, ldvl, pdatmp, pdwork, pdwi, pdvl
c
      ldvl  = 1
      lwork = 4*n + 4*n
c
c     double workspace
      pdatmp = 1
c     pdatemp needs n*n
      pdwi = pdatmp + n*n
c     pdwi needs n
      pdvl = pdwi + n
c     pdwl needs n
      pdwork = pdvl + n
c     dwork needds n*n + 6*n + 4*n
c
c     matrix copy
      CALL YM(n, n, Q, dwork(pdatmp))
c
c     DGEEV computes for an N-by-N real nonsymmetric matrix
c     the eigenvalues and eigenvectors
      CALL DGEEV('N', 'V', n, dwork(pdatmp), n, D, dwork(pdwi),
     &           dwork(pdvl), ldvl, U, n, dwork(pdwork), lwork, info)
      RETURN
      END

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
C
C	$Id: gelms.f,v 1.4 2008-07-27 00:20:57 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
        SUBROUTINE GELMS(IOS,STATUS)
C
C  This routine invokes the processors for the valid
C  graphical elements.
C
C  Supported graphical elements:
C       POLYLINE
C       POLYMARKER
C       TEXT
C       POLYGON
C       CELL ARRAY
C       GENERALIZED DRAWING PRIMITIVE
C
      include 'trinst.h'
      include 'trcode.h'
C
      INTEGER IOS, STATUS
C
      IF (OPID .EQ. GPELPL) THEN
C
C  POLYLINE
C
        CALL GPUTPT(11)
C
      ELSE IF (OPID .EQ. GPELPM) THEN
C
C  POLYMARKER
C
        CALL GPUTPT(12)
C
      ELSE IF (OPID .EQ. GPELPG) THEN
C
C  POLYGON
C
        CALL GPUTPT(14)
C
      ELSE IF (OPID .EQ. GPELTX) THEN
C
C  TEXT
C
        CALL GTXDRW(IOS,STATUS)
C
      ELSE IF (OPID .EQ. GPELCA) THEN
C
C  CELL ARRAY
C
        CALL GCELDR(IOS,STATUS)
C
      ELSE IF (OPID.EQ.GPELGP) THEN
C
C  GENERALIZED DRAWING PRIMITIVE (Currently not supported).
C
        CALL GSKPOP(8,LEN,IOS,STATUS)
      END IF
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
! { dg-do compile }
! Tests the fix for PR32464, where the use associated procedure would
! mess up the check for "grandparent" host association.
!
! Contributed by Harald Anlauf <anlauf@gmx.de>
!

module gfcbug64_mod1
  implicit none

  public :: inverse

  interface inverse
     module procedure copy
  end interface

contains

  function copy (d) result (y)
    real, intent(in) :: d(:)
    real             :: y(size (d))     ! <- this version kills gfortran
!    real, intent(in) :: d
!    real             :: y
    y = d
  end function copy

end module gfcbug64_mod1

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

module gfcbug64_mod2
  implicit none
contains

  subroutine foo (x_o)
    real, intent(in) :: x_o(:)

    integer          :: s(size (x_o))           ! <- this line kills gfortran

  contains

    subroutine bar ()
      use gfcbug64_mod1, only: inverse          ! <- this line kills gfortran
    end subroutine bar

  end subroutine foo
end module gfcbug64_mod2
! { dg-final { cleanup-modules "gfcbug64_mod1 gfcbug64_mod2" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
module greatcircle
  !!Author   : Masashi Ogiso (masashi.ogiso@gmail.com)
  !!Copyright: (c) Masashi Ogiso 2020
  !!License  : MIT License https://opensource.org/licenses/MIT
  private
  public :: greatcircle_dist, latctog, latgtoc

  contains

  subroutine greatcircle_dist(lat0, lon0, lat1, lon1, delta_out, distance, azimuth, backazimuth)
    use nrtype, only : fp
    use constants, only : pi, r_earth, r2r1, deg2rad, rad2deg
    implicit none

    !!Calculate distance, azimuth, backazimuth using Vincenty (1975, Survey review)
    !!Input: (lat0, lon0) (point A) , (lat1, lon1) (point B)
    !!output: delta(rad), dist(km), az(rad, point A -> point B), baz(rad, point B -> point A)

    real(kind = fp), intent(in) :: lat0, lon0, lat1, lon1
    real(kind = fp), intent(out), optional :: delta_out, distance, azimuth, backazimuth

    real(kind = fp) :: c, b, delta, lat_0, lat_1, sin_lat0, cos_lat0, sin_lat1, cos_lat1, lon_diff, &
    &                  cos_ldiff, sin_ldiff

    !!geographical latitude to geocentral latitude
    call latgtoc(lat0, lat_0)
    call latgtoc(lat1, lat_1)

    sin_lat0 = sin(lat_0)
    cos_lat0 = cos(lat_0)
    sin_lat1 = sin(lat_1)
    cos_lat1 = cos(lat_1)
    lon_diff = (lon1 - lon0) * deg2rad
    cos_ldiff = cos(lon_diff)
    sin_ldiff = sin(lon_diff)
    delta = atan2(sqrt(cos_lat1 * sin_ldiff * cos_lat1 * sin_ldiff &
    &                  + (cos_lat0 * sin_lat1 - sin_lat0 * cos_lat1 * cos_ldiff) &
    &                  * (cos_lat0 * sin_lat1 - sin_lat0 * cos_lat1 * cos_ldiff)), &
    &                  (sin_lat0 * sin_lat1 + cos_lat0 * cos_lat1 * cos_ldiff))

    c = atan2(cos_lat1 * sin_ldiff, cos_lat0 * sin_lat1 - sin_lat0 * cos_lat1 * cos_ldiff)
    if(c .lt. 0.0_fp) c = c + 2.0_fp * pi
    b = pi + atan2(cos_lat0 * sin_ldiff, -sin_lat0 * cos_lat1 + cos_lat0 * sin_lat1 * cos_ldiff)
    if(b .ge. 2.0_fp * pi) b = b - 2.0_fp * pi

    if(present(delta_out)) delta_out = delta
    if(present(distance)) distance = r_earth * delta
    if(present(azimuth)) azimuth = c
    if(present(backazimuth)) backazimuth = b

    return
  end subroutine greatcircle_dist

  subroutine latgtoc(lat_in, lat_out)
    use nrtype, only : fp
    use constants, only : r2r1, deg2rad
    real(kind = fp), intent(in) :: lat_in    !!degree
    real(kind = fp), intent(out) :: lat_out  !!radian

    lat_out = atan(r2r1 * r2r1 * tan(lat_in * deg2rad))
    return
  end subroutine latgtoc

  subroutine latctog(lat_in, lat_out)
    use nrtype, only : fp
    use constants, only : r2r1, deg2rad, rad2deg
    real(kind = fp), intent(in) :: lat_in    !!radian
    real(kind = fp), intent(out) :: lat_out  !!degree

    lat_out = atan(tan(lat_in) / (r2r1 * r2r1)) * rad2deg
    return
  end subroutine latctog

end module greatcircle

max_stars_repo_name: mogiso/AmplitudeSourceLocation

max_issues_repo_name: mogiso/AmplitudeSourceLocation

max_forks_repo_name: mogiso/AmplitudeSourceLocation


text:
MODULE integrators
USE environment_settings
USE rad_decay_class
USE rad_decay_chain
IMPLICIT NONE
PRIVATE

	TYPE, PUBLIC :: Integrator
		TYPE(DecayChain)             :: chain
		REAL(kind=rknd), ALLOCATABLE :: ynew(:)
		REAL(kind=rknd), ALLOCATABLE :: y(:)
		REAL(kind=rknd), ALLOCATABLE :: yp(:)
		REAL(kind=rknd)              :: relerr, abserr
		REAL(kind=rknd)              :: t, ti, tf, dt
		INTEGER(kind=iknd)           :: stepmax, steps, neqn
		CHARACTER(LEN=100)           :: intg_type
		CONTAINS
			procedure :: init => initialize_integrator
			procedure :: f    => function
			procedure :: step => integration_step
	END TYPE

CONTAINS

	SUBROUTINE initialize_integrator(this, chain, ti, tf, dt, stepmax, relerr, abserr, intg_type)
		CLASS(Integrator)  :: this
		TYPE(DecayChain)   :: chain
		REAL(kind=rknd)    :: relerr, abserr
		REAL(kind=rknd)    :: ti, tf, dt
		INTEGER(kind=iknd) :: stepmax
		CHARACTER(LEN=*)   :: intg_type
		INTEGER(kind=iknd) :: k,j,i
		this%chain     = chain
		this%t         = ti
		this%ti        = ti
		this%tf        = tf
		this%dt        = dt
		this%stepmax   = stepmax
		this%relerr    = relerr
		this%abserr    = abserr
		this%intg_type = intg_type
		this%neqn      = this%chain%nisos
		ALLOCATE(this%ynew(this%neqn))
		ALLOCATE(this%y(this%neqn))
		ALLOCATE(this%yp(this%neqn))
		this%yp   = 0.0d0
		this%ynew = 0.0d0
		DO k = 1, this%neqn
			this%y(k) = this%chain%isotope_list(k)%c0
		ENDDO
		this%steps = 0
	END SUBROUTINE initialize_integrator

	SUBROUTINE function(this)
	IMPLICIT NONE
		CLASS(Integrator)  :: this
		INTEGER(KIND=iknd) :: k, j, kd, jd
		REAL(KIND=rknd)    :: fx, lm
		this%yp = 0.0d0
		DO k = 1, this%neqn
			IF (this%chain%isotope_list(k)%t_half .le. 0.0d0) CYCLE
			kd = this%chain%isotope_list(k)%ID
			lm = this%chain%isotope_list(k)%lambda
			fx = lm * this%y(kd)
			this%yp(kd) = this%yp(kd) - fx
			DO j = 1, this%chain%isotope_list(k)%nprods
				jd = this%chain%isotope_list(k)%prodID(j)
				this%yp(jd) = this%yp(jd) + fx
			ENDDO
		ENDDO
	END SUBROUTINE function

	SUBROUTINE integration_step(this)
		CLASS(Integrator)  :: this
		INTEGER(KIND=iknd) :: k,j,i
		CALL this%f()
		IF (this%tf - this%t .lt. this%dt) this%dt = this%tf - this%t
		SELECT CASE (TRIM(this%intg_type))

		CASE('forward_euler')
			CALL forward_euler(this%neqn, this%y, this%yp, this%t, this%dt)
		END SELECT

		DO k = 1, this%chain%nisos
			j = this%chain%isotope_list(k)%ID
			this%chain%isotope_list(k)%c = this%y(j)
		ENDDO
		this%steps = this%steps + 1
	END SUBROUTINE integration_step

	SUBROUTINE forward_euler(n,y,yp,t,dt)
	IMPLICIT NONE
		INTEGER(KIND=iknd), intent(IN) :: n
		REAL(KIND=rknd), intent(INOUT) :: y(n), t
		REAL(KIND=rknd),    intent(IN) :: yp(n), dt
		INTEGER(KIND=iknd)             :: i
		t = t + dt
		DO i = 1, n
			y(i) = y(i) + dt*yp(i)
		ENDDO
	END SUBROUTINE forward_euler

END MODULE integrators

max_stars_repo_name: ebrooker/RadDecay

max_issues_repo_name: ebrooker/RadDecay

max_forks_repo_name: ebrooker/RadDecay


text:
       module compfrac
!***********************************************************************
! $Id: compfrac.f,v 1.1 2006/05/17 15:23:20 zvd Exp $
!***********************************************************************
!  Copyright, 2002, 2004,  The  Regents of the University of California.
!  This  program  was  prepared by  the  Regents  of the  University  of
!  California  at Los Alamos National Laboratory (the  University) under
!  contract No. W-7405-ENG-36 with the U. S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is  granted to the  public to copy  and  use this software
!  without  charge,  provided that this  Notice  and  any  statement  of
!  authorship are reproduced on all copies. Neither the U. S. Government
!  nor  the  University  makes any  warranty,  express  or  implied,  or
!  assumes any liability or responsibility for the use of this software.
!***********************************************************************
!     Module containing variables for parallel fracture dispersion
!     interpolation calculations. Developed originally for FEHM.
!
!     Identifier      Type     Description
!
!     numparams       integer  Number of parameters defining type curves
!     nump1           integer  Number of param1 values used to define
!                                transfer function curves.
!     nump2           integer  Number of param2 values used to define
!                                transfer function curves.
!     nump3           integer  Number of param3 values used to define
!                                transfer function curves.
!     d4              integer  Fracture-matrix flow interaction flag.
!     curve_structure integer  If 1, free format structure of curves
!     log_flag        integer  If 0, use value itself in interpolation,
!                                if 1, take log of value first
!     normal_param    real*4   Array of values for normalizing the
!                                parameters before performing the
!                                nearest neighbor or multi-curve
!                                interpolation.
!     nump            integer  Number of time/concentration pairs in
!                                each curve, array of size
!     wt_curve        real*4   array of size
!     weight_factor   real*4   Optional weight factor used when
!                                curve_structure > 1.
!     param1          real*8   Parameter 1 values used to define
!                                transfer function curves, array of
!                                size nump1
!     param2          real*8   Parameter 2 values used to define
!                                transfer function curves, array of
!                                size nump2
!     param3          real*8   Parameter 3 values used to define
!                                transfer function curves, array of
!                                size nump3
!     dtime           real*8   Transfer function curve delay times,
!                                array of size (nump1, nump2, nump3,
!                                d4, nump).
!     conc            real*8   Transfer function curve concentrations,
!                                array of size (nump1, nump2, nump3,
!                                d4, nump).
!     param_density   real*8   array of size
!     pfrac_read      logical  Flag indicating if data has been read
!     ipartout        integer
!***********************************************************************

      implicit none
      save

      integer :: numparams
      integer :: nump1
      integer :: nump2
      integer :: nump3
      integer :: d4
      integer :: curve_structure
      integer :: log_flag(100)
      real :: normal_param(100)
      integer, allocatable, dimension(:,:,:,:) :: nump
      integer, allocatable, dimension(:,:,:) :: itf_curve
      real, allocatable, dimension(:,:,:) :: wt_curve
      real weight_factor
      real*8 sigma_low, sigma_high, omega_low, omega_high
      real*8, allocatable, dimension(:) :: param1
      real*8, allocatable, dimension(:) :: param2
      real*8, allocatable, dimension(:) :: param3
      real*8, allocatable, dimension(:,:,:,:,:) :: dtime
      real*8, allocatable, dimension(:,:,:,:,:) :: conc
      integer, allocatable, dimension(:,:,:) :: param_density

      logical :: pfrac_read = .false.
      logical :: error_flag = .false.
      integer :: ipartout = 0

      end module compfrac

max_stars_repo_name: lanl/PLUMECALC

max_issues_repo_name: lanl/PLUMECALC

max_forks_repo_name: lanl/PLUMECALC


! { dg-do run }
!
! Tests the fix for PR82622 comment #1, where the declaration of
! 'x' choked during initialization. Once fixed, it was found that
! IO was not working correctly for PDT array components.
!
! Contributed by Gerhard Steinmetz  <gscfq@t-online.de>
!
program p
   character(120) :: buffer
   integer :: i(4)
   type t(a)
      integer, len :: a
   end type
   type t2(b)
      integer, len :: b
      type(t(1)) :: r(b)
   end type
   type(t2(3)) :: x
   write (buffer,*) x
   read (buffer,*) i
   if (any (i .ne. [3,1,1,1])) STOP 1
end


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Copyright (c) 2018-2019, Lawrence Livermore National Security, LLC.
! Produced at the Lawrence Livermore National Laboratory
!
! Created by David Beckingsale, david@llnl.gov
! LLNL-CODE-747640
!
! All rights reserved.
!
! This file is part of Umpire.
!
! For details, see https://github.com/LLNL/Umpire
! Please also see the LICENSE file for MIT license.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program fortran_test
      use umpire_mod
      implicit none
      logical ok

      integer(C_INT), pointer :: array(:)
      type(UmpireAllocator) allocator
      type(UmpireResourceManager) rm

      rm = rm%get_instance()
      allocator = rm%get_allocator_by_id(0)

      call allocator%allocate(array, [ 10 ])

      write(10,*) "Allocated array of size ", 10

      call allocator%deallocate(array)

      write(10,*) "deallocated."
end program fortran_test

max_stars_repo_name: aaroncblack/Umpire

max_issues_repo_name: aaroncblack/Umpire

max_forks_repo_name: aaroncblack/Umpire


!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006-2018
!        Salvatore Filippone
!        Alfredo Buttari
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
program s_file_spmv
  use psb_base_mod
  use psb_util_mod
  implicit none

  ! input parameters
  character(len=40) :: kmethd, ptype, mtrx_file, rhs_file

  ! sparse matrices
  type(psb_sspmat_type) :: a
  type(psb_lsspmat_type) :: aux_a

  ! dense matrices
  real(psb_spk_), allocatable, target ::  aux_b(:,:), d(:)
  real(psb_spk_), allocatable , save  :: x_col_glob(:), r_col_glob(:)
  real(psb_spk_), pointer  :: b_col_glob(:)
  type(psb_s_vect_type)    :: b_col, x_col, r_col

  ! communications data structure
  type(psb_desc_type):: desc_a

  type(psb_ctxt_type) :: ctxt
  integer(psb_ipk_) :: iam, np

  ! solver paramters
  integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,&
       & methd, istopc, irst, nr
  integer(psb_lpk_) :: lnp
  integer(psb_epk_) :: amatsize, descsize, annz, nbytes
  real(psb_spk_)   :: err, eps,cond

  character(len=5)   :: afmt
  character(len=20)  :: name
  character(len=2)   :: filefmt
  integer(psb_ipk_), parameter :: iunit=12
  integer(psb_ipk_), parameter :: times=20
  integer(psb_ipk_) :: iparm(20)

  ! other variables
  integer(psb_lpk_) :: i,j,m_problem
  integer(psb_ipk_) :: internal, m,ii,nnzero,info
  real(psb_dpk_) :: t1, t2, r_amax, b_amax,&
       &scale,resmx,resmxp, flops, bdwdth
  real(psb_dpk_) :: tt1, tt2, tflops
  integer(psb_ipk_) :: nrhs, nrow, n_row, dim, nv, ne
  integer(psb_ipk_), allocatable :: ivg(:), ipv(:)

  call psb_init(ctxt)
  call psb_info(ctxt,iam,np)

  if (iam < 0) then
    ! This should not happen, but just in case
    call psb_exit(ctxt)
    stop
  endif

  name='s_file_spmv'
  if(psb_get_errstatus() /= 0) goto 9999
  info=psb_success_
  call psb_set_errverbosity(2)
  !
  ! Hello world
  !
  if (iam == psb_root_) then
    write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_
    write(*,*) 'This is the ',trim(name),' sample program'
    read(psb_inp_unit,*) mtrx_file
    read(psb_inp_unit,*) filefmt
    read(psb_inp_unit,*) ipart
  end if
  call psb_bcast(ctxt,mtrx_file)
  call psb_bcast(ctxt,filefmt)
  call psb_bcast(ctxt,ipart)
  rhs_file = 'NONE'
  afmt     = 'CSR'
  call psb_barrier(ctxt)
  t1 = psb_wtime()
  ! read the input matrix to be processed and (possibly) the rhs
  nrhs = 1

  if (iam==psb_root_) then
    select case(psb_toupper(filefmt))
    case('MM')
      ! For Matrix Market we have an input file for the matrix
      ! and an (optional) second file for the RHS.
      call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file)
      if (info == psb_success_) then
        if (rhs_file /= 'NONE') then
          call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file)
        end if
      end if

    case ('HB')
      ! For Harwell-Boeing we have a single file which may or may not
      ! contain an RHS.
      call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file)

    case default
      info = -1
      write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt
    end select
    if (info /= psb_success_) then
      write(psb_err_unit,*) 'Error while reading input matrix '
      call psb_abort(ctxt)
    end if

    m_problem = aux_a%get_nrows()
    call psb_bcast(ctxt,m_problem)

    ! At this point aux_b may still be unallocated
    if (psb_size(aux_b,dim=1)==m_problem) then
      ! if any rhs were present, broadcast the first one
      write(psb_err_unit,'("Ok, got an rhs ")')
      b_col_glob =>aux_b(:,1)
    else
      write(psb_out_unit,'("Generating an rhs...")')
      write(psb_out_unit,'(" ")')
      call psb_realloc(m_problem,1,aux_b,ircode)
      if (ircode /= 0) then
        call psb_errpush(psb_err_alloc_dealloc_,name)
        goto 9999
      endif

      b_col_glob => aux_b(:,1)
      do i=1, m_problem
        b_col_glob(i) = 1.d0
      enddo
    endif

  else

    call psb_bcast(ctxt,m_problem)
    b_col_glob =>aux_b(:,1)

  end if

  ! switch over different partition types
  write(psb_out_unit,'("Number of processors : ",i0)')np
  if (ipart == 0) then
    call psb_barrier(ctxt)
    if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
    allocate(ivg(m_problem),ipv(np))
    do i=1,m_problem
      call part_block(i,m_problem,np,ipv,nv)
      ivg(i) = ipv(1)
    enddo
    call psb_matdist(aux_a, a, ctxt,desc_a,info,fmt=afmt,vg=ivg)

  else if (ipart == 2) then
    if (iam==psb_root_) then
      write(psb_out_unit,'("Partition type: graph")')
      write(psb_out_unit,'(" ")')
      !      write(psb_err_unit,'("Build type: graph")')
      call build_mtpart(aux_a,lnp)

    endif
    call psb_barrier(ctxt)
    call distr_mtpart(psb_root_,ctxt)
    call getv_mtpart(ivg)
    call psb_matdist(aux_a, a, ctxt, desc_a,info,fmt=afmt,vg=ivg)

  else
    if (iam==psb_root_) write(psb_out_unit,'("Partition type: default block")')
    call psb_matdist(aux_a, a,  ctxt, desc_a,info,fmt=afmt,parts=part_block)
  end if

  call psb_geall(x_col,desc_a,info)
  call x_col%set(sone)
  call psb_geasb(x_col,desc_a,info)
  call psb_geall(b_col,desc_a,info)
  call x_col%zero()
  call psb_geasb(b_col,desc_a,info)
  t2 = psb_wtime() - t1

  call psb_amx(ctxt, t2)

  if (iam==psb_root_) then
    write(psb_out_unit,'(" ")')
    write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2
    write(psb_out_unit,'(" ")')
  end if

  call psb_barrier(ctxt)
  t1 = psb_wtime()
  do i=1,times
    call psb_spmm(sone,a,x_col,szero,b_col,desc_a,info,'n')
  end do
  call psb_barrier(ctxt)
  t2 = psb_wtime() - t1
  call psb_amx(ctxt,t2)

  ! FIXME: cache flush needed here
  call psb_barrier(ctxt)
  tt1 = psb_wtime()
  do i=1,times
    call psb_spmm(sone,a,x_col,szero,b_col,desc_a,info,'t')
  end do
  call psb_barrier(ctxt)
  tt2 = psb_wtime() - tt1
  call psb_amx(ctxt,tt2)

  nr       = desc_a%get_global_rows()
  annz     = a%get_nzeros()
  amatsize = psb_sizeof(a)
  descsize = psb_sizeof(desc_a)
  call psb_sum(ctxt,annz)
  call psb_sum(ctxt,amatsize)
  call psb_sum(ctxt,descsize)

  if (iam==psb_root_) then
    flops = 2.d0*times*annz
    tflops=flops
    write(psb_out_unit,'("Matrix: ",a)') mtrx_file
    write(psb_out_unit,'("Test on                          : ",i20," processors")') np
    write(psb_out_unit,'("Size of matrix                   : ",i20,"           ")') nr
    write(psb_out_unit,'("Number of nonzeros               : ",i20,"           ")') annz
    write(psb_out_unit,'("Memory occupation                : ",i20,"           ")') amatsize
    write(psb_out_unit,'("Number of flops (",i0," prod)        : ",F20.0,"           ")') times,flops
    flops = flops / (t2)
    tflops = tflops / (tt2)
    write(psb_out_unit,'("Time for ",i0," products (s)         : ",F20.3)')times, t2
    write(psb_out_unit,'("Time per product    (ms)         : ",F20.3)') t2*1.d3/(1.d0*times)
    write(psb_out_unit,'("MFLOPS                           : ",F20.3)') flops/1.d6

    write(psb_out_unit,'("Time for ",i0," products (s) (trans.): ",F20.3)') times,tt2
    write(psb_out_unit,'("Time per product    (ms) (trans.): ",F20.3)') tt2*1.d3/(1.d0*times)
    write(psb_out_unit,'("MFLOPS                   (trans.): ",F20.3)') tflops/1.d6

    !
    ! This computation is valid for CSR
    !
    nbytes = nr*(2*psb_sizeof_sp + psb_sizeof_ip)+ &
         & annz*(psb_sizeof_sp + psb_sizeof_ip)
    bdwdth = times*nbytes/(t2*1.d6)
    write(psb_out_unit,*)
    write(psb_out_unit,'("MBYTES/S                         : ",F20.3)') bdwdth
    bdwdth = times*nbytes/(tt2*1.d6)
    write(psb_out_unit,'("MBYTES/S                  (trans): ",F20.3)') bdwdth

  end if

  call psb_gefree(b_col, desc_a,info)
  call psb_gefree(x_col, desc_a,info)
  call psb_spfree(a, desc_a,info)
  call psb_cdfree(desc_a,info)

  call psb_exit(ctxt)
  stop

9999 call psb_error(ctxt)

  stop

end program s_file_spmv