text:
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     https://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

module test_helpers

  use, intrinsic :: iso_c_binding, only: c_double, c_int, c_bool
  use helpers, only: &
       WIGGLE, cross_product, bbox, wiggle_interval, contains_nd, &
       vector_close, in_interval, in_sorted, convex_hull, polygon_collide, &
       solve2x2
  use types, only: dp
  use unit_test_helpers, only: MACHINE_EPS, print_status
  implicit none
  private &
       test_cross_product, test_bbox, test_wiggle_interval, &
       test_contains_nd, test_vector_close, test_in_interval, &
       test_in_sorted, test_convex_hull, test_polygon_collide, test_solve2x2
  public helpers_all_tests

contains

  subroutine helpers_all_tests(success)
    logical(c_bool), intent(inout) :: success

    call test_cross_product(success)
    call test_bbox(success)
    call test_wiggle_interval(success)
    call test_contains_nd(success)
    call test_vector_close(success)
    call test_in_interval(success)
    call test_in_sorted(success)
    call test_convex_hull(success)
    call test_polygon_collide(success)
    call test_solve2x2(success)

  end subroutine helpers_all_tests

  subroutine test_cross_product(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: vec0(2)
    real(c_double) :: vec1(2)
    real(c_double) :: result_
    integer :: case_id
    character(13) :: name

    case_id = 1
    name = "cross_product"

    ! CASE 1: Identical vector.
    vec0 = [1.0_dp, 7.0_dp] / 8
    vec1 = [-11.0_dp, 24.0_dp] / 32
    call cross_product(vec0, vec1, result_)
    case_success = (result_ == 101.0_dp / 256)
    call print_status(name, case_id, case_success, success)

  end subroutine test_cross_product

  subroutine test_bbox(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: nodes1(2, 2), nodes2(2, 6)
    real(c_double) :: left, right, bottom, top
    integer :: case_id
    character(4) :: name

    case_id = 1
    name = "bbox"

    ! CASE 1: Simple case (just two nodes).
    nodes1(:, 1) = [0.0_dp, 5.0_dp]
    nodes1(:, 2) = [1.0_dp, 3.0_dp]
    call bbox(2, nodes1, left, right, bottom, top)
    case_success = ( &
         left == 0.0_dp .AND. right == 1.0_dp &
         .AND. bottom == 3.0_dp .AND. top == 5.0_dp)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: "Many" nodes.
    nodes2(:, 1) = [1.0_dp, 0.0_dp]
    nodes2(:, 2) = [2.0_dp, 1.0_dp]
    nodes2(:, 3) = [-1.0_dp, 2.0_dp]
    nodes2(:, 4) = [5.0_dp, -3.0_dp]
    nodes2(:, 5) = [4.0_dp, 4.0_dp]
    nodes2(:, 6) = [0.0_dp, 0.0_dp]
    call bbox(6, nodes2, left, right, bottom, top)
    case_success = ( &
         left == -1.0_dp .AND. right == 5.0_dp &
         .AND. bottom == -3.0_dp .AND. top == 4.0_dp)
    call print_status(name, case_id, case_success, success)

  end subroutine test_bbox

  subroutine test_wiggle_interval(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: value1, value2, value3, value4, value5
    real(c_double) :: result1, result2, result3, result4, result5
    logical(c_bool) :: &
         wiggle_success1, wiggle_success2, wiggle_success3, &
         wiggle_success4, wiggle_success5
    integer :: case_id
    character(15) :: name

    case_id = 1
    name = "wiggle_interval"

    ! CASE 1: **Exactly** at endpoint.
    value1 = 0.0_dp
    call wiggle_interval(value1, result1, wiggle_success1)
    case_success = (wiggle_success1 .AND. result1 == 0.0_dp)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Near endpoint.
    value1 = 1.0_dp + 0.5_dp**20
    call wiggle_interval(value1, result1, wiggle_success1)
    case_success = (.NOT. wiggle_success1)
    call print_status(name, case_id, case_success, success)

    ! CASE 3: "Below" the interval.
    value1 = -0.25_dp
    call wiggle_interval(value1, result1, wiggle_success1)
    case_success = (.NOT. wiggle_success1)
    call print_status(name, case_id, case_success, success)

    ! CASE 4: "Above" the interval.
    value1 = 1.5_dp
    call wiggle_interval(value1, result1, wiggle_success1)
    case_success = (.NOT. wiggle_success1)
    call print_status(name, case_id, case_success, success)

    ! CASE 5: Inside the interval.
    value1 = 0.25_dp
    call wiggle_interval(value1, result1, wiggle_success1)
    case_success = (wiggle_success1 .AND. result1 == value1)
    call print_status(name, case_id, case_success, success)

    ! CASE 6: Wiggle "below" the interval.
    value1 = -0.5_dp**60
    call wiggle_interval(value1, result1, wiggle_success1)
    case_success = (wiggle_success1 .AND. result1 == 0.0_dp)
    call print_status(name, case_id, case_success, success)

    ! CASE 7: Wiggle "above" the interval.
    value1 = 1.0_dp + MACHINE_EPS
    call wiggle_interval(value1, result1, wiggle_success1)
    case_success = (wiggle_success1 .AND. result1 == 1.0_dp)
    call print_status(name, case_id, case_success, success)

    ! CASE 8: Test "lower outer" boundary, i.e. values in (-epsilon, 0).
    value1 = -WIGGLE + MACHINE_EPS * WIGGLE
    value2 = -WIGGLE + MACHINE_EPS * WIGGLE / 2
    value3 = -WIGGLE
    value4 = -WIGGLE - MACHINE_EPS * WIGGLE
    call wiggle_interval(value1, result1, wiggle_success1)
    call wiggle_interval(value2, result2, wiggle_success2)
    call wiggle_interval(value3, result3, wiggle_success3)
    call wiggle_interval(value4, result4, wiggle_success4)
    case_success = ( &
         wiggle_success1 .AND. result1 == 0.0_dp .AND. &
         wiggle_success2 .AND. result2 == 0.0_dp .AND. &
         .NOT. wiggle_success3 .AND. .NOT. wiggle_success4)
    call print_status(name, case_id, case_success, success)

    ! CASE 9: Test "upper outer" boundary, i.e. values in (1, 1 + epsilon).
    value1 = 1.0_dp + WIGGLE - 2 * MACHINE_EPS
    value2 = 1.0_dp + WIGGLE - MACHINE_EPS
    value3 = 1.0_dp + WIGGLE
    value4 = 1.0_dp + WIGGLE + MACHINE_EPS
    call wiggle_interval(value1, result1, wiggle_success1)
    call wiggle_interval(value2, result2, wiggle_success2)
    call wiggle_interval(value3, result3, wiggle_success3)
    call wiggle_interval(value4, result4, wiggle_success4)
    case_success = ( &
         wiggle_success1 .AND. result1 == 1.0_dp .AND. &
         wiggle_success2 .AND. result2 == 1.0_dp .AND. &
         .NOT. wiggle_success3 .AND. .NOT. wiggle_success4)
    call print_status(name, case_id, case_success, success)

    ! CASE 10: Test "lower inner" boundary, i.e. values in (0, epsilon).
    value1 = WIGGLE - WIGGLE * MACHINE_EPS
    value2 = WIGGLE - WIGGLE * MACHINE_EPS / 2
    value3 = WIGGLE
    value4 = WIGGLE + WIGGLE * MACHINE_EPS
    call wiggle_interval(value1, result1, wiggle_success1)
    call wiggle_interval(value2, result2, wiggle_success2)
    call wiggle_interval(value3, result3, wiggle_success3)
    call wiggle_interval(value4, result4, wiggle_success4)
    case_success = ( &
         wiggle_success1 .AND. result1 == 0.0_dp .AND. &
         wiggle_success2 .AND. result2 == 0.0_dp .AND. &
         wiggle_success3 .AND. result3 == value3 .AND. &
         wiggle_success4 .AND. result4 == value4)
    call print_status(name, case_id, case_success, success)

    ! CASE 11: Test "uper inner" boundary, i.e. values in (1 - epsilon, 1).
    value1 = 1.0_dp - WIGGLE - MACHINE_EPS
    value2 = 1.0_dp - WIGGLE - MACHINE_EPS / 2
    value3 = 1.0_dp - WIGGLE
    value4 = 1.0_dp - WIGGLE + MACHINE_EPS / 2
    value5 = 1.0_dp - WIGGLE + MACHINE_EPS
    call wiggle_interval(value1, result1, wiggle_success1)
    call wiggle_interval(value2, result2, wiggle_success2)
    call wiggle_interval(value3, result3, wiggle_success3)
    call wiggle_interval(value4, result4, wiggle_success4)
    call wiggle_interval(value5, result5, wiggle_success5)
    case_success = ( &
         wiggle_success1 .AND. result1 == value1 .AND. &
         wiggle_success2 .AND. result2 == value2 .AND. &
         wiggle_success3 .AND. result3 == value3 .AND. &
         wiggle_success4 .AND. result4 == 1.0_dp .AND. &
         wiggle_success5 .AND. result5 == 1.0_dp)
    call print_status(name, case_id, case_success, success)

  end subroutine test_wiggle_interval

  subroutine test_contains_nd(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: nodes1(2, 3), nodes2(3, 2), nodes3(4, 2)
    real(c_double) :: point1(2), point2(3), point3(4)
    logical(c_bool) :: is_contained
    integer :: case_id
    character(11) :: name

    case_id = 1
    name = "contains_nd"

    ! CASE 1: Below bounding box.
    nodes1(:, 1) = [0.0_dp, 1.0_dp]
    nodes1(:, 2) = [0.5_dp, 0.0_dp]
    nodes1(:, 3) = [1.0_dp, 2.0_dp]
    point1 = [-0.5_dp, 1.0_dp]
    call contains_nd(3, 2, nodes1, point1, is_contained)
    case_success = (.NOT. is_contained)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Above bounding box.
    nodes2(:, 1) = [0.0_dp, -4.0_dp, 2.0_dp]
    nodes2(:, 2) = [-1.0_dp, 1.0_dp, 3.0_dp]
    point2 = [-0.5_dp, 2.0_dp, 2.5_dp]
    call contains_nd(2, 3, nodes2, point2, is_contained)
    case_success = (.NOT. is_contained)
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Inside bounding box.
    nodes3(:, 1) = [0.0_dp, 1.0_dp, 2.0_dp, 3.0_dp]
    nodes3(:, 2) = [1.0_dp, -2.0_dp, -4.0_dp, 1.0_dp]
    point3 = [0.5_dp, 0.0_dp, 0.0_dp, 2.0_dp]
    call contains_nd(2, 4, nodes3, point3, is_contained)
    case_success = (is_contained)
    call print_status(name, case_id, case_success, success)

  end subroutine test_contains_nd

  subroutine test_vector_close(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: eps
    real(c_double) :: vec1(2)
    real(c_double) :: vec2(2)
    logical(c_bool) :: is_close
    integer :: case_id
    character(12) :: name

    eps = 0.5_dp**40
    case_id = 1
    name = "vector_close"

    ! CASE 1: Identical vector.
    vec1 = [0.5_dp, 4.0_dp]
    is_close = vector_close(2, vec1, vec1, eps)
    case_success = (is_close)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Far apart vectors.
    vec1 = [0.0_dp, 6.0_dp]
    vec2 = [1.0_dp, -4.0_dp]
    is_close = vector_close(2, vec1, vec2, eps)
    case_success = (.NOT. is_close)
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Close but different.
    vec1 = [2.25_dp, -3.5_dp]
    vec2 = vec1 + 0.5_dp**43 * [-5.0_dp, 12.0_dp]
    is_close = vector_close(2, vec1, vec2, eps)
    case_success = (is_close)
    call print_status(name, case_id, case_success, success)

    ! CASE 4: Custom epsilon.
    vec1 = [3.0_dp, 4.0_dp]
    vec2 = [2.0_dp, 5.0_dp]
    is_close = vector_close(2, vec1, vec2, 0.5_dp)
    case_success = (is_close .AND. .NOT. vector_close(2, vec1, vec2, eps))
    call print_status(name, case_id, case_success, success)

    ! CASE 5: Near zero.
    vec1 = [0.0_dp, 0.0_dp]
    vec2 = 0.5_dp**45 * [3.0_dp, 4.0_dp]
    is_close = vector_close(2, vec1, vec2, eps)
    case_success = (is_close)
    call print_status(name, case_id, case_success, success)

    ! CASE 6: Near zero failure (i.e. not near enough).
    vec1 = 0.5_dp**20 * [1.0_dp, 0.0_dp]
    vec2 = [0.0_dp, 0.0_dp]
    is_close = vector_close(2, vec1, vec2, eps)
    case_success = (.NOT. is_close)
    call print_status(name, case_id, case_success, success)

  end subroutine test_vector_close

  subroutine test_in_interval(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    logical(c_bool) :: is_inside
    integer :: case_id
    character(11) :: name

    case_id = 1
    name = "in_interval"

    ! CASE 1: Interior value.
    is_inside = in_interval(1.5_dp, 1.0_dp, 2.0_dp)
    case_success = (is_inside)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Barely inside.
    is_inside = in_interval(1.0_dp + MACHINE_EPS, 1.0_dp, 2.0_dp)
    case_success = (is_inside)
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Barely outside.
    is_inside = in_interval(1.0_dp - MACHINE_EPS / 2, 1.0_dp, 2.0_dp)
    case_success = (.NOT. is_inside)
    call print_status(name, case_id, case_success, success)

    ! CASE 4: Exterior value.
    is_inside = in_interval(-1.0_dp, 1.0_dp, 2.0_dp)
    case_success = (.NOT. is_inside)
    call print_status(name, case_id, case_success, success)

  end subroutine test_in_interval

  subroutine test_in_sorted(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    integer(c_int) :: values(5)
    integer :: case_id
    character(9) :: name

    case_id = 1
    name = "in_sorted"

    ! CASE 1: Inside the interval.
    values = [0, 5, 8, 12, 17]
    case_success = ( &
         .NOT. in_sorted(5, values, 4) .AND. &
         in_sorted(5, values, 5) .AND. &
         .NOT. in_sorted(5, values, 6) .AND. &
         .NOT. in_sorted(5, values, 7) .AND. &
         in_sorted(5, values, 8) .AND. &
         .NOT. in_sorted(5, values, 9) .AND. &
         .NOT. in_sorted(5, values, 10) .AND. &
         .NOT. in_sorted(5, values, 11) .AND. &
         in_sorted(5, values, 12) .AND. &
         .NOT. in_sorted(5, values, 13))
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Near the left endpoint.
    values(:3) = [9, 11, 220]
    case_success = ( &
         .NOT. in_sorted(3, values(:3), 7) .AND. &
         .NOT. in_sorted(3, values(:3), 8) .AND. &
         in_sorted(3, values(:3), 9) .AND. &
         .NOT. in_sorted(3, values(:3), 10))
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Near the right endpoint.
    values(:3) = [16, 18, 20]
    case_success = ( &
         .NOT. in_sorted(3, values(:3), 19) .AND. &
         in_sorted(3, values(:3), 20) .AND. &
         .NOT. in_sorted(3, values(:3), 21) .AND. &
         .NOT. in_sorted(3, values(:3), 22))
    call print_status(name, case_id, case_success, success)

  end subroutine test_in_sorted

  subroutine test_convex_hull(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    integer :: case_id
    real(c_double) :: points1(2, 5), points2(2, 2), points3(2, 100)
    integer(c_int) :: polygon_size, i, j
    real(c_double) :: polygon(2, 100)
    character(11) :: name

    case_id = 1
    name = "convex_hull"

    ! CASE 1: Triangle with centroid (i.e. inside) and a repeated corner.
    points1(1, :) = [0.0_dp, 0.0_dp, 1.0_dp, 3.0_dp, 0.0_dp]
    points1(2, :) = [0.0_dp, 3.0_dp, 1.0_dp, 0.0_dp, 3.0_dp]
    call convex_hull(5, points1, polygon_size, polygon(:, :5))
    case_success = ( &
         polygon_size == 3 .AND. &
         all(polygon(1, :3) == [0.0_dp, 3.0_dp, 0.0_dp]) .AND. &
         all(polygon(2, :3) == [0.0_dp, 0.0_dp, 3.0_dp]))
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Degenerate case (2 points, i.e. the "polygon" is just a line).
    points2(1, :) = [0.0_dp, 1.0_dp]
    points2(2, :) = [0.0_dp, 0.0_dp]
    call convex_hull(2, points2, polygon_size, polygon(:, :2))
    case_success = ( &
         polygon_size == 2 .AND. &
         all(polygon(1, :2) == [0.0_dp, 1.0_dp]) .AND. &
         all(polygon(2, :2) == [0.0_dp, 0.0_dp]))
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Same as CASE 2, but points are in reverse order.
    points2(1, :) = [1.0_dp, 0.0_dp]
    points2(2, :) = [0.0_dp, 0.0_dp]
    call convex_hull(2, points2, polygon_size, polygon(:, :2))
    case_success = ( &
         polygon_size == 2 .AND. &
         all(polygon(1, :2) == [0.0_dp, 1.0_dp]) .AND. &
         all(polygon(2, :2) == [0.0_dp, 0.0_dp]))
    call print_status(name, case_id, case_success, success)

    ! CASE 4: 10-by-10 grid of points. Also covers case where ``polygon`` is
    !         allocated but not large enough.
    do i = 1, 10
       do j = 1, 10
          points3(1, 10 * (i - 1) + j) = j
          points3(2, 10 * (i - 1) + j) = i
       end do
    end do
    call convex_hull(100, points3, polygon_size, polygon)
    case_success = ( &
         polygon_size == 4 .AND. &
         all(polygon(1, :4) == [1.0_dp, 10.0_dp, 10.0_dp, 1.0_dp]) .AND. &
         all(polygon(2, :4) == [1.0_dp, 1.0_dp, 10.0_dp, 10.0_dp]))
    call print_status(name, case_id, case_success, success)

    ! CASE 5: No points (degenerate case).
    call convex_hull(0, points1(:, :0), polygon_size, polygon(:, :0))
    case_success = (polygon_size == 0)
    call print_status(name, case_id, case_success, success)

    ! CASE 6: One point (degenerate case).
    points1(:, 1) = [2.0_dp, 3.0_dp]
    call convex_hull(1, points1(:, :1), polygon_size, polygon(:, :1))
    case_success = ( &
         polygon_size == 1 .AND. &
         all(polygon(1, :1) == [2.0_dp]) .AND. &
         all(polygon(2, :1) == [3.0_dp]))
    call print_status(name, case_id, case_success, success)

    ! CASE 7: Almost linear. In a previous implementation, this caused
    !         an illegal memory access because the middle point of the
    !         "line" was placed in both the upper and lower hull (which used
    !         4 points for the hull when only 3 were allocated).
    points1(1, 1) = -2320058418903481.0_dp * 0.5_dp**54
    points1(1, 2) = -6216142699324449.0_dp * 0.5_dp**56
    points1(1, 3) = -6304103446069947.0_dp * 0.5_dp**57
    points1(2, 1) = -1911920978562235.0_dp * 0.5_dp**55
    points1(2, 2) = -461442255594807.0_dp * 0.5_dp**57
    points1(2, 3) = 6724799403059327.0_dp * 0.5_dp**57

    call convex_hull(3, points1(:, :3), polygon_size, polygon(:, :3))
    case_success = ( &
         polygon_size == 3 .AND. &
         all(polygon(:, :3) == points1(:, :3)))
    call print_status(name, case_id, case_success, success)

  end subroutine test_convex_hull

  subroutine test_polygon_collide(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    integer :: case_id
    real(c_double) :: polygon1(2, 3), polygon2(2, 4), polygon3(2, 4)
    logical(c_bool) :: collision
    character(15) :: name

    case_id = 1
    name = "polygon_collide"

    ! CASE 1: Triangle and square that do not collide.
    polygon1(1, :) = [1.0_dp, 3.0_dp, -2.0_dp]
    polygon1(2, :) = [1.0_dp, 4.0_dp, 2.0_dp]
    polygon2(1, :) = [3.5_dp, 6.5_dp, 6.5_dp, 3.5_dp]
    polygon2(2, :) = [3.0_dp, 3.0_dp, 7.0_dp, 7.0_dp]
    call polygon_collide(3, polygon1, 4, polygon2, collision)
    case_success = (.NOT. collision)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Same as CASE 1 but with polygons swapped.
    call polygon_collide(4, polygon2, 3, polygon1, collision)
    case_success = (.NOT. collision)
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Triangles that do collide.
    polygon1(1, :) = [0.0_dp, 3.0_dp, 0.0_dp]
    polygon1(2, :) = [0.0_dp, 0.0_dp, 3.0_dp]
    polygon2(1, :3) = [1.0_dp, 4.0_dp, 1.0_dp]
    polygon2(2, :3) = [1.0_dp, 1.0_dp, 4.0_dp]
    call polygon_collide(3, polygon1, 3, polygon2(:, :3), collision)
    case_success = (collision)
    call print_status(name, case_id, case_success, success)

    ! CASE 4: Two squares where an edge other than edge 1 in polygon 1 is
    !         the separating line.
    polygon2(1, :) = [1.0_dp, 1.0_dp, 0.0_dp, 0.0_dp]
    polygon2(2, :) = [0.0_dp, 1.0_dp, 1.0_dp, 0.0_dp]
    polygon3(1, :) = [2.0_dp, 3.0_dp, 3.0_dp, 2.0_dp]
    polygon3(2, :) = [0.0_dp, 0.0_dp, 1.0_dp, 1.0_dp]
    call polygon_collide(4, polygon2, 4, polygon3, collision)
    case_success = (.NOT. collision)
    call print_status(name, case_id, case_success, success)

    ! CASE 5: Square and triangle where none of the edges in polygon 1 are
    !         separating lines and an edge other than edge 1 in polygon 2
    !         **is** the separating line.
    polygon2(1, :) = [0.0_dp, 2.0_dp, 2.0_dp, 0.0_dp]
    polygon2(2, :) = [0.0_dp, 0.0_dp, 2.0_dp, 2.0_dp]
    polygon1(1, :) = [1.0_dp, 4.0_dp, 4.0_dp]
    polygon1(2, :) = [4.0_dp, 1.0_dp, 4.0_dp]
    call polygon_collide(4, polygon2, 3, polygon1, collision)
    case_success = (.NOT. collision)
    call print_status(name, case_id, case_success, success)

  end subroutine test_polygon_collide

  subroutine test_solve2x2(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    integer :: case_id
    real(c_double) :: lhs(2, 2), rhs(2), x_val, y_val
    logical(c_bool) :: singular
    character(8) :: name

    case_id = 1
    name = "solve2x2"

    ! CASE 1: Solve without a row-swap.
    lhs(:, 1) = [2.0_dp, 1.0_dp]
    lhs(:, 2) = [3.0_dp, 2.0_dp]
    rhs = [31.0_dp, 19.0_dp]
    call solve2x2(lhs, rhs, singular, x_val, y_val)
    case_success = ( &
         .NOT. singular .AND. &
         x_val == 5.0_dp .AND. &
         y_val == 7.0_dp)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Solve with a row-swap.
    lhs(:, 1) = [1.0_dp, 4.0_dp]
    lhs(:, 2) = [0.0_dp, 1.0_dp]
    rhs = [3.0_dp, 13.0_dp]
    call solve2x2(lhs, rhs, singular, x_val, y_val)
    case_success = ( &
         .NOT. singular .AND. &
         x_val == 3.0_dp .AND. &
         y_val == 1.0_dp)
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Singular due to first column all zero.
    lhs(:, 1) = 0
    call solve2x2(lhs, rhs, singular, x_val, y_val)
    case_success = singular
    call print_status(name, case_id, case_success, success)

    ! CASE 4: Singular without a row-swap.
    lhs(:, 1) = [2.0_dp, 1.0_dp]
    lhs(:, 2) = [4.0_dp, 2.0_dp]
    call solve2x2(lhs, rhs, singular, x_val, y_val)
    case_success = singular
    call print_status(name, case_id, case_success, success)

    ! CASE 5: Singular with a row-swap.
    lhs(:, 1) = [3.0_dp, 12.0_dp]
    lhs(:, 2) = [1.0_dp, 4.0_dp]
    call solve2x2(lhs, rhs, singular, x_val, y_val)
    case_success = singular
    call print_status(name, case_id, case_success, success)

  end subroutine test_solve2x2

end module test_helpers

max_stars_repo_name: dibir-magomedsaygitov/bezier

max_issues_repo_name: dibir-magomedsaygitov/bezier

max_forks_repo_name: dibir-magomedsaygitov/bezier


text:
The Nachtmann Analytical Laboratory (NAL) is equipped to run water and sediment analysis and sampling. They work with government agencies and private industries.

NAL was established in 1979, and has been certified through the State Department of Health (ELAP) laboratory certification program since 1982.  Nachtman continuously participates in EPAs annul Performance Evaluation studies for both Drinking water and Effluent.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C * * * * * * * * * * * * * * * * * * * * * * * * *
C --- DRIVER FOR RADAU5 ON BEAM PROBLEM
C * * * * * * * * * * * * * * * * * * * * * * * * *
        IMPLICIT REAL*8 (A-H,O-Z)
C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN)
        PARAMETER (ND=80,LWORK=4*ND*ND+12*ND+20,LIWORK=3*ND+20)
        COMMON/NNNN/NCOM,NNCOM,NSQ,NQUATR,DELTAS
        DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK)
C -------- END PARAMETER LIST --------
        REAL*4 TARRAY(2),TRESULT
        EXTERNAL FTIGE,SOLOUT
c ------ FILE DE DONNEES ----------
        OPEN(8,FILE='res_rad5')
        REWIND 8
C --- LOOP FOR DIFFERENT TOLERANCES
        NTOLMN=2
        NTOLMX=10
        NTOLDF=4
        NRLOOP=(NTOLMX-NTOLMN)*NTOLDF+1
        TOLST=0.1D0**NTOLMN
        TOLFC=0.1D0**(1.D0/NTOLDF)
        DO 30 NTOL=1,NRLOOP
C ---------- CONSTANTS --------------
        N=40
        NN=2*N
        NCOM=N
        NSQ=N*N
        NQUATR=NSQ*NSQ
        NNCOM=NN
        AN=N
        DELTAS=1.D0/AN
C --- SET DEFAULT VALUES
        DO 12 I=1,20
        WORK(I)=0.D0
  12    IWORK(I)=0
C --------- TUNED VALUES OF WORK PARAMETERS
        WORK(3)=0.1
        WORK(4)=0.3
        WORK(5)=0.99
        WORK(6)=2.0
C --------- INITIAL VALUES -------------
        T=0.D0
        DO 1 I=1,NN
    1   Y(I)=0.D0
        TEND=5.D0
        H=1.D-3
C --- REQUIRED TOLERANCE
        RTOL=TOLST
        ATOL=RTOL
        ITOL=0
C --------- SWITCH FOR OUTPUT --------
        IOUT=0
C --- SECOND ORDER OPTION
           IJAC=0
           IWORK(9)=N
           MLJAC=N
C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
        IMAS=0
        CALL DTIME(TARRAY,TRESULT)
C --- CALL OF THE SUBROUTINE RADAU5
        CALL RADAU5(NN,FTIGE,T,Y,TEND,H,
     &                  RTOL,ATOL,ITOL,
     &                  FTIGE,IJAC,MLJAC,MUJAC,
     &                  FTIGE,IMAS,MLMAS,MUMAS,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
        CALL DTIME(TARRAY,TRESULT)
C --- PRINT SOLUTION
        WRITE(8,9921)(Y(I),I=1,NN)
 9921   FORMAT(1X,F22.16)
        WRITE(8,*)TARRAY(1)
        WRITE (8,*) (IWORK(I),I=14,20)
        WRITE(6,*)' ***** TOL=',RTOL,'  ELAPSED TIME=',TARRAY(1),' ****'
        WRITE (6,91) (IWORK(I),I=14,20)
 91     FORMAT(' fcn=',I5,' jac=',I4,' step=',I4,
     &        ' accpt=',I4,' rejct=',I3,' dec=',I4,
     &        ' sol=',I5)
C -------- NEW TOLERANCE ---
        IF(TARRAY(1).GT.1000.)STOP
        TOLST=TOLST*TOLFC
 30     CONTINUE
        STOP
        END
C
        SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN)
        IMPLICIT REAL*8 (A-H,O-Z)
        DIMENSION Y(N),CONT(LRC)
           WRITE (6,99) X,Y(1),Y(2),NR-1
 99     FORMAT(1X,'X =',F6.2,'    Y =',2F18.7,'    NSTEP =',I4)
        RETURN
        END

max_stars_repo_name: cpmech/hwode

max_issues_repo_name: cpmech/hwode

max_forks_repo_name: cpmech/hwode


text:
subroutine  sslope(sinsolelslope,solelslope,effdec,slb,hrangl)
  use simsphere_mod
  implicit none

  real :: effdec, slb, hrangl
  real :: sinsolelslope,solelslope,dipan,dipaz

!      INCLUDE 'modvars.h'

! subroutine calculates solar elevation angle and azimuth
! for sloping terrain when slope is non-zero.
! elevation for northwest corner of square box is ZNW, etc.
! grid size (called xmeshl) is in same units as ZNW

! calculate dip angle and azimuth angle of sloping surface plane

!      dipaz = 0.0
!      dipan = 0.0

!      edip = (ZNW + ZSW - ZNE - ZSE)/2.0
!      sdip = (ZNW - ZSW + ZNE - ZSE)/2.0

!      if ( edip. ne. 0 ) then

! dipan = abs (atan (edip / xmeshl / (cos(atan(sdip/edip)))))
! dipaz = 3. * 3.14159/ 2. + atan( sdip / edip)

! if ( edip. lt. 0. ) then
!   dipaz = dipaz - 3.14159
! end if

!      else

! dipan = abs (atan (sdip / xmeshl ))
! dipaz = 3.14159

!   if ( sdip. gt. 0. ) then
!     dipaz = 0.0
!   end if

! end if

  dipaz = ASPECT / 57.2958

! if (dipaz. gt. 6.2835) then
!   azimuthangle = azimuthangle  - 360.
! endif

  dipan = SLOPE / 57.2958

!  compute solar elevation angle function for sloping terrain

  sinsolelslope = cos(dipan) * ( sin(slb) * sin(effdec) + cos(slb) *    &
                  cos(effdec) * cos(hrangl)) + sin(dipan) *             &
                  ( cos(dipaz) * ( tan(slb) * ( sin(slb) * sin(effdec) +&
                  cos(slb) * cos(effdec) * cos(hrangl)) - sin(effdec) / &
                  cos(slb)) + sin(dipaz) * cos(effdec) * sin(hrangl))

  solelslope = asin(sinsolelslope)

  if ( solelslope <= 0.01) then
    solelslope = 0.01
  end if

  call albedo(sinsolelslope)

!   note that solar elevation can be less than zero for non-zero slope
!   although greater than zero for flat slope.
!   we use solar elevation for slope in albedo but not for path
!   which depends on solar angle with respect to flat plane

  return
end

max_stars_repo_name: derekmorr/simsphere

max_issues_repo_name: derekmorr/simsphere

max_forks_repo_name: derekmorr/simsphere


text:
!=========================================================================!
!                                                                         !
!                   Module for time integration with diagonal             !
!                      implicite Rosenbrock-Wanner-Methods                !
!                                                                         !
!=========================================================================!
!
MODULE Integration_Mod
  !
  USE Kind_Mod
  USE MPI_Mod
  USE Control_Mod
  USE IO_Mod
  USE Reac_Mod
  USE ChemKinInput_Mod, ONLY: Density
  USE Sparse_Mod, ONLY: Jac_CC
  USE Rosenbrock_Mod
  USE NetCDF_Mod
  USE Meteo_Mod, ONLY: Temp, cv
  IMPLICIT NONE
  !
  TYPE PI_Param
    REAL(dp) :: Kp
    REAL(dp) :: KI
    REAL(dp) :: ThetaMAX
    REAL(dp) :: rho
  END TYPE PI_Param

  INTEGER :: iBar=-1              ! waitbar increment
  TYPE(PI_Param) :: PI_norm, PI_rej
  !
  CONTAINS
  !
  !=======================================================================
  !===================    Time Integration Routine  ======================
  !=======================================================================
  SUBROUTINE Integrate(y_iconc, Temperature, h0, Tspan, Atol, RtolRow, method)
    !--------------------------------------------------------------------
    ! Input:
    !   - y0 ............. Initial vector
    !   - R0 ............. reaction rates at time=Tspan(1)
    !   - Tspan .......... (/ SimulationTimeStart , SimulationTimeEnd /)
    !   - Atol ........... abs. tolerance for gas spc
    !   - RtolRow ........ rel. tolerance for Rosenbrock-Wanner-Method
    !   - method.......... Rosenbrock-Wanner-Method
    !   - PrintSpc ....... print spc PrintSpc(1:3)
    REAL(dp), INTENT(INOUT) :: y_iconc(nspc)
    REAL(dp), INTENT(INOUT) :: h0, Temperature
    REAL(dp) :: Tspan(2)
    REAL(dp) :: Atol(2)
    REAL(dp) :: RtolROW
    CHARACTER(*) :: method
    !-------------------------------------------------------------------
    ! Output:
    !   - Output.......... struct Out (above) contains output statistics
    !   - Unit............ .simul data
    !-------------------------------------------------------------------
    ! Temporary variables:
    !
    REAL(dp) :: Y0(nDIM), Y(nDIM)       ! old, current y vector
    !
    REAL(dp) :: t             ! current time
    REAL(dp) :: timepart
    REAL(dp) :: time_int = 0.0d0

    REAL(dp) :: h, tnew, tmp, hOld
    REAL(dp) :: error, errorOld
    REAL(dp) :: tmp_tb
    REAL(dp) :: actLWC, zen
    REAL(dp) :: wetRad(nFrac)
    INTEGER  :: ierr(1,1)
    REAL(dp) :: ErrVals(nspc)
    !
    INTEGER :: i, k

    INTEGER, PARAMETER :: maxnsteps = 50000

    LOGICAL :: done=.FALSE.
    LOGICAL :: failed
    !
    CHARACTER(10) :: swMethod
    !
    !-----------------------------
    ! LSODE Parameter
    INTEGER :: ITOL , ITASK , ISTATE, IOPT, MF, LIW, LRW
    INTEGER,  ALLOCATABLE :: IWORK(:)
    REAL(dp), ALLOCATABLE :: RWORK(:)
    !
    REAL(dp) :: RTOL1 , ATOL1(nDIM)

    Y0(1:nspc) = y_iconc(:)
    Y(1:nspc)  = y_iconc(:)

    IF ( Teq ) THEN
      !--- initial temperature
      Y0(nDIM) = Temperature     ! = 750 [K] aus speedchem debug
      Y(nDIM)  = Temperature     ! = 750 [K]

      Y0 = MAX( ABS(Y0) , eps ) * SIGN( ONE , Y0 )    ! |y| >= eps
      Y  = MAX( ABS(Y)  , eps ) * SIGN( ONE , Y  )    ! |y| >= eps
    END IF

    t = Tspan(1)
    tnew = Tspan(1)
    h = h0

    ! this is for the waitbar
    tmp_tb = (tEnd-tBegin) * 0.01_dp

    SELECT CASE (TRIM(method(1:7)))

      CASE('METHODS')

        !IF (MPI_master) WRITE(*,'(10X,A)',ADVANCE='NO') 'Start Integration.............      '
        time_int = MPI_WTIME()

        MAIN_LOOP_ROSENBROCK_METHOD: DO

          h = MIN( maxStp, MAX( minStp , h ) )

          !-- Stretch the step if within 5% of tfinal-t.
          IF ( 1.05_dp * h >= Tspan(2) - t ) THEN
            h = ABS(Tspan(2) - t)
            done = .TRUE.
          END IF

          DO    ! Evaluate the formula.

            !-- LOOP FOR ADVANCING ONE STEP.
            failed = .FALSE.      ! no failed attempts

            ! Rosenbrock Timestep
            CALL Rosenbrock(  Y             & ! new concentration
            &               , error         & ! error value
            &               , ierr          & ! max error component
            &               , Y0            & ! current concentration
            &               , t             & ! current time
            &               , h             & ! stepsize
            &               , Euler=.FALSE. ) ! new concentration

            tnew  = t + h

            IF (done) THEN
              tnew  = Tspan(2)    ! Hit end point exactly.
              h     = tnew-t      ! Purify h.
            END IF
            Out%ndecomps   = Out%ndecomps   + 1
            Out%nRateEvals = Out%nRateEvals + ROS%nStage
            Out%nSolves    = Out%nSolves    + ROS%nStage

            failed = (error > ONE)

            IF (failed) THEN      ! failed step
              ! Accept the solution only if the weighted error is no more than the
              ! tolerance one.  Estimate an h that will yield an error of rtol on
              ! the next step or the next try at taking this step, as the case may be,
              ! and use 0.8 of this value to avoid failures.

              Out%nfailed  = Out%nfailed+1
              IF ( h <= minStp ) THEN
                CALL FinishMPI()
                STOP '....Integration_Mod '
              END IF

              h    = MAX( minstp , h * MAX( rTEN , 0.8_dp * error**(-ROS%pow) ) )
              done = .FALSE.
            ELSE                  ! successful step
              EXIT
            END IF
          END DO
          Out%nsteps = Out%nsteps + 1

          ! --- save to NetCDF file
          IF ( NetCdfPrint ) THEN
            IF ( done .OR. StpNetCDF < ZERO ) THEN
              TimeNetCDFA = MPI_WTIME()
              IF (ChemKin) Temperature = Y(nDIM)
              CALL SetOutputNCDF( NetCDF, tnew , h , Y , Temperature )
              CALL StepNetCDF( NetCDF )
              WRITE(999,'(*(Es14.6))') tnew, ( Y(Diag_Index(i)), i=1,SIZE(Diag_Index) )
              TimeNetCDF  = TimeNetCDF + (MPI_WTIME() - TimeNetCDFA)
            END IF
          END IF

          !-- Call progress bar.
          IF ( WaitBar .AND. t-tBegin > iBar*tmp_tb ) CALL Progress(iBar)

          !-- If there were no failures compute a new h.
          tmp = 1.25_dp * error**ROS%pow
          IF ( TWO * tmp > ONE ) THEN
            h = h / tmp
          ELSE
            h = h * TWO
          END IF

          !-- Advance the integration one step.
          t  = tnew
          Y0 = Y

          !-- for PI stepsize control
          errorOld  = error
          hOld      = h

          !-- Termination condition for the main loop.
          IF ( done ) EXIT  MAIN_LOOP_ROSENBROCK_METHOD

        END DO MAIN_LOOP_ROSENBROCK_METHOD  ! MAIN LOOP

      CASE('LSODE')

        time_int=MPI_WTIME()

        LIW    = 20 + nDIM
        LRW    = 22 + 9*nDIM + nDIM**2
        ITOL   = 2              ! 2 if atol is an array
        RTOL1  = RtolRow        ! relative tolerance parameter
        ATOL1(1:nspc) = Atol(1)  ! atol dimension nspc+1
        ATOL1(nDim)   = Atol(2)  ! temperature tolerance

        ITASK  = 1             ! 1 for normal computation of output values of y at t = TOUT
        ISTATE = 1             ! This is the first call for a problem
        IOPT   = 1             ! optional inputs are used (see Long Desciption Part 1.)
        MF     = 22            ! Stiff method, internally generated full Jacobian.

        ALLOCATE(IWORK(LIW))
        ALLOCATE(RWORK(LRW))
        IWORK  = 0
        RWORK  = ZERO

        !RWORK(5) = 1.0d-6    ! first step size, if commented -> calculate h0
        RWORK(6) = maxStp
        RWORK(7) = minStp

        t        = Tspan(1)
        timepart = (Tspan(2)-Tspan(1)) / REAL(nOutP - 1, dp)
        tnew     = t + timepart

        MAIN_LOOP_LSODE: DO k = 1 , nOutP-1

          IWORK(6) = maxnsteps

          CALL DLSODE ( FRhs  , nDIM  , Y     , t     , tnew           &
          &           , ITOL  , RTOL1 , ATOL1 , ITASK , ISTATE , IOPT  &
          &           , RWORK , LRW   , IWORK , LIW   , dummy  , MF    )

          Out%nsteps     = Out%nsteps     + IWORK(11)
          Out%nRateEvals = Out%nRateEvals + IWORK(12)
          Out%npds       = Out%npds       + IWORK(13)
          Out%ndecomps   = Out%ndecomps   + IWORK(21)

          ! save data
          IF ( NetCdfPrint ) THEN
            TimeNetCDFA = MPI_WTIME()
            ! Time to save a step?
            IF ( done .OR. StpNetCDF < ZERO ) THEN
              IF ( ChemKin ) Temperature = Y(nDIM)
              CALL SetOutputNCDF( NetCDF, t , h , Y , Temperature )
              CALL StepNetCDF( NetCDF )
            END IF
            TimeNetCDF  = TimeNetCDF + (MPI_WTIME() - TimeNetCDFA)
          END IF

          !-- Call progress bar.
          IF ( WaitBar .AND. t-Tspan(1) >= iBar*tmp_tb ) THEN
            iBar = iBar + 1         ! iBar runs from 0 to 100
            CALL Progress(iBar)
          END IF

          t    = tnew
          tnew = tnew + timepart

        END DO MAIN_LOOP_LSODE

        CALL Progress(100) ! last * needs an extra call

      CASE('LSODES')

        time_int=MPI_WTIME()

        LIW    = 30
        LRW    = 20 + 21 * nDIM + 4*(BAT%nnz+A%nnz+2*nDIM)
        ITOL   = 2              ! 2 if atol is an array
        RTOL1  = RtolRow        ! relative tolerance parameter
        IF (Teq) THEN
          ATOL1(1:nspc) = Atol(1)  ! atol dimension nspc+1
          ATOL1(nDim)   = Atol(2)  ! temperature tolerance
        ELSE
          ATOL1(1:ns_Gas)  = Atol(1)  ! atol dimension nspc+1
          ATOL1(ns_gas+1:) = Atol(2)  ! temperature tolerance
        END IF

        ITASK  = 1             ! 1 for normal computation of output values of y at t = TOUT
        ISTATE = 1             ! This is the first call for a problem
        IOPT   = 1             ! optional inputs are used (see Long Desciption Part 1.)
        MF     = 222            ! Stiff method, internally generated full Jacobian.

        ALLOCATE(IWORK(LIW))
        ALLOCATE(RWORK(LRW))
        IWORK  = 0
        RWORK  = ZERO

        !RWORK(5) = 1.0d-6    ! first step size, if commented -> calculate h0
        RWORK(6) = maxStp
        RWORK(7) = minStp

        t        = Tspan(1)
        timepart = (Tspan(2)-Tspan(1)) / REAL(nOutP - 1, dp)
        tnew     = t + timepart

        MAIN_LOOP_LSODES: DO k = 1 , nOutP-1

          IWORK(6) = maxnsteps

          CALL DLSODES( FRhs  , nDIM  , Y     , t     , tnew           &
          &           , ITOL  , RTOL1 , ATOL1 , ITASK , ISTATE , IOPT  &
          &           , RWORK , LRW   , IWORK , LIW   , dummy  , MF    )

          Out%nsteps     = Out%nsteps     + IWORK(11)
          Out%nRateEvals = Out%nRateEvals + IWORK(12)
          Out%npds       = Out%npds       + IWORK(13)
          Out%ndecomps   = Out%ndecomps   + IWORK(21)

          ! save data
          IF ( NetCdfPrint ) THEN
            TimeNetCDFA = MPI_WTIME()
            IF ( done .OR. StpNetCDF < ZERO ) THEN
              IF ( ChemKin ) Temperature = Y(nDIM)
              CALL SetOutputNCDF( NetCDF, t , h , Y , Temperature )
              CALL StepNetCDF( NetCDF )
            END IF
            TimeNetCDF  = TimeNetCDF + (MPI_WTIME() - TimeNetCDFA)
          END IF

          !-- Call progress bar.
          IF ( WaitBar .AND. t-Tspan(1) >= iBar*tmp_tb ) THEN
            iBar = iBar + 1         ! iBar runs from 0 to 100
            CALL Progress(iBar)
          END IF

          t    = tnew
          tnew = tnew + timepart

        END DO MAIN_LOOP_LSODES

        CALL Progress(100) ! last * needs an extra call

      CASE ('bwEuler')

        h = maxStp

        BACKWARD_EULER: DO

          !-- Stretch the step if within 10% of tfinal-t.
          IF ( 1.05_dp*h >= Tspan(2)-t ) THEN
            h    = Tspan(2) - t
            done = .TRUE.
          END IF

          CALL Rosenbrock(  Y             & ! new concentration
          &               , error         & ! error value
          &               , ierr          & ! max error component
          &               , Y0            & ! current concentration
          &               , t             & ! current time
          &               , h             & ! stepsize
          &               , Euler=.TRUE. )  ! new concentration

          tnew = t + h
          IF (done) THEN
            tnew = Tspan(2)  ! Hit end point exactly.
            h    = tnew - t  ! Purify h.
          END IF

          Out%nRateEvals = Out%nRateEvals + 1
          Out%nSolves    = Out%nSolves + 1
          Out%nsteps     = Out%nsteps + 1

          IF ( NetCdfPrint ) THEN
            TimeNetCDFA = MPI_WTIME()
            IF ( done .OR. StpNetCDF < ZERO ) THEN
              IF ( ChemKin ) Temperature = Y(nDIM)
              CALL SetOutputNCDF( NetCDF, t , h , Y , Temperature )
              CALL StepNetCDF( NetCDF )
            END IF
            TimeNetCDF  = TimeNetCDF + (MPI_WTIME() - TimeNetCDFA)
          END IF

          !-- Call progress bar.
          IF ( WaitBar .AND. t-Tspan(1) >= iBar*tmp_tb ) THEN
            iBar = iBar + 1         ! iBar runs from 0 to 100
            CALL Progress(iBar)
          END IF

          !-- Termination condition for the main loop.
          IF ( done ) EXIT

          !-- Advance the integration one step.
          t   = tnew
          Y0  = Y

        END DO BACKWARD_EULER

      CASE DEFAULT

        WRITE(*,*) ' No other methods implemented jet. Use Rosenbrock, LSODE[S] or backward Euler'
        CALL DropOut()

    END SELECT
    ! save values for next initial vector
    IF ( ChemKin ) Temperature = Y(nDIM)
    y_iconc = Y
    h0 = h

    TimeIntegration = TimeIntegration + MPI_WTIME() - time_int

  END SUBROUTINE Integrate

  ! The Progressbar
  SUBROUTINE Progress(j)
    !
    INTEGER(4)    :: j,k
    CHARACTER(69) :: bar="          Start Integration...........    ???% |                    |"

    iBar = iBar + 1         ! iBar runs from 0 to 100
    !
    IF (MPI_ID==0) THEN
      WRITE(bar(43:45),'(I3)') j
      !
      DO k=1,j/5
        bar(48+k:48+k)="*"
      END DO
      ! print the progress bar.
      WRITE(*,'(A1,A69,$)') char(13), bar
    END IF
  END SUBROUTINE Progress
  !
  !
  SUBROUTINE PI_StepsizeControl(hnew,Tol,er,erOld,h,hOld,PI_Par,RCo)
    REAL(dp) :: hnew
    !
    !
    REAL(dp) :: Tol       ! rel tol
    REAL(dp) :: er        ! lokal error step n
    REAL(dp) :: erOld     ! lokal error step n-1
    REAL(dp) :: h         ! stepsize n
    REAL(dp) :: hOld      ! stepsize n-1
    TYPE(PI_Param) :: PI_Par    ! PI control parameter
    TYPE(RosenbrockMethod_T)  :: RCo
    !
    REAL(dp) :: htmp
    !
    !
    htmp = ( Tol/er )**(0.7d0*ROS%pow)* ( erOld/Tol )**(0.4d0*ROS%pow) * h
    !htmp = ( Tol/er )**PI_Par%KI * ( erOld/Tol )**PI_Par%Kp * h
    !
    ! limitation term
    IF (htmp>PI_Par%ThetaMAX*h) THEN
      hnew = PI_Par%ThetaMAX*h
    ELSE
      hnew = htmp
    END IF
  END SUBROUTINE PI_StepsizeControl

  SUBROUTINE FRhs(NEQ1, T, Y, YDOT)
    USE Rates_Mod
    INTEGER        :: NEQ1
    REAL(dp) :: T , Y(NEQ1) , YDOT(NEQ1)
    !
    REAL(dp) :: dCdt(nspc), Emiss(nspc)
    REAL(dp) :: Rate(neq) , DRate(neq)
    REAL(dp) :: Tarr(10)
    REAL(dp) :: U(nspc) , dUdT(nspc)
    REAL(dp) :: c_v

    IF (Teq) THEN         ! OUT:   IN:
      ! MASS CONSERVATION
      CALL ReactionRates( T , Y , Rate , DRate)
      CALL UpdateEmission(Emiss,T)
      dCdt = BAT * Rate + Emiss   ! [mol/cm3/s]

      YDOT(1:nspc) = dCdt

      ! ENERGY CONSERVATION
      Tarr = UpdateTempArray ( Y(NEQ1) )
      CALL InternalEnergy    (  U   , Tarr )
      CALL DiffInternalEnergy( dUdT , Tarr)
      CALL MassAveMixSpecHeat( c_v  , dUdT               &
                             &      , MoleConc=Y(1:nspc) &
                             &      , rho=rho          )

      YDOT(NEQ1) = - SUM( U * dCdt ) * rRho / c_v   ! rRho=mega/rho
    ELSE
      ! MASS CONSERVATION
      CALL ReactionRates( T , Y , Rate )
      CALL UpdateEmission(Emiss,T)
      dCdt = BAT * Rate + Emiss   ! [mol/cm3/s]

      YDOT(1:nspc) = dCdt
    END IF

  END SUBROUTINE FRhs

  SUBROUTINE Jacobian_LSODES (NEQ1, T, Y, J, IA, JA, PDJ)
    INTEGER  ::  NEQ1, J
    REAL(dp) :: T, Y(NEQ1), PDJ(NEQ1)
    INTEGER  :: IA(:), JA(:)

 !1   PDJ(1) = -RK1
     !PDJ(2) = RK1
     !RETURN
 !2   PDJ(2) = -RK3*Y(3) - RK15*Y(12) - RK2
     !PDJ(3) = RK2 - RK3*Y(3)
     !PDJ(4) = RK3*Y(3)
     !PDJ(5) = RK15*Y(12)
     !PDJ(12) = -RK15*Y(12)
     !RETURN
 !3   PDJ(2) = -RK3*Y(2)
     !PDJ(3) = -RK5 - RK3*Y(2) - RK7*Y(10)
     !PDJ(4) = RK3*Y(2)
     !PDJ(6) = RK7*Y(10)
     !PDJ(10) = RK5 - RK7*Y(10)
     !RETURN
 !4   PDJ(2) = RK11*RK14
     !PDJ(3) = RK11*RK14
     !PDJ(4) = -RK11*RK14 - RK4
     !PDJ(9) = RK4
     !RETURN
 !5   PDJ(2) = RK19*RK14
     !PDJ(5) = -RK19*RK14 - RK16
     !PDJ(9) = RK16
     !PDJ(12) = RK19*RK14
     !RETURN
 !6   PDJ(3) = RK12*RK14
     !PDJ(6) = -RK12*RK14 - RK8
     !PDJ(9) = RK8
     !PDJ(10) = RK12*RK14
     !RETURN
 !7   PDJ(7) = -RK20*RK14 - RK18
     !PDJ(9) = RK18
     !PDJ(10) = RK20*RK14
     !PDJ(12) = RK20*RK14
     !RETURN
 !8   PDJ(8) = -RK13*RK14 - RK10
     !PDJ(10) = RK13*RK14
     !PDJ(11) = RK10
 !9   RETURN
 !10  PDJ(3) = -RK7*Y(3)
     !PDJ(6) = RK7*Y(3)
     !PDJ(7) = RK17*Y(12)
     !PDJ(8) = RK9
     !PDJ(10) = -RK7*Y(3) - RK17*Y(12) - RK6 - RK9
     !PDJ(12) = RK6 - RK17*Y(12)
 !11  RETURN
 !12  PDJ(2) = -RK15*Y(2)
     !PDJ(5) = RK15*Y(2)
     !PDJ(7) = RK17*Y(10)
     !PDJ(10) = -RK17*Y(10)
     !PDJ(12) = -RK15*Y(2) - RK17*Y(10)
     !RETURN
  END SUBROUTINE Jacobian_LSODES

  ! dummy routine for ode solver LSODE
  SUBROUTINE dummy()
    IMPLICIT NONE
  END SUBROUTINE dummy

END MODULE Integration_Mod

max_stars_repo_name: KarlJohnsonnn/AtCSol

max_issues_repo_name: KarlJohnsonnn/AtCSol

max_forks_repo_name: KarlJohnsonnn/AtCSol


text:
!==========================================================================
elemental function gsw_rho_ice (t, p)
!==========================================================================
!
!  Calculates in-situ density of ice from in-situ temperature and pressure.
!  Note that the output, rho_ice, is density, not density anomaly;  that
!  is, 1000 kg/m^3 is not subracted from it.
!
!  t   =  in-situ temperature (ITS-90)                            [ deg C ]
!  p   =  sea pressure                                             [ dbar ]
!         ( i.e. absolute pressure - 10.1325 dbar )
!
!  rho_ice  =  in-situ density of ice (not density anomaly)      [ kg/m^3 ]
!--------------------------------------------------------------------------

use gsw_mod_toolbox, only : gsw_gibbs_ice

use gsw_mod_kinds

implicit none

real (r8), intent(in) :: t, p

real (r8) :: gsw_rho_ice

gsw_rho_ice = 1.0_r8/gsw_gibbs_ice(0,1,t,p)

return
end function

!--------------------------------------------------------------------------

max_stars_repo_name: minsukji/ci-debug

max_issues_repo_name: minsukji/ci-debug

max_forks_repo_name: minsukji/ci-debug


text:
   subroutine loadcurve_label2num
!
! Convert loadcurve labels to their internal number for options that
! do not call an initialization routine (usually in initlz).
!
!     use mod_assembly
!
     implicit none
!
   integer, parameter:: singI = kind(0)
   type fric_damp_data_T
   integer   :: lc_damp_dev   = 0
   end type fric_damp_data_T
!
   type part_T
   type(fric_damp_data_T),pointer :: fric_damp_data => Null()
   end type part_T
!
   type assembly_T
   type(part_T), dimension(:),pointer,contiguous :: part => Null()
   end type assembly_T
!
   type(assembly_T),save,target :: assembly
   integer(singI) :: lc_num
     type(part_T), pointer :: part
!
   type(fric_damp_data_T), pointer :: tmp
!
     part => assembly%part(1)
       tmp = part%fric_damp_data
!     lc_num = part%fric_damp_data%lc_damp_dev
!
   end subroutine loadcurve_label2num

max_stars_repo_name: ouankou/rose

max_issues_repo_name: WildeGeist/rose

max_forks_repo_name: WildeGeist/rose


PROGRAM FIZZBUZZ

  IMPLICIT NONE

  INTEGER :: i, max_num

  max_num = 100  ! Classic FizzBuzz runs from 1 to 100

  DO i = 1, max_num

    ! Start by considering divides by 3 AND 5
    ! The MOD function returns the remainder of division
    ! SO MOD(4, 3) -> 1 and MOD(7, 5) -> 2 etc
    ! For the first case, we could also do:
    ! IF(MOD(i, 3) == 0 .AND. MOD(i, 5) == 0) THEN
    IF(MOD(i, 15) == 0) THEN
      ! Divides by 15 i.e. by 5 and by 3
      PRINT*, "Fizzbuzz"
    ELSE IF(MOD(i, 3) == 0) THEN
      !Divdes by 3 but NOT 5
      PRINT*, "Fizz"
    ELSE IF(MOD(i, 5) == 0) THEN
      !Divdes by 5 but NOT 3
      PRINT*, "Buzz"
    ELSE
      ! Format as an integer with up to 8 digits
      PRINT'(I8)', i
    END IF

  ENDDO

END PROGRAM


text:
c*************************************************************************
c                            RMVS_STEP.F
c*************************************************************************
c This subroutine takes a step in helio coord.
c both massive and test particles
c INCLUDES close approuches between test particles and planets
c
c             Input:
c                 i1st          ==>  = 0 if first step; = 1 not (int scalar)
c                 time          ==>  current time (real scalar)
c                 nbod          ==>  number of massive bodies (int scalar)
c                 ntp            ==>  number of massive bodies (int scalar)
c                 mass          ==>  mass of bodies (real array)
c                 j2rp2,j4rp4   ==>  J2*radii_pl^2 and  J4*radii_pl^4
c                                     (real scalars)
c                 xh,yh,zh      ==>  initial planet position in helio coord
c                                    (real arrays)
c                 vxh,vyh,vzh   ==>  initial planet velocity in helio coord
c                                    (real arrays)
c                 xht,yht,zht    ==>  initial tp position in helio coord
c                                      (real arrays)
c                 vxht,vyht,vzht ==>  initial tp velocity in helio coord
c                                        (real arrays)
c                 istat           ==>  status of the test paricles
c                                      (2d integer array)
c                                      istat(i,1) = 0 ==> active:  = 1 not
c                                      istat(i,2) = -1 ==> Danby did not work
c                 rstat           ==>  status of the test paricles
c                                      (2d real array)
c                 dt            ==>  time step
c             Output:
c                 xh,yh,zh      ==>  final planet position in helio coord
c                                       (real arrays)
c                 vxh,vyh,vzh   ==>  final planet velocity in helio coord
c                                       (real arrays)
c                 xht,yht,zht    ==>  final tp position in helio coord
c                                       (real arrays)
c                 vxht,vyht,vzht ==>  final tp position in helio coord
c                                       (real arrays)
c
c
c Remarks: Adopted from martin's nbwh.f program
c Authors:  Hal Levison
c Date:    2/19/93
c Last revision: 2/24/94

      subroutine rmvs_step(i1st,time,nbod,ntp,mass,j2rp2,j4rp4,
     &     xh,yh,zh,vxh,vyh,vzh,xht,yht,zht,vxht,vyht,vzht,
     &     istat,rstat,dt)

      include '../swift.inc'
      include 'rmvs.inc'

c...  Inputs Only:
      integer nbod,ntp,i1st
      real*8 mass(nbod),dt,time,j2rp2,j4rp4

c...  Inputs and Outputs:
      integer istat(NTPMAX,NSTAT)
      real*8 rstat(NTPMAX,NSTATR)
      real*8 xh(nbod),yh(nbod),zh(nbod)
      real*8 vxh(nbod),vyh(nbod),vzh(nbod)
      real*8 xht(ntp),yht(ntp),zht(ntp)
      real*8 vxht(ntp),vyht(ntp),vzht(ntp)

c...  Internals
      integer i1sttp,i1sto,icflg,i,j,np
      real*8 xtmp(NPLMAX,NTPENC),ytmp(NPLMAX,NTPENC)
      real*8 ztmp(NPLMAX,NTPENC),peri(NTPMAX)
      real*8 xbeg(NPLMAX),ybeg(NPLMAX),zbeg(NPLMAX)
      real*8 vxbeg(NPLMAX),vybeg(NPLMAX),vzbeg(NPLMAX)
      real*8 xend(NPLMAX),yend(NPLMAX),zend(NPLMAX)
      real*8 vxend(NPLMAX),vyend(NPLMAX),vzend(NPLMAX)
      Integer nenco(NPLMAX),itpenco(NTPMAX,NPLMAX),nenci(NPLMAX)
      integer itpenci(NTPMAX,NPLMAX),ienci(NTPMAX),ienco(NTPMAX)
      integer istattmp(NTPMAX,NSTAT),isperi(NTPMAX)

c----
c...  Executable code

      i1sttp = i1st
      i1sto = i1st

c...  are there any encounters?
      call rmvs_chk(nbod,ntp,mass,xh,yh,zh,vxh,vyh,vzh,xht,yht,zht,
     &       vxht,vyht,vzht,istat,dt,icflg,nenco,itpenco,nenci,
     &       itpenci,ienci,ienco)

c...  keep track of encounters in the inner region
      call rmvs_elog(ntp,icflg,ienci,istat)

c.... if not just do a normal step and leave
      if(icflg.eq.0) then
         call step_kdk(i1st,time,nbod,ntp,mass,j2rp2,j4rp4,xh,yh,zh,
     &        vxh,vyh,vzh,xht,yht,zht,vxht,vyht,vzht,istat,rstat,dt)

         do i=1,ntp
            if(istat(i,1).eq.0) then
               istat(i,2) = 0
            endif
         enddo

         return       !  NOTE AN EXIT
      endif

c...  ENCOUNTER STUFF FROM HERE ON!!!!!

c...  save initial x and v of planets if there are planocentric enc
      do i=1,nbod
            xbeg(i) = xh(i)
            ybeg(i) = yh(i)
            zbeg(i) = zh(i)
            vxbeg(i) = vxh(i)
            vybeg(i) = vyh(i)
            vzbeg(i) = vzh(i)
      enddo

c...  first do the outer region integration, save the planet position
      call rmvs_step_out(i1st,nbod,ntp,mass,j2rp2,j4rp4,xh,yh,zh,
     &              vxh,vyh,vzh,xht,yht,zht,vxht,vyht,vzht,istat,
     &              icflg,ienco,dt,xtmp,ytmp,ztmp)

c...  save the final position and velocity of planets
      do i=1,nbod
         xend(i) = xh(i)
         yend(i) = yh(i)
         zend(i) = zh(i)
         vxend(i) = vxh(i)
         vyend(i) = vyh(i)
         vzend(i) = vzh(i)
      enddo

c...  do the plaocentric stuff if necessary
      if(icflg.eq.-1) then
         call rmvs_step_in(i1sttp,nbod,ntp,mass,j2rp2,j4rp4,xtmp,
     &           ytmp,ztmp,xbeg,ybeg,zbeg,vxbeg,vybeg,vzbeg,
     &           xend,yend,zend,vxend,vyend,vzend,xht,yht,zht,
     &           vxht,vyht,vzht,istat,nenci,itpenci,isperi,peri,dt)
         do i=1,ntp
            if(istat(i,1).eq.0) then
               istat(i,2) = 0
            endif
         enddo
      endif

c...  As of this point all the test particles that are involved in an
c...  encounter have been moved.  But not the ones that have not.
c...  so move those,  BUT NOT the onces in the encounter

c...  make a temporary istat array ao only they are active
      do i=1,ntp
         if(istat(i,1).eq.0) then
            if( (ienco(i).ne.0) .or. (ienci(i).ne.0) )  then
               istattmp(i,1) = 1 ! don't integrate it
            else
               istattmp(i,1) = 0 ! integrate it
            endif
            do j=2,NSTAT
               istattmp(i,j) = 0
            enddo
         else
            istattmp(i,1) = 1   ! don't integrate it
         endif
       enddo

c...  do a full step
      i1sto = 0      ! we need to recalculate accel arrays
      call step_kdk_tp(i1sto,nbod,ntp,mass,j2rp2,j4rp4,
     &              xbeg,ybeg,zbeg,xend,yend,zend,
     &              xht,yht,zht,vxht,vyht,vzht,istattmp,dt)

c...  fix up the istat array
      do i=1,ntp
         if(istattmp(i,2) .ne. 0) then   ! danby screwed up
            istat(i,1) = 1
            do j=2,NSTAT
               istat(i,j) = istattmp(i,j)
            enddo
          endif
       enddo

c...  put the enc info into istat
      do i=1,ntp
         if(istat(i,1).eq.0) then
            if(ienci(i).ne.0) then
               istat(i,2) = -ienci(i)
               istat(i,3) = ienci(i)
               if(isperi(i).eq.0) then
                  rstat(i,3) = peri(i)
                  np = NSTATP + ienci(i) - 1
                  if(rstat(i,np).eq.0) then
                     rstat(i,np) = peri(i)
                  else
                     rstat(i,np) = min(rstat(i,np),peri(i))
                  endif
               else
                  rstat(i,3) = 0.0d0
               endif
            else if(ienco(i).ne.0) then
               istat(i,2) = ienco(i)
            else
               istat(i,2) = 0
            endif
         endif
      enddo

c...  we MUST make sure that the saved accel arrays are ok
c...  calculate them again
       i1st = 0

      return
      end   ! step_enc
c------------------------------------------------------------------------

max_stars_repo_name: Simske/exostriker


text:
module math_utilities
!*Brief Description:* This module contains solvers required for lung problems
!
!*LICENSE:*
!
!
!
!*Full Description:*
!
!
  use arrays, only: dp
  use diagnostics, only: enter_exit

  implicit none
  private
  public ax_cr,diagonal_pointer_cr,ilu_cr,lus_cr,mult_givens,rearrange_cr
  public sort_integer_list
  public sort_real_list

contains
!
!###########################################################################
!
!*ax_cr:* Computes A*x for a matrix stored in sparse compressed row form
  subroutine ax_cr ( n, ia, ja, a, x, w )
    implicit none

    integer ( kind = 4 ) n !the order of the system
    integer ( kind = 4 ) ia(*) !ia(n+1) row indices
    integer ( kind = 4 ) ja(*) !ja(nz_num) column indices
    real ( kind = 8 ) a(*) !a(nz_num) Matrix values
    real ( kind = 8 ) x(*) !x(n) Vector to be multiplied by A
    real ( kind = 8 ) w(*) !w(n) Value of A*x

    integer ( kind = 4 ) i
    integer ( kind = 4 ) k1
    integer ( kind = 4 ) k2

    w(1:n) = 0.0D+00

    do i = 1, n
       k1 = ia(i)
       k2 = ia(i+1) - 1
       w(i) = w(i) + dot_product ( a(k1:k2), x(ja(k1:k2)) )
    end do

    return
  end subroutine ax_cr
!
!##############################################################################
!
! *ILU_CR:* computes the incomplete LU factorization of a matrix. For a matrix
! stored in compressed row format.
    !    Input, integer ( kind = 4 ) UA(N), the index of the diagonal element
    !    of each row.
    !    Output, real ( kind = 8 ) L(NZ_NUM), the ILU factorization of A.
  subroutine ilu_cr ( n, nz_num, ia, ja, a, ua, l )
    integer ( kind = 4 ) n
    integer ( kind = 4 ) nz_num
    integer ( kind = 4 ) ia(*) !ia(n+1)
    integer ( kind = 4 ) ja(*) !ja(nz_num)
    real ( kind = 8 ) a(*) !a(nz_num)
    integer ( kind = 4 ) ua(*) !ua(n)
    real ( kind = 8 ) l(*) !l(nz_num)

    integer ( kind = 4 ) i
    integer ( kind = 4 ) iw(n)
    integer ( kind = 4 ) j
    integer ( kind = 4 ) jj
    integer ( kind = 4 ) jrow
    integer ( kind = 4 ) jw
    integer ( kind = 4 ) k
    real ( kind = 8 ) tl

    !  Copy A.
    l(1:nz_num) = a(1:nz_num)

    do i = 1, n ! for each row, up to max number of rows
       !  IW points to the nonzero entries in row I.
       iw(1:n) = -1
       do k = ia(i), ia(i+1) - 1 !for each
          iw(ja(k)) = k
       end do
       do j = ia(i), ia(i+1) - 1
          jrow = ja(j)
          if ( i <= jrow ) then
             exit
          end if
          tl = l(j) * l(ua(jrow))
          l(j) = tl
          do jj = ua(jrow) + 1, ia(jrow+1) - 1
             jw = iw(ja(jj))
             if ( jw /= -1 ) then
                l(jw) = l(jw) - tl * l(jj)
             end if
          end do
       end do
       ua(i) = j
       if ( jrow /= i ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'ILU_CR - Fatal error!'
          write ( *, '(a)' ) '  JROW ~= I'
          write ( *, '(a,i8)' ) '  JROW = ', jrow
          write ( *, '(a,i8)' ) '  I    = ', i
          stop
       end if
       if ( l(j) == 0.0D+00 ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'ILU_CR - Fatal error!'
          write ( *, '(a,i8)' ) '  Zero pivot on step I = ', i
          write ( *, '(a,i8,a)' ) '  L(', j, ') = 0.0'
          stop
       end if
       l(j) = 1.0D+00 / l(j)
    end do

    l(ua(1:n)) = 1.0D+00 / l(ua(1:n))

    return
  end subroutine ilu_cr
!
!##############################################################################
!
!*DIAGONAL_POINTER_CR:* finds diagonal entries in a sparse compressed row matrix.
    !    The array UA can be used to locate the diagonal elements of the matrix.
    !    It is assumed that every row of the matrix includes a diagonal element,
    !    and that the elements of each row have been ascending sorted.
subroutine diagonal_pointer_cr ( n, ia, ja, ua )
    integer ( kind = 4 ) n
    integer ( kind = 4 ) ia(*) !ia(n+1)
    integer ( kind = 4 ) ja(*) !ja(nz_num)
    integer ( kind = 4 ) ua(*) !ua(n)

    integer ( kind = 4 ) i
    integer ( kind = 4 ) k

    ua(1:n) = -1

    do i = 1, n
       do k = ia(i), ia(i+1) - 1
          if ( ja(k) == i ) then
             ua(i) = k
          end if
       end do
    end do
    return
  end subroutine diagonal_pointer_cr

  !*****************************************************************************80

  subroutine lus_cr ( n, ia, ja, l, ua, r, z )
!!! LUS_CR applies the incomplete LU preconditioner.
    !    The linear system M * Z = R is solved for Z.  M is the incomplete
    !    LU preconditioner matrix, and R is a vector supplied by the user.
    !    So essentially, we're solving L * U * Z = R.
    !    Input, integer ( kind = 4 ) UA(N), the index of the diagonal element
    !    of each row.
    !    Input, real ( kind = 8 ) R(N), the right hand side.
    !    Output, real ( kind = 8 ) Z(N), the solution of the system M * Z = R.
    implicit none

    integer ( kind = 4 ) n
    integer ( kind = 4 ) ia(*) !ia(n+1)
    integer ( kind = 4 ) ja(*) !ja(nz_num)
    real ( kind = 8 ) l(*) !l(nz_num)
    integer ( kind = 4 ) ua(*) !ua(n)
    real ( kind = 8 ) r(*) !r(n)

    integer ( kind = 4 ) i
    integer ( kind = 4 ) j
    real ( kind = 8 ) w(n)
    real ( kind = 8 ) z(n)

    !  Copy R in.
    w(1:n) = r(1:n)

    !  Solve L * w = w where L is unit lower triangular.
    do i = 2, n
       do j = ia(i), ua(i) - 1
          w(i) = w(i) - l(j) * w(ja(j))
       end do
    end do

    !  Solve U * w = w, where U is upper triangular.
    do i = n, 1, -1
       do j = ua(i) + 1, ia(i+1) - 1
          w(i) = w(i) - l(j) * w(ja(j))
       end do
       w(i) = w(i) / l(ua(i))
    end do

    !  Copy Z out.
    z(1:n) = w(1:n)

    return
  end subroutine lus_cr

  !*****************************************************************************80
  subroutine mult_givens ( c, s, k, g )
!!! MULT_GIVENS applies a Givens rotation to two successive entries of a vector.
    !    In order to make it easier to compare this code with the Original C,
    !    the vector indexing is 0-based.
    !    Input, real ( kind = 8 ) C, S, the cosine and sine of a Givens
    !    rotation.
    !
    !    Input, integer ( kind = 4 ) K, indicates the location of the first
    !    vector entry.
    !
    !    Input/output, real ( kind = 8 ) G(1:K+1), the vector to be modified.
    !    On output, the Givens rotation has been applied to entries G(K) and G(K+1).

    implicit none

    real ( kind = 8 ) c
    real ( kind = 8 ) s
    integer ( kind = 4 ) k
    real ( kind = 8 ) g(*) !g(1:k+1)

    real ( kind = 8 ) g1
    real ( kind = 8 ) g2

    g1 = c * g(k) - s * g(k+1)
    g2 = s * g(k) + c * g(k+1)

    g(k)   = g1
    g(k+1) = g2

    return
  end subroutine mult_givens

  !*****************************************************************************80
  subroutine rearrange_cr ( n, ia, ja, a )
!!! REARRANGE_CR sorts a sparse compressed row matrix.
    !    This routine guarantees that the entries in the CR matrix
    !    are properly sorted.
    !
    !    After the sorting, the entries of the matrix are rearranged in such
    !    a way that the entries of each column are listed in ascending order
    !    of their column values.
    !    Input, integer ( kind = 4 ) N, the order of the system.
    !
    !    Input, integer ( kind = 4 ) NZ_NUM, the number of nonzeros.
    !
    !    Input, integer ( kind = 4 ) IA(N+1), the compressed row indices.
    !
    !    Input/output, integer ( kind = 4 ) JA(NZ_NUM), the column indices.
    !    On output, these may have been rearranged by the sorting.
    !
    !    Input/output, real ( kind = 8 ) A(NZ_NUM), the matrix values.  On output,
    !    the matrix values may have been moved somewhat because of the sorting.
    !
    implicit none

    integer ( kind = 4 ) n
    integer ( kind = 4 ) ia(*) !ia(n+1)
    integer ( kind = 4 ) ja(*) !ja(nz_num)
    real ( kind = 8 ) a(*) !a(nz_num)

    integer ( kind = 4 ) i
    integer ( kind = 4 ) i4temp
    integer ( kind = 4 ) k
    integer ( kind = 4 ) l
    real ( kind = 8 ) r8temp

    do i = 1, n

       do k = ia(i), ia(i+1) - 2
          do l = k + 1, ia(i+1) - 1

             if ( ja(l) < ja(k) ) then
                i4temp = ja(l)
                ja(l)  = ja(k)
                ja(k)  = i4temp

                r8temp = a(l)
                a(l)   = a(k)
                a(k)   = r8temp
             end if

          end do
       end do

    end do

    return
  end subroutine rearrange_cr

!######################################################################

  !
  !*sort_integer_list:* sorts a list of integer values into a non-decreasing order.
  ! sorts N integer IDATA values into a non-decreasing sequence using IHEAPSORT
  ! (N>50) or ISHELLSORT (N>50) and then  removes all duplicates from the list. On
  ! exit N contains the number of unique elements in the list.
  !
  subroutine sort_integer_list(N,IDATA)
    !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_SORT_INTEGER_LIST" :: SORT_INTEGER_LIST
    integer :: IDATA(:),N

    !Local Variables
    integer :: count1,index,itemp,nolist
    logical :: continue

    character(len=60) :: sub_name

    sub_name = 'sort_integer_list'
    call enter_exit(sub_name,1)

    !order the array non-decreasing
    do nolist=2,N
       count1=0
       continue=.true.
       do while(continue)
          if(IDATA(nolist-count1).lt.IDATA(nolist-count1-1))then
             itemp=IDATA(nolist-1)
             IDATA(nolist-1)=IDATA(nolist)
             IDATA(nolist)=itemp
             count1=count1+1
             if(nolist-count1-1.eq.0) continue=.false.
          else
             continue=.false.
          endif
       enddo !while
    enddo !N

    !eliminate duplicate entries
    index=0
    do nolist=2,N
       if(IDATA(nolist).eq.IDATA(nolist-1)) then
          index=index+1
       else
          IDATA(nolist-index)=IDATA(nolist)
       endif
    enddo !nolist

    N=N-index

    call enter_exit(sub_name,2)

  end subroutine sort_integer_list

!!!#########################################################################
  !*sort_real_list:* sorts a list of real values into a non-decreasing order
  ! using a bubble sort algorithm.

  subroutine sort_real_list(n,RDATA,INDEX)
    !DEC$ ATTRIBUTES DLLEXPORT,ALIAS:"SO_SORT_REAL_LIST" :: SORT_REAL_LIST

    integer :: INDEX(*),n
    real(dp) :: RDATA(*)

    !Local Variables
    integer :: FLAG,i,ITEMP,j,k
    real(dp) :: TEMP

    character(len=60) :: sub_name

    sub_name = 'sort_real_list'
    call enter_exit(sub_name,1)

    if(N.LE.1) then
    else
       FLAG=n
       do i=1,n
          k=FLAG-1
          FLAG=0
          do j=1,k
             if(RDATA(j).gt.RDATA(j+1)) then
                TEMP=RDATA(j)
                RDATA(j)=RDATA(j+1)
                RDATA(j+1)=TEMP
                ITEMP=INDEX(j)
                INDEX(j)=INDEX(j+1)
                INDEX(j+1)=ITEMP
                FLAG=j
             endif
          enddo
          if(FLAG.eq.0) then
             write(*,*) 'warning in rsort'
          endif
       enddo
    endif

    call enter_exit(sub_name,2)

  end subroutine sort_real_list

end module math_utilities

max_stars_repo_name: AnnikaGspandl/lungsim

max_issues_repo_name: AnnikaGspandl/lungsim

max_forks_repo_name: AnnikaGspandl/lungsim


! { dg-do run }
! { dg-options "-std=gnu" }
  integer(kind=4) :: i4, j4
  integer(kind=8) :: i8, j8

  i4 = mclock()
  i8 = mclock8()
  j4 = mclock()
  j8 = mclock8()

  if (i4 > j4 .or. i8 > j8 .or. i4 > i8 .or. j4 > j8) STOP 1

  end


text:
      SUBROUTINE CODGET
C     GET THE CODE OF THE NEXT CHARACTER IN THE INPUT
C
      INCLUDE 'PARS.f'
      INCLUDE 'CODS.f'
      INCLUDE 'LINS.f'
C
      CALL ZZCI (LINE,FLIN,SC)
      FLIN = FLIN + 1
      SC = SC + 1
      SC = CODTAB(SC)
      RETURN
      END

max_stars_repo_name: gfis/jextra

max_issues_repo_name: gfis/jextra

max_forks_repo_name: gfis/jextra


!========================================================================
!
!                   S P E C F E M 2 D  Version 7 . 0
!                   --------------------------------
!
!     Main historical authors: Dimitri Komatitsch and Jeroen Tromp
!                        Princeton University, USA
!                and CNRS / University of Marseille, France
!                 (there are currently many more authors!)
! (c) Princeton University and CNRS / University of Marseille, April 2014
!
! This software is a computer program whose purpose is to solve
! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
! using a spectral-element method (SEM).
!
! This software is governed by the CeCILL license under French law and
! abiding by the rules of distribution of free software. You can use,
! modify and/or redistribute the software under the terms of the CeCILL
! license as circulated by CEA, CNRS and Inria at the following URL
! "http://www.cecill.info".
!
! As a counterpart to the access to the source code and rights to copy,
! modify and redistribute granted by the license, users are provided only
! with a limited warranty and the software's author, the holder of the
! economic rights, and the successive licensors have only limited
! liability.
!
! In this respect, the user's attention is drawn to the risks associated
! with loading, using, modifying and/or developing or reproducing the
! software by the user in light of its specific status of free software,
! that may mean that it is complicated to manipulate, and that also
! therefore means that it is reserved for developers and experienced
! professionals having in-depth computer knowledge. Users are therefore
! encouraged to load and test the software's suitability as regards their
! requirements in conditions enabling the security of their systems and/or
! data to be ensured and, more generally, to use and operate it in the
! same conditions as regards security.
!
! The full text of the license is available in file "LICENSE".
!
!========================================================================

  subroutine define_derivation_matrices()

  use specfem_par, only: xigll,zigll,wxgll,wzgll,hprime_xx,hprime_zz,hprimewgll_xx,hprimewgll_zz
  implicit none

  include "constants.h"

  double precision, parameter :: alphaGLL = 0.d0, betaGLL = 0.d0

! function for calculating derivatives of Lagrange polynomials
  double precision, external :: lagrange_deriv_GLL

  integer i1,i2,k1,k2

! set up coordinates of the Gauss-Lobatto-Legendre points
  call zwgljd(xigll,wxgll,NGLLX,alphaGLL,betaGLL)
  call zwgljd(zigll,wzgll,NGLLZ,alphaGLL,betaGLL)

! if number of points is odd, the middle abscissa is exactly zero
  if(mod(NGLLX,2) /= 0) xigll((NGLLX-1)/2+1) = ZERO
  if(mod(NGLLZ,2) /= 0) zigll((NGLLZ-1)/2+1) = ZERO

! calculate derivatives of the Lagrange polynomials
! and precalculate some products in double precision
! hprime(i,j) = h'_j(xigll_i) by definition of the derivation matrix
  do i1=1,NGLLX
    do i2=1,NGLLX
      hprime_xx(i2,i1) = lagrange_deriv_GLL(i1-1,i2-1,xigll,NGLLX)
      hprimewgll_xx(i2,i1) = wxgll(i2) * hprime_xx(i2,i1)
    enddo
  enddo

  do k1=1,NGLLZ
    do k2=1,NGLLZ
      hprime_zz(k2,k1) = lagrange_deriv_GLL(k1-1,k2-1,zigll,NGLLZ)
      hprimewgll_zz(k2,k1) = wzgll(k2) * hprime_zz(k2,k1)
    enddo
  enddo

  end subroutine define_derivation_matrices

!____________________________________________________________________________________
!
!subroutine define_GLJ_derivation_matrix(xiglj,wxglj,hprimeBar_xx,hprimeBarwglj_xx)
! Calculate all that we need for the GLJ quadrature on axial elements :
! Weights, GLJ points and derivatives of polynomials at GLJ points.
!____________________________________________________________________________________
!

  subroutine define_GLJ_derivation_matrix()

  use specfem_par, only: xiglj,wxglj,hprimeBar_xx,hprimeBarwglj_xx
  implicit none

  include "constants.h"

  double precision, parameter    :: alphaGLJ=0.d0,betaGLJ=1.d0

! function for calculating derivatives of GLJ polynomials
  double precision, external :: poly_deriv_GLJ

  integer i1,i2

! set up coordinates of the Gauss-Lobatto-Jacobi points
  call zwgljd(xiglj,wxglj,NGLJ,alphaGLJ,betaGLJ)

! calculate derivatives of the GLJ quadrature polynomials
! and precalculate some products in double precision
! hprimeBar(i,j) = hBar'_j(xiglj_i) by definition of the derivation matrix
  do i1=1,NGLJ
    do i2=1,NGLJ
      hprimeBar_xx(i2,i1) = poly_deriv_GLJ(i1-1,i2-1,xiglj,NGLJ)
      hprimeBarwglj_xx(i2,i1) = wxglj(i2) * hprimeBar_xx(i2,i1)
    enddo
  enddo

  end subroutine define_GLJ_derivation_matrix


program rk_4
  use tools

  implicit none
  real(dp) :: x0, x_point, y0, h

  x0 = 0.0; y0 = 1.0; h = 0.05
  x_point = 0.25 !the value of x at which y is to be calculated

  write (*, form_header3) 'itn', 'x', 'y'
  write (*, line3) '-------', '-------', '-------'

  call rungkut_4(x0, y0, x_point, h)

  print *, ''
  print *, 'result:'
  write (*, form_result2r) x0, '-- given value of x'
  write (*, form_result2r) y0, '-- the required value of y'
end program rk_4


      subroutine calsoft_hyd_bfr

      use hru_module, only : cn2, hru, hru_init
      use soil_module
      use plant_module
      use hydrograph_module
      use ru_module
      use aquifer_module
      use channel_module
      use hru_lte_module
      use sd_channel_module
      use basin_module
      use maximum_data_module
      use calibration_data_module
      use conditional_module
      use reservoir_module
      use organic_mineral_mass_module
      use time_module

      implicit none

      integer :: iter_all      !none      |counter
      integer :: iterall       !none      |counter

      ! calibrate hydrology
        iter_all = 1

      do iterall = 1, iter_all

        ! calibrate cn3_swf for surface runoff
        call calsoft_hyd_bfr_et

        ! calibrate cn3_swf for surface runoff
        call calsoft_hyd_bfr_surq

        ! calibrate latq_co for water yield
        call calsoft_hyd_bfr_latq

        ! calibrate perco for water yield
        call calsoft_hyd_bfr_perc

        ! calibrate cn3_swf for surface runoff
        call calsoft_hyd_bfr_surq

      end do    ! iter_all loop

      !cal_codes%hyd_hru = "n"

	  return
      end subroutine calsoft_hyd_bfr


text:
      subroutine setprob
      implicit double precision (a-h,o-z)
      common /comaux/ z1,c1,z2,c2
c
c     # Set the material parameters for the acoustic equations
c
      !open(unit=7,file='setprob.data',status='old',form='formatted')
      call opendatafile(7, 'setprob.data')
c
c     # Piecewise constant medium with single interface as specified
c     # in setaux.f

c     # Impedance and sound speed in the two materials:

      read(7,*) z1
      read(7,*) c1
      read(7,*) z2
      read(7,*) c2

      return
      end

max_stars_repo_name: ClimateImpactLab/classic

max_issues_repo_name: ClimateImpactLab/classic

max_forks_repo_name: cjvogl/classic


text:
!===============================================================================
!> Add gravity to any solver
!===============================================================================
MODULE gravity_mod
  USE patch_mod
  USE units_mod
  USE scaling_mod
  USE link_mod
  USE io_mod
  USE kinds_mod
  USE trace_mod
  implicit none
  private
  type, public:: gravity_t
    procedure(const_grav), pointer :: grav
    real(kind=KindScalarVar)       :: const_grav = 0.0
    real(kind=KindScalarVar)       :: mass = -1.0
    real(8)                        :: minrad = 0.0
    real(8)                        :: position(3) = 0.0                    ! centre of gravity
    integer                        :: axis = 3
    class(patch_t), pointer:: patch
    logical:: on=.false.
  contains
    procedure:: init
    procedure:: pre_update
  end type
  type(gravity_t), public:: gravity
CONTAINS

!===============================================================================
!> Initialize gravity
!===============================================================================
SUBROUTINE init (self, link)
  class(gravity_t):: self
  class(link_t), pointer :: link
  ! ............................................................................
  class(patch_t), pointer:: patch
  integer                           :: iostat, gn(3)
  real(kind=KindScalarVar), save    :: constant = 0.0
  logical, save                     :: first_time=.true.
  logical, save :: on          = .false.
  real(kind=KindScalarVar), save:: mass = -1.0
  real(8),                  save:: minrad = 0.0
  real(8), save :: position(3) = 0.0                    ! centre of gravity
  integer, save :: axis  = 3
  namelist /gravity_params/ on, constant, axis, position, mass, minrad
  !-----------------------------------------------------------------------------
  call trace%begin ('selfgravity_t%init')
  !$omp critical (input_cr)
  if (first_time) then
    first_time = .false.
    rewind (io%input)
    read (io%input, gravity_params, iostat=iostat)
    if (io%master) write (io%output, gravity_params)
    if (constant/=0.0) constant = constant * scaling%t**2 / scaling%l
  end if
  !$omp end critical (input_cr)
  if (on) then
    self%patch => task2patch(link%task)
    self%on = on
    self%const_grav = constant
    self%axis = axis
    self%position   = position
    self%mass       = mass
    self%minrad     = minrad
    if (constant /= 0.0) then
      self%grav => const_grav
    else if (mass > -1.0) then
      self%grav => var_grav
    else
      call io%abort("gravity_mod%init:: neither constant nor mass is set.")
    end if
    if (.not.allocated(self%patch%force_per_unit_mass)) then
      gn = self%patch%gn
      allocate (self%patch%force_per_unit_mass(gn(1),gn(2),gn(3),3))
      self%patch%force_per_unit_mass = 0.0
    end if
  end if
  call trace%end()
END SUBROUTINE init

!===============================================================================
!===============================================================================
SUBROUTINE pre_update(self)
  class(gravity_t) :: self
  ! ----------------------------------------------------------------------------
  call self%grav()
END SUBROUTINE pre_update

!===============================================================================
!===============================================================================
SUBROUTINE const_grav(self)
  class(gravity_t) :: self
  ! ----------------------------------------------------------------------------
  self%patch%force_per_unit_mass(:,:,:,self%axis) = &
  self%patch%force_per_unit_mass(:,:,:,self%axis) + self%const_grav
END SUBROUTINE const_grav

!===============================================================================
!===============================================================================
SUBROUTINE var_grav(self)
  class(gravity_t) :: self
  ! ----------------------------------------------------------------------------
  integer:: ix, iy, iz, i
  real(8):: rp
  real(8):: x, y, z, xs, ys, zs, xp, yp, zp, ds2(3)
  real(kind=KindScalarVar):: newton
  ! ............................................................................
  newton = self%mass*scaling%grav
  associate (m=>self%patch%mesh)
  associate (r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
  do i=1,3
    ds2(i) = m(i)%h(self%patch%idx%px+i-1)  *m(i)%d
  end do
  do iz=m(3)%lb,m(3)%ub
    z = m(3)%p + r3(iz)
    zp = z  - self%position(3)
    zs = zp + ds2(3)
    do iy=m(2)%lb,m(2)%ub
      y = m(2)%p + r2(iy)
      yp = y  - self%position(2)
      ys = yp + ds2(2)
      do ix=m(1)%lb,m(1)%ub
        x = m(1)%p + r1(ix)
        xp = x  - self%position(1)
        xs = xp + ds2(1)
        !-----------------------------------------------------------------------
        ! The force of gravity is included explicitly, and is oriented
        ! away from the center.
        !-----------------------------------------------------------------------
        rp = sqrt(xs**2+yp**2+zp**2)
        rp = max(rp,self%minrad)
        self%patch%force_per_unit_mass(ix,iy,iz,1) = &
        self%patch%force_per_unit_mass(ix,iy,iz,1) - newton*xs/rp**3
        rp = sqrt(xp**2+ys**2+zp**2)
        rp = max(rp,self%minrad)
        self%patch%force_per_unit_mass(ix,iy,iz,2) = &
        self%patch%force_per_unit_mass(ix,iy,iz,2) - newton*ys/rp**3
        rp = sqrt(xp**2+yp**2+zs**2)
        rp = max(rp,self%minrad)
        self%patch%force_per_unit_mass(ix,iy,iz,3) = &
        self%patch%force_per_unit_mass(ix,iy,iz,3) - newton*zs/rp**3
      end do
    end do
  end do
  end associate
  end associate
END SUBROUTINE var_grav

END MODULE gravity_mod

max_stars_repo_name: applejwjcat/dispatch

max_issues_repo_name: applejwjcat/dispatch

max_forks_repo_name: applejwjcat/dispatch


text:
program ex

  use gabriel
  use MPI
  use iso_fortran_env
  implicit none

  namelist/halo_exchange/array_size,nproc_row,halo_width,warmloop,timeloop,vars

  integer :: array_size=50
  integer :: nproc_row=8
  integer :: halo_width=2
  integer :: warmloop=100
  integer :: timeloop=1000
  integer :: vars=4
  real*8  :: gbpers

  real,dimension(:,:,:),allocatable :: a,b,c,e
  real :: dummyreal
  integer ierr,rank,right,left,mpisize,i,j,k
  integer hor,ver,pow2,s,realsize,realtype
  real*8 t0, t1,t,datasize

  type(box) :: bo
  type(distribution) :: d
  call MPI_Init(ierr)
  call MPI_Comm_rank(MPI_COMM_WORLD,rank,ierr)
  call MPI_Comm_size(MPI_COMM_WORLD,mpisize,ierr)
  right=rank+1
  left=rank-1
  if (right.ge.mpisize)right=0
  if (left.lt.0)left=mpisize-1
!  print*,'left,rank,right=',left,rank,right

  open(11, file='input')
  read(11, nml = halo_exchange)
  s=array_size
  if(rank.eq.0)then
    write(output_unit,nml=halo_exchange)
    print*,'MPI size: ',mpisize
  endif
  if (vars.lt.1.or.vars.gt.4)then
    print*,'Vars should be between 1 and 4'
    call MPI_Abort(MPI_COMM_WORLD,1,ierr)
  endif

  hor=mod(rank,nproc_row)
  ver=rank/nproc_row
  allocate(a(hor*s-halo_width:hor*s+s+halo_width-1,ver*s-halo_width+1:ver*s+s+halo_width,s))
  if(vars.ge.2)allocate(b(hor*s-halo_width:hor*s+s+halo_width-1,ver*s-halo_width+1:ver*s+s+halo_width,s))
  if(vars.ge.3)allocate(c(hor*s-halo_width:hor*s+s+halo_width-1,ver*s-halo_width+1:ver*s+s+halo_width,s))
  if(vars.ge.4)allocate(e(hor*s-halo_width:hor*s+s+halo_width-1,ver*s-halo_width+1:ver*s+s+halo_width,s))
  a=rank
  if(vars.ge.2)b=10+rank
  call zet_array(a,halo_width)
  if(vars.ge.2)call zet_array(b,halo_width)
  if(vars.ge.3)call zet_array(c,halo_width)
  if(vars.ge.4)call zet_array(e,halo_width)

!  write(*,'(a,4i3)')'Size,Rank,hor,ver=',mpisize,rank,hor,ver

  call gabriel_init

  call bo%init(a,(/hor*s,ver*s+1,1/),(/hor*s+s-1,ver*s+s,s/),MPI_COMM_WORLD,periodic=(/.true.,.true.,.true./))
  call d%halo(bo)
!  call d%autocreate(a,(/hor*s,ver*s+1,1/),(/hor*s+s-1,ver*s+s,6/),MPI_COMM_WORLD,periodic=(/.true.,.false.,.false./))
  if(vars.ge.2) then
    call d%joined(vars)
    call d%joined_add(a)
    call d%joined_add(b)
    if(vars.ge.3)call d%joined_add(c)
    if(vars.ge.4)call d%joined_add(e)
  endif
  call d%create()

  if(rank.eq.0)print*,'Warming up..'
  if (vars.ge.2) then
    do i=1,warmloop
      call d%update()
    enddo
  else
    do i=1,warmloop
      call d%update(a)
    enddo
  endif

  if(rank.eq.0)print*,'Apply distribution..'
  call gabriel_disable_checking

  t=0.0
  pow2=1
  if (vars.ge.2) then
  do i=1,timeloop
    if (rank.eq.0.and.i.eq.pow2) then
      write(*,'(i5,a1)')i,' '
      pow2=pow2*2
    endif
    call MPI_Barrier(MPI_COMM_WORLD,ierr)
    t0=MPI_Wtime()
    call d%update()

    t1=MPI_Wtime()
    t=t+(t1-t0)
!    a=0.0
    call zet_array(a,halo_width)
    call zet_array(b,halo_width)
    if(vars.ge.3)call zet_array(c,halo_width)
    if(vars.ge.4)call zet_array(e,halo_width)
  enddo
  else
  do i=1,timeloop
    if (rank.eq.0.and.i.eq.pow2) then
      write(*,'(i5,a1)')i,' '
      pow2=pow2*2
    endif
    call MPI_Barrier(MPI_COMM_WORLD,ierr)
    t0=MPI_Wtime()
    call d%update(a)

    t1=MPI_Wtime()
    t=t+(t1-t0)
!    a=0.0
    call zet_array(a,halo_width)
  enddo
  endif

  call MPI_Type_create_f90_real(precision(dummyreal),exponent(dummyreal),realtype,ierr)
  call MPI_Type_size(realtype,realsize,ierr)
!          vars                                      directions*type_size)/kb
  datasize=vars*(halo_width*s*s+halo_width*halo_width*s)*4.d0*realsize/1024.d0
  if(rank.eq.0)then
    call MPI_Reduce(MPI_IN_PLACE,t,1,MPI_DOUBLE_PRECISION,MPI_MAX,0,MPI_COMM_WORLD,ierr)
  else
    call MPI_Reduce(t,t,1,MPI_DOUBLE_PRECISION,MPI_MAX,0,MPI_COMM_WORLD,ierr)
  endif
  if(rank.eq.0)then
    print*,'Bytes per value         :',realsize
    print*,'Data size per comm (kB) :',datasize
    print*,'Comm. time (s)          :',t
    gbpers=datasize*timeloop/(1024.d0*1024.d0*t)
    write(*,'(a,2f10.2)')'Comm. speed (GB/s, GBit/s)      :',gbpers,gbpers*8
  endif

!  write(*,'(a,i3,a,i4,4f13.3)')' AFTER Rank',rank,' data=',ver*s,a(:,ver*s,1)
!  write(*,'(a,i3,a,i4,4f13.3)')' AFTER Rank',rank,' data=',ver*s+1,a(:,ver*s+1,1)
!  write(*,'(a,i3,a,i4,4f13.3)')' AFTER Rank',rank,' data=',ver*s+2,a(:,ver*s+2,1)
!  write(*,'(a,i3,a,i4,4f13.3)')' AFTER Rank',rank,' data=',ver*s+3,a(:,ver*s+3,1)
  deallocate(a)
  if(vars.ge.2)deallocate(b)
  if(vars.ge.3)deallocate(c)
  if(vars.ge.4)deallocate(e)
  call MPI_Finalize(ierr)

contains

  subroutine zet_array(a,halo_width)
    real,dimension(:,:,:),allocatable :: a
    integer halo_width
    integer ii,j,k

    a=0.0
    do k=lbound(a,3)+halo_width,ubound(a,3)-halo_width
    do j=lbound(a,2)+halo_width,ubound(a,2)-halo_width
    do ii=lbound(a,1)+halo_width,ubound(a,1)-halo_width
      a(ii,j,k)=a(ii,j,k)+ii*0.1+j*0.01+k*0.001
    enddo
    enddo
    enddo
  end subroutine zet_array

  logical function check_array_halo(a,halo_width)
    real,dimension(:,:,:),allocatable :: a,b
    integer :: nprocs,nproc_row,row,col,halo_width
    logical :: check
    integer ii,j,k

    check=.true.
    do k=lbound(a,3),ubound(a,3)
    do j=lbound(a,2),ubound(a,2)
    do ii=lbound(a,1),ubound(a,1)
      if (ii.lt.lbound(a,1)+halo_width .or.ii.gt.lbound(a,1)-halo_width.or. &
     &    j.lt.lbound(a,2)+halo_width .or.j.gt.lbound(a,2)-halo_width.or.&
     &    k.lt.lbound(a,3)+halo_width .or.k.gt.lbound(a,3)-halo_width)then
!         point is part of the halo
        check=check.and.(a(ii,j,k).eq.ii*0.1+j*0.01+k*0.001)
      endif
    enddo
    enddo
    enddo
    check_array_halo=check
  end function check_array_halo
end

max_stars_repo_name: jdonners/gabriel

max_issues_repo_name: jdonners/gabriel

max_forks_repo_name: jdonners/gabriel


text:
! { dg-do run }
! { dg-options "-fdefault-integer-8" }
!
! NOTE: This test is identical to negative_unit.f except -fdefault-integer-8
!
! PR libfortran/20660 and other bugs (not filed in bugzilla) relating
! to negative units
! PR 33055 Runtime error in INQUIRE unit existance with -fdefault-integer-8
! Test case update by Jerry DeLisle <jvdelisle@gcc.gnu.org>
!
! Bugs submitted by Walt Brainerd
      integer i
      integer, parameter ::ERROR_BAD_UNIT = 5005
      logical l

      i = -1
! gfortran created a 'fort.-1' file and wrote "Hello" in it
      write (unit=i, fmt=*, iostat=i) "Hello"
      if (i <= 0) STOP 1

      i = -11
      open (unit=i, file="xxx", iostat=i)
      if (i <= 0) STOP 2

      i = -42
      inquire (unit=i, exist=l)
      if (l) STOP 3

      i = 2_8*huge(0_4)+20_8
! This one is nasty
      inquire (unit=i, exist=l, iostat=i)
      if (l) STOP 4
      if (i.ne.0) STOP 5

      end


text:
Program testkind

! int8, int16, int32, int64
!    Parámetros de tipo Kind para especificar un tipo "INTEGER" con un tamaño de
! almacenamiento de 8, 16, 32 y 64 bits. Es negativo si una plataforma de
! destino no es compatible con el tipo particular. (Fortran 2008 o posterior).

! real32, real64, real128
!    Parámetros de tipo Kind para especificar un tipo "REAL" con un tamaño de
! almacenamiento de 32, 64 y 128 bits. Es negativo si una plataforma de
! destino no es compatible con el tipo particular. (Fortran 2008 o posterior).

use ISO_FORTRAN_ENV

implicit none

integer(kind = int8), parameter :: is = int8, id = int16, il = int32, ix = int64
integer(kind = int8), parameter :: rs = real32, rd = real64, rl = real128
! s : Short, d : Double, l : Large, x : eXtralarge

integer(kind = is) :: el_entero
integer(is) :: otro_entero
real(kind = rl) :: los_reales
real(rl) :: otros_reales
complex(kind = rl) :: el_complejado
complex(rl) :: otro_complejado

el_entero = 8
los_reales = 10.00_rl
otros_reales = 355.00_rl / 113.00_rl ! Almost pi
el_complejado = (355.00_rl, 113.00_rl)
otro_complejado = (3.00_rl, 4.00_rl)

write(*,*) el_entero, los_reales
write(*,*) otro_entero, otros_reales
write(*,*) el_complejado, otro_complejado
write(*,*) real(el_complejado)/imag(el_complejado), abs(otro_complejado)

write(*,*) is, id, il, ix
write(*,*) rs, rd, rl

end program testkind

max_stars_repo_name: EdgardoBonzi/Fortran-Examples

max_issues_repo_name: EdgardoBonzi/Fortran-Examples

max_forks_repo_name: EdgardoBonzi/Fortran-Examples


text:
program test_union
use iso_fortran_env, only:int8, int16, int32, int64, real32, real64, real128
use mod_testing, only:assert, initialize_tests, report_tests
use mod_functional

implicit none

logical, dimension(:), allocatable :: tests
logical :: test_failed
integer :: n, norder, ntests
integer, parameter :: stdout = 6

n = 1
ntests = 15
call initialize_tests(tests, ntests)

tests(n) = assert(all(union([1_int8, 2_int8], [2_int8, 3_int8]) == [1, 2, 3]), &
                  'union,  int8')
n = n + 1

tests(n) = assert(all(union([1_int16, 2_int16], [2_int16, 3_int16]) == [1, 2, 3]), &
                  'union,  int16')
n = n + 1

tests(n) = assert(all(union([1_int32, 2_int32], [2_int32, 3_int32]) == [1, 2, 3]), &
                  'union,  int32')
n = n + 1

tests(n) = assert(all(union([1_int64, 2_int64], [2_int64, 3_int64]) == [1, 2, 3]), &
                  'union,  int64')
n = n + 1

tests(n) = assert(all(union([1._real32, 2._real32], [2._real32, 3._real32]) == [1, 2, 3]), &
                  'union,  real32')
n = n + 1

tests(n) = assert(all(union([1._real64, 2._real64], [2._real64, 3._real64]) == [1, 2, 3]), &
                  'union,  real64')
n = n + 1

tests(n) = assert(all(union([1._real128, 2._real128], [2._real128, 3._real128]) == [1, 2, 3]), &
                  'union,  real128')
n = n + 1

tests(n) = assert(all(union([cmplx(1._real32, 0._real32), cmplx(2._real32, 0._real32)], &
                            [cmplx(2._real32, 0._real32), cmplx(3._real32, 0._real32)])&
                      == [cmplx(1._real32, 0._real32), &
                          cmplx(2._real32, 0._real32), &
                          cmplx(3._real32, 0._real32)]), &
                  'union,  complex real32')
n = n + 1

tests(n) = assert(all(union([cmplx(1._real64, 0._real64), cmplx(2._real64, 0._real64)], &
                            [cmplx(2._real64, 0._real64), cmplx(3._real64, 0._real64)])&
                      == [cmplx(1._real64, 0._real64), &
                          cmplx(2._real64, 0._real64), &
                          cmplx(3._real64, 0._real64)]), &
                  'union,  complex real64')
n = n + 1

tests(n) = assert(all(union([cmplx(1._real128, 0._real128), cmplx(2._real128, 0._real128)], &
                            [cmplx(2._real128, 0._real128), cmplx(3._real128, 0._real128)])&
                      == [cmplx(1._real128, 0._real128), &
                          cmplx(2._real128, 0._real128), &
                          cmplx(3._real128, 0._real128)]), &
                  'union,  complex real128')
n = n + 1

tests(n) = assert(union('Hello', 'world') == 'Helowrd', &
                  'union, character string')
n = n + 1

tests(n) = assert(all(union(arange(1, 0), arange(1, 0)) == arange(1, 0)), &
                  'union of empty arrays is an empty array')
n = n + 1

tests(n) = assert(all(union([1, 2, 2, 3], arange(1, 0)) == set([1, 2, 2, 3])), &
                  'union(x, []) == set(x)')
n = n + 1

tests(n) = assert(all(union(arange(1, 0), [1, 2, 2, 3]) == set([1, 2, 2, 3])), &
                  'union([], x) == set(x)')
n = n + 1

tests(n) = assert(all(union([1, 2], [3, 4]) == ([1, 2].union.[3, 4])), &
                  'union operator,  a.union.b')
n = n + 1

test_failed = .false.
call report_tests(tests, test_failed)
if(test_failed)stop 1

end program test_union

max_stars_repo_name: econwang/functional-fortran

max_issues_repo_name: econwang/functional-fortran

max_forks_repo_name: econwang/functional-fortran


text:
c-----------------------------------------------------------------------------
c
c Copyright (c) 2000 - 2014, Lawrence Livermore National Security, LLC
c Produced at the Lawrence Livermore National Laboratory
c LLNL-CODE-442911
c All rights reserved.
c
c This file is  part of VisIt. For  details, see https://visit.llnl.gov/.  The
c full copyright notice is contained in the file COPYRIGHT located at the root
c of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
c
c Redistribution  and  use  in  source  and  binary  forms,  with  or  without
c modification, are permitted provided that the following conditions are met:
c
c  - Redistributions of  source code must  retain the above  copyright notice,
c    this list of conditions and the disclaimer below.
c  - Redistributions in binary form must reproduce the above copyright notice,
c    this  list of  conditions  and  the  disclaimer (as noted below)  in  the
c    documentation and/or other materials provided with the distribution.
c  - Neither the name of  the LLNS/LLNL nor the names of  its contributors may
c    be used to endorse or promote products derived from this software without
c    specific prior written permission.
c
c THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT  HOLDERS AND CONTRIBUTORS "AS IS"
c AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO, THE
c IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR  PURPOSE
c ARE  DISCLAIMED. IN  NO EVENT  SHALL LAWRENCE  LIVERMORE NATIONAL  SECURITY,
c LLC, THE  U.S.  DEPARTMENT OF  ENERGY  OR  CONTRIBUTORS BE  LIABLE  FOR  ANY
c DIRECT,  INDIRECT,   INCIDENTAL,   SPECIAL,   EXEMPLARY,  OR   CONSEQUENTIAL
c DAMAGES (INCLUDING, BUT NOT  LIMITED TO, PROCUREMENT OF  SUBSTITUTE GOODS OR
c SERVICES; LOSS OF  USE, DATA, OR PROFITS; OR  BUSINESS INTERRUPTION) HOWEVER
c CAUSED  AND  ON  ANY  THEORY  OF  LIABILITY,  WHETHER  IN  CONTRACT,  STRICT
c LIABILITY, OR TORT  (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY  WAY
c OUT OF THE  USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
c DAMAGE.
c
c-----------------------------------------------------------------------------

      program main
      implicit none
      include "silo.inc"
      integer dbfile, ierr

c Create the Silo file
      ierr = dbcreate("fucdvar2d.silo", 14, DB_CLOBBER, DB_LOCAL,
     .                "Unstructured 2D mesh with variable", 34,
     .                 DB_HDF5, dbfile)
      if(dbfile.eq.-1) then
          write (6,*) 'Could not create Silo file!\n'
          goto 10000
      endif
c Add other Silo calls here.
      call write_ucd2d(dbfile)
      call write_ucdvars(dbfile)
c Close the Silo file.
      ierr = dbclose(dbfile)
10000 stop
      end

c
c   *---*---*--*   nodes (5,6,7,8)
c   |   |\5/ 3 |
c   | 1 |4*----*   nodes (3,4)
c   |   |/  2  |
c   *---*------*   nodes(0,1,2)
c
c   cell 1 quad(0,1,6,5)
c   cell 2 quad(1,2,4,3)
c   cell 3 quad(3,4,8,7)
c   cell 4 tri(1,3,6)
c   cell 5 tri(3,7,6)
c
      subroutine write_ucd2d(dbfile)
      implicit none
      integer dbfile
      include "silo.inc"
      integer err, ierr, ndims, nshapetypes, nnodes, nzones
c Node coordinates
      real x(9) /0., 2., 5., 3., 5., 0., 2., 4., 5./
      real y(9) /0., 0., 0., 3., 3., 5., 5., 5., 5./
c Connectivity
      integer LNODELIST
      parameter (LNODELIST = 18)
      integer nodelist(LNODELIST) /2,4,7,
     .  4,8,7,
     .  1,2,7,6,
     .  2,3,5,4,
     .  4,5,9,8/
c Shape type 1 has 3 nodes (tri), shape type 2 is quad
      integer shapesize(2) /3, 4/
c We have 2 tris and 3 quads
      integer shapecounts(2) /2, 3/
      nshapetypes = 2
      nnodes = 9
      nzones = 5
      ndims = 2
c Write out connectivity information.
      err = dbputzl(dbfile, "zonelist", 8, nzones, ndims, nodelist,
     . LNODELIST, 1, shapesize, shapecounts, nshapetypes, ierr)
c Write an unstructured mesh
      err = dbputum(dbfile, "mesh", 4, ndims, x, y, DB_F77NULL,
     . "X", 1, "Y", 1, DB_F77NULL, 0, DB_FLOAT, nnodes, nzones,
     . "zonelist", 8, DB_F77NULL, 0, DB_F77NULL, ierr)
      end

      subroutine write_ucdvars(dbfile)
      implicit none
      integer dbfile
      include "silo.inc"
      integer err, ierr, NNODES, NZONES, optlistid
      parameter (NNODES = 9)
      parameter (NZONES = 5)
      real nodal(NNODES) /1.,2.,3.,4.,5.,6.,7.,8.,9./
      real zonal(NZONES) /1.,2.,3.,4.,5./
      err = dbmkoptlist(1, optlistid)
      err = dbaddcopt(optlistid, DBOPT_UNITS, "g/cc", 4)
c Write a zone-centered variable.
      err = dbputuv1(dbfile, "zonal", 5, "mesh", 4, zonal, NZONES,
     . DB_F77NULL, 0, DB_FLOAT, DB_ZONECENT, optlistid, ierr)
      err = dbfreeoptlist(optlistid)
c Write a node-centered variable.
      err = dbputuv1(dbfile, "nodal", 5, "mesh", 4, nodal, NNODES,
     . DB_F77NULL, 0, DB_FLOAT, DB_NODECENT, DB_F77NULL, ierr)
      end