!utilities and look up tables for lj fitting trialanine
MODULE utils
use graphf
implicit none

character(len=4), allocatable, dimension(:) :: ordering_array, lj_species
character(len=4), allocatable, dimension(:) :: opt_species, o_f_species, all_o_f
real*8, allocatable, dimension(:,:) :: chff_lj_params
real*8, allocatable, dimension(:,:) :: o_f_array
real*8, allocatable, dimension(:,:,:) :: crd_data, dist_array
real*8, allocatable, dimension(:) :: ref_energies
character(len=200), allocatable, dimension(:) :: crd_names
integer :: natoms, nfiles, nspecies, nbonds, nonefour, nopt, npep_atoms
integer :: num_one_four
integer, allocatable, dimension(:,:) :: bond_array
integer, allocatable, dimension(:,:) :: excl_array
integer, allocatable, dimension(:) :: stan_lj_index, spec_lj_index
logical, allocatable, dimension(:) :: is_opt_arr
real*8 :: r_off, r_on, upper_bound, lower_bound
character(len=3), allocatable, dimension(:) :: print_helper

public ordering_array, lj_species, bond_array, nbonds, o_f_species, dist_array
public natoms, nfiles, nspecies, chff_lj_params, nopt, crd_data, npep_atoms
public excl_array, o_f_array, r_on, r_off, opt_species, ref_energies, crd_names
public stan_lj_index, spec_lj_index, is_opt_arr, upper_bound, lower_bound
public nonefour, num_one_four, all_o_f, print_helper

contains
! initializes parameters of system

subroutine init_params(num_files, cut_on, cut_off, num_pep_atoms, n_onefour,&
                        l_bound, u_bound)
    implicit none
    integer :: num_atoms, num_files, num_pep_atoms, num_species, n_onefour
    real*8 :: cut_on, cut_off, l_bound, u_bound
    nfiles = num_files
    r_on = cut_on
    r_off = cut_off
    npep_atoms = num_pep_atoms
    num_one_four = n_onefour
    upper_bound = u_bound
    lower_bound = l_bound

end subroutine init_params

subroutine read_crd_file(file_name, crd_conf)
    implicit none
    character(len=100) :: file_name, junk
    real*8, allocatable, dimension(:,:) :: crd_conf
    integer :: i, num_at
    if(.NOT. allocated(crd_conf)) allocate(crd_conf(natoms, 3))
    open(69, file=file_name, status = 'old')
    read(69,*)
    read(69,*)
    read(69,*)
    read(69,*) num_at, junk
    do i = 1, num_at
        read(69,*) junk, junk, junk, junk, crd_conf(i,1), crd_conf(i,2), crd_conf(i,3)
    end do

    close(69)
end subroutine read_crd_file

subroutine load_data(crd_name_path, ref_file)
    implicit none
    character(len=100) :: crd_name_path, enum, file_name, ref_file
    real*8, allocatable, dimension(:,:) :: crd_conf
    integer :: i, iatom, idir
    allocate(crd_data(natoms, 3, nfiles))
    allocate(crd_names(nfiles))
    open(69, file=crd_name_path, status = 'old')
    do i = 1, nfiles
        read(69,'(A)') crd_names(i)
    end do
    crd_data = 0
    do i = 1, nfiles
        file_name = trim(crd_names(i))
        call read_crd_file(file_name, crd_conf)
        do iatom = 1, natoms
            do idir = 1, 3
                crd_data(iatom, idir, i) = crd_conf(iatom, idir)
            end do
        end do
    end do
    close(69)
    open(69, file=ref_file, status = 'old')
    allocate(ref_energies(nfiles))
    do i = 1, nfiles
        read(69, *) ref_energies(i)
    end do
    close(69)
    call crt_dist_array()
end subroutine load_data

! loads relevant parameters like atom ordering and lj species
subroutine load_lj_params(file_name)
    implicit none
    character(len=100) :: file_name
    integer :: i
    real*8 :: junk
    open(69, file=file_name, status = 'old')
    read(69, *) nspecies
    if (.not. allocated(lj_species)) allocate(lj_species(nspecies))
    if (.not. allocated(chff_lj_params)) allocate(chff_lj_params(nspecies, 2))
    do i = 1, nspecies
       read(69,*) lj_species(i), junk, chff_lj_params(i,1), chff_lj_params(i,2)
    end do
    close(69)
end subroutine load_lj_params

subroutine load_one_four_params(file_name, name_file)
    implicit none
    character(len=100) :: file_name, junk, name_file
    integer :: i
    open(69, file=file_name, status = 'old')
    read(69, *) junk
    if (.not. allocated(o_f_array)) allocate(o_f_array(nspecies, 2))
    if (.not. allocated(all_o_f)) allocate(all_o_f(nspecies))
    do i = 1, nspecies
        read(69, *) all_o_f(i), o_f_array(i, 1), o_f_array(i, 2)
    end do
    close(69)
    open(69, file=name_file, status = 'old')
    read(69,*) nonefour
    if (.not. allocated(o_f_species)) allocate(o_f_species(nonefour))
    do i = 1, nonefour
        read(69,*) o_f_species(i)
    end do
    close(69)
end subroutine load_one_four_params

subroutine load_psf(file_name)
    implicit none
    character(len=100) :: file_name, junk
    integer :: i
    open(69, file=file_name, status = 'old')
    do i = 1, 6
        read(69,*)
    end do
    read(69, *) natoms
    if (.not. allocated(ordering_array)) allocate(ordering_array(natoms))
    do i = 1, natoms
        read(69,*) junk, junk, junk, junk, junk, ordering_array(i)
    end do
    close(69)
end subroutine load_psf

subroutine load_opt_spec(file_name)
    implicit none
    character(len=100) file_name
    integer :: i
    open(69, file=file_name, status = 'old')
    read(69, *) nopt
    allocate(opt_species(nopt))
    do i = 1, nopt
        read(69, *) opt_species(i)
    end do
end subroutine load_opt_spec

! converts atom index in psf file to its LJ label
character(len=3) function at_to_label(iatom)
    integer :: iatom
    at_to_label = ordering_array(iatom)
end function at_to_label

! converts atom LJ label to index in LJ array
integer function lab_to_lj_i(label)
    character(len=3) :: label
    integer :: i
    do i = 1, nspecies
        if (trim(label) == lj_species(i)) then
            lab_to_lj_i = i
        end if
    end do
end function lab_to_lj_i

! calculates euclidian distance between crd1 and crd2
real*8 function get_distance(crd1, crd2)
    real*8, dimension(3) :: crd1, crd2
    get_distance = NORM2(crd1 - crd2)
end function  get_distance

! calculates lj interaction energy between a pair of atoms at distance dist
real*8 function calc_lj_pair(eps1, rmin1, eps2, rmin2, dist)
    real*8 :: eps1, rmin1, eps2, rmin2, epsij, rminij, dist

    epsij = SQRT(eps1 * eps2)
    rminij = rmin1 + rmin2
    calc_lj_pair = epsij*((rminij/dist)**12 - 2*((rminij/dist)**6))
end function calc_lj_pair

logical function is_opt(iatom)
    integer :: iatom, i
    do i=1, nopt
        if (at_to_label(iatom) == opt_species(i)) then
            is_opt = .true.
            return
        else
            is_opt = .false.
        end if
    end do
end function is_opt

logical function is_o_f(iatom)
    integer :: iatom, i
    do i = 1, nonefour
        if (at_to_label(iatom) == o_f_species(i)) then
            is_o_f = .true.
            return
        else
            is_o_f = .false.
        end if
    end do
end function is_o_f

integer function label_to_opt_index(label)
    character(len=3) :: label
    integer :: i
    do i = 1, nopt
        if(trim(label) == opt_species(i)) then
            label_to_opt_index = i
            return
        end if
    end do
end function label_to_opt_index

integer function label_to_x_vec_o_f(label)
    character(len=3) :: label
    integer :: i
    do i = nopt+1, nopt+num_one_four
        if (label == print_helper(i)) then
            label_to_x_vec_o_f = i
        end if
    end do

end function label_to_x_vec_o_f

subroutine calc_look_ups()
    implicit none
    integer :: iatom, lj_index, nopt_atoms
    integer :: i, j
    if (.not. allocated(is_opt_arr)) allocate(is_opt_arr(natoms))
    do iatom = 1, natoms
       is_opt_arr(iatom) = is_opt(iatom)
    end do
    if (.not. allocated(stan_lj_index)) allocate(stan_lj_index(natoms))
    do iatom = 1, natoms
        if (is_opt_arr(iatom) .eqv. .true.) then
            stan_lj_index(iatom) = label_to_opt_index(at_to_label(iatom))
        else
            stan_lj_index(iatom) = lab_to_lj_i(at_to_label(iatom))
        end if
    end do
end subroutine calc_look_ups

real*8 function get_eps_stand(iatom, x)
    integer :: iatom, lj_index
    logical :: optimize
    real*8, dimension(2*(nopt+num_one_four)) :: x
    optimize = is_opt_arr(iatom)
    if (optimize .eqv. .true.) then
        lj_index = stan_lj_index(iatom)
        get_eps_stand = -1.0d0*abs(x(2*lj_index-1))
        return
    else
        lj_index = stan_lj_index(iatom)
        get_eps_stand = chff_lj_params(lj_index, 1)
        return
    end if
end function get_eps_stand

real*8 function get_rmin_stand(iatom, x)
    integer :: iatom, lj_index
    logical :: optimize
    real*8, dimension(2*(nopt+num_one_four)) :: x
    optimize = is_opt_arr(iatom)
    if (optimize .eqv. .true.) then
        lj_index = stan_lj_index(iatom)
        get_rmin_stand = abs(x(2*lj_index))
        return
    else
        lj_index = stan_lj_index(iatom)
        get_rmin_stand = chff_lj_params(lj_index, 2)
        return
    end if
end function get_rmin_stand

real*8 function get_eps_spec(iatom, x)
    integer :: iatom, lj_index
    logical :: optimize
    real*8, dimension(2*(nopt+num_one_four)) :: x
    optimize = is_opt_arr(iatom)
    if ((optimize .eqv. .true.) .and. (is_o_f(iatom) .eqv. .true.)) then
        lj_index = label_to_x_vec_o_f(at_to_label(iatom))
        get_eps_spec = -1*abs(x(2*lj_index-1))
        return
    else if ((optimize .eqv. .true.) .and. (is_o_f(iatom) .eqv. .false.)) then
        lj_index = stan_lj_index(iatom)
        get_eps_spec = -1*abs(x(2*lj_index-1))
        return
    else
       lj_index = stan_lj_index(iatom)
        get_eps_spec = o_f_array(lj_index, 1)
    end if
end function get_eps_spec

real*8 function get_rmin_spec(iatom, x)
    integer :: iatom, lj_index
    logical :: optimize
    real*8, dimension(2*(nopt+num_one_four)) :: x
    optimize = is_opt_arr(iatom)
    if ((optimize .eqv. .true.) .and. (is_o_f(iatom) .eqv. .true.)) then
        lj_index = label_to_x_vec_o_f(at_to_label(iatom))
        get_rmin_spec = abs(x(2*lj_index))
        return
    else if ((optimize .eqv. .true.) .and. (is_o_f(iatom) .eqv. .false.)) then
        lj_index = stan_lj_index(iatom)
        get_rmin_spec = abs(x(2*lj_index))
        return
    else
        lj_index = stan_lj_index(iatom)
        get_rmin_spec = o_f_array(lj_index, 2)
    end if
end function get_rmin_spec

subroutine get_lj_energy(ifile, energy, x)
    implicit none
    real*8 :: energy, dist_ij, curr_lj_energy, eps1, eps2, rmin1, rmin2
    integer :: iatom, jatom, bond_dist, i, ifile
    real*8, dimension(2*(nopt+num_one_four)) :: x, curr_sol
    real*8 :: t1, t2
    energy = 0.0
    curr_lj_energy = 0.0
    do iatom = 1, npep_atoms
        do jatom = iatom + 1, natoms
            bond_dist = excl_array(iatom, jatom)
            dist_ij = dist_array(iatom, jatom, ifile)
            if (bond_dist .ge. 4) then
                eps1 = get_eps_stand(iatom, x)
                eps2 = get_eps_stand(jatom, x)
                rmin1 = get_rmin_stand(iatom, x)
                rmin2 = get_rmin_stand(jatom, x)
                if ((eps1 .gt. 0.0d0) .or. (eps2 .gt. 0.0d0)) then
                do i = 1, nopt+5
                    write(*,*) "ERROR: EPS larger than 0"
                    write(*,*) x(2*i-1), x(2*i)
                end do
                end if
                curr_lj_energy = calc_lj_pair(eps1, rmin1, eps2, rmin2, dist_ij)
                curr_lj_energy = curr_lj_energy*vswitch(dist_ij, r_on, r_off)
            else if (bond_dist == 3) then
                eps1 = get_eps_spec(iatom, x)
                eps2 = get_eps_spec(jatom, x)
                rmin1 = get_rmin_spec(iatom, x)
                rmin2 = get_rmin_spec(jatom, x)
                curr_lj_energy = calc_lj_pair(eps1, rmin1, eps2, rmin2, dist_ij)
                curr_lj_energy = curr_lj_energy*vswitch(dist_ij, r_on, r_off)
            end if
            energy = energy + curr_lj_energy
            curr_lj_energy = 0.0
        end do
    end do
end subroutine get_lj_energy

real*8 function vswitch(dist, r_on, r_off)
    real*8 :: dist, r_on, r_off

    if (dist .le. r_on) then
        vswitch = 1.0
        return
    else if ((dist .gt. r_on) .and. (dist .le. r_off)) then
        vswitch = (r_off**2 - dist**2)**2*(r_off**2 + 2*dist**2 - 3*r_on**2)
        vswitch = vswitch / (r_off**2 - r_on **2)**3
        return
    else if (dist .gt. r_off) then
        vswitch = 0.0
        return
    end if
end function vswitch

subroutine crt_dist_array()
    implicit none
    integer :: iatom, jatom, ifile
    if (.not. allocated(dist_array)) allocate(dist_array(natoms, natoms, nfiles))
      do ifile = 1, nfiles
        do iatom = 1, natoms
          do jatom = 1, natoms
            dist_array(iatom, jatom, ifile) = get_distance(crd_data(iatom, :, ifile),&
                                                           crd_data(jatom, :, ifile))
          end do
        end do
      end do
end subroutine crt_dist_array

subroutine get_excl_array(bond_file)
    implicit none
    type(graph) :: G
    real*8 :: distance
    character(len=100) :: bond_file
    integer :: i, j, start, sto
    logical :: bonded

    call G%setOrder(natoms)
    open(69, file=bond_file, status = 'old')
    read(69,*) nbonds
    if(.not. allocated(bond_array)) allocate(bond_array(nbonds, 2))
    do i = 1, nbonds
        read(69, *) bond_array(i, 1), bond_array(i, 2)
    end do
    do i = 1, nbonds
        call G%addEdge(bond_array(i, 1), bond_array(i, 2))
    end do
    if (.not. allocated(excl_array)) allocate(excl_array(natoms, natoms))
    do i = 1, natoms
        do j = 1, natoms
            call exampleDijkstra(G, i, j, distance, bonded)
            if (bonded .eqv. .TRUE.) then
                excl_array(i,j) = INT(distance)
            else
                excl_array(i,j) = 1000
            end if
        end do
    end do
end subroutine get_excl_array

subroutine init_print_helper()
    implicit none

    integer :: i, j, k

    if (.not. allocated(print_helper)) allocate(print_helper(nopt+num_one_four))
    print_helper(1:nopt) = opt_species
    k = 1
    do i = nopt+1, nopt+num_one_four
        do j = k, nopt
            if (check_o_f(print_helper(j)) .eqv. .true.) then
                print_helper(i) = print_helper(j)
                k = j+1
                exit
            end if
        end do
    end do
end subroutine init_print_helper

logical function check_o_f(label)
    integer :: i
    character(len=3) :: label
    do i = 1, nonefour
        if (label == o_f_species(i)) then
            check_o_f = .true.
            return
        else
            check_o_f = .false.
        end if
    end do
end function check_o_f
end module utils


text:
!
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!

program testieee18flushz
use ieee_arithmetic
real*4 a, ral, rone
real*4 aa1(140), aa2(140)
logical ll1(140), le1(140)
logical luflow
logical lfsav(5), lfset(5)

lfset = .false.
call ieee_get_halting_mode(ieee_all, lfsav)
call ieee_set_halting_mode(ieee_all, lfset)

call ieee_get_underflow_mode(luflow)
call ieee_set_underflow_mode(.false.)

rone = 1.0
rlogb = 0.0
a = rone
do j = 1, 140
  ral = ieee_logb(a)
  aa1(j) = ral
  if (j .le. 127) then
    aa2(j) = rlogb
  else
    aa2(j) = ieee_value(rlogb, ieee_negative_inf)
  end if
  i = int(ieee_rint(ral))
  a = a * 0.5
  rlogb = rlogb - 1.0
end do

do j = 1, 140
  ll1(j) = aa1(j) .eq. aa2(j)
end do

le1 = .true.

call check(ll1, le1, 140)

call ieee_set_underflow_mode(luflow)
call ieee_set_halting_mode(ieee_all, lfsav)
end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
!* Copyright (c) 1999, NVIDIA CORPORATION.  All rights reserved.
!*
!* 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
!*
!*     http://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.

!* Attribute oriented initializations using intrinsic functions

program e10

   parameter(NTEST=76)
   integer :: result(NTEST)
   integer :: expect(NTEST) = (/  &
   !idx1
       5,		&
   !idx2
       5,		&
   !idx3
       6,		&
   !idx4
       6,		&
   !idx5
       5,		&
   !idx6
       5,		&
   !idx7
       5,		&
   !idx8
       5,		&
   !idx9
       2,		&
   !idxarry1
       2, 3, 4, 5,		&
   !idxarry2
       2, 3,		&
   !lentrim1
       6,		&
   !lentrim2
       13,		&
   !lentrim3
       6,		&
   !lentrim4
       3,		&
   !lentrimarry1
       4, 4, 5, 6,		&
   !lentrimarry2
       3, 4,		&
   !scan1
       5,		&
   !scan2
       6,		&
   !scan3
       6,		&
   !scan4
       6,		&
   !scan5
       5,		&
   !scan6
       6,		&
   !scan7
       6,		&
   !scan8
       5,		&
   !scan9
       2,		&
   !scanarry1
       2, 3, 4, 5,		&
   !scanarry2
       2, 4,		&
   !vrfy1
       0,		&
   !vrfy2
       10,		&
   !vrfy3
       1,		&
   !vrfy4
       1,		&
   !vrfy5
       1,		&
   !vrfy6
       3,		&
   !vrfy7
       4,		&
   !vrfy8
       1,		&
   !vrfy9
       1,		&
   !vrfyarry1
       4, 0, 0, 0,		&
   !vrfyarry2
       4, 1,		&
   !char1
       65,		&
   !char2
       66,		&
   !chararry1
       48, 49, 50, 51, 52, 53,  &
    !index kind
       2,4,8,           &
    !scan kind
       2,4,8,           &
    !verify kind
       2,4,8,           &
    !len_trim kind
       2,4,8,           &
    !char kind
       1 &
   /)

   integer, parameter :: iparam = 8
   integer, parameter :: iparamarry(5) = (/(i,i=1,5)/)
   character(len=6), parameter :: strparamarry1(4) = &
                    (/"abcz  ", " abc  ", "  abc ", "   abc" /)
   character(len=6), parameter :: strparam1 = "   abc"
   character(len=2), parameter :: strparam2 = "bc"

   character(len=6), parameter :: strparam3 = "   abc"

  type t1
    character(len=6) :: str1
  end type
  type (t1), parameter:: t1_param1 = t1( adjustl("   abc") )

  integer :: idx1 = index("xxyybcwwvv", "bc", KIND=4)
  integer :: idx2 = index("xxyybcwwvv", "bc", .TRUE., KIND=4)
  integer :: idx3 = index("zxxyybcwwvvzz", "bc", .FALSE., KIND=4)
  integer :: idx4 = index("zxxyybcwwvvzz", strparam2, .FALSE., KIND=4)
  integer :: idx5 = index(strparam1, "bc", KIND=4)
  integer :: idx6 = index(strparam1, "bc", .TRUE., KIND=4)
  integer :: idx7 = index(strparam1, strparam2, .TRUE., KIND=4)
  integer :: idx8 = index(strparam1, "bc", .FALSE., KIND=4)
  integer :: idx9 = index(t1_param1%str1, "bc", .FALSE., KIND=8)
  integer :: idxarry1(4) = index(strparamarry1, "bc", KIND=4)
  integer :: idxarry2(2) = index((/"abc   ", " abc  "/), "bc", KIND=4)

  integer :: scan1 = scan("xxyybcwwvv", "bc", KIND=4)
  integer :: scan2 = scan("xxyybcwwvv", "bc", .TRUE., KIND=4)
  integer :: scan3 = scan("zxxyybcwwvvzz", "bc", .FALSE., KIND=4)
  integer :: scan4 = scan("zxxyybcwwvvzz", strparam2, .FALSE., KIND=4)
  integer :: scan5 = scan(strparam1, "bc", KIND=4)
  integer :: scan6 = scan(strparam1, "bc", .TRUE.,KIND=4)
  integer :: scan7 = scan(strparam1, strparam2, .TRUE.,KIND=4)
  integer :: scan8 = scan(strparam1, "bc", .FALSE.,KIND=4)
  integer :: scan9 = scan(t1_param1%str1, "bc", .FALSE.,KIND=4)
  integer :: scanarry1(4) = scan(strparamarry1, "bc",KIND=4)
  integer :: scanarry2(2) = scan((/"abc   ", " aac  "/), "bc",KIND=4)

  integer :: vrfy1 = verify("bbbbc", "bc", KIND=4)
  integer :: vrfy2 = verify("xxyybcwwvv", "xywzbc", .TRUE.,KIND=4)
  integer :: vrfy3 = verify("zxxyybcwwvvzz", "bc", .FALSE., KIND=4)
  integer :: vrfy4 = verify("zxxyybcwwvvzz", strparam2, .FALSE., KIND=4)
  integer :: vrfy5 = verify(strparam1, "abc", KIND=4)
  integer :: vrfy6 = verify(strparam1, "bac", .TRUE.,KIND=4)
  integer :: vrfy7 = verify(strparam1, strparam2, .TRUE.,KIND=4)
  integer :: vrfy8 = verify(strparam1, "abc", .FALSE.,KIND=4)
  integer :: vrfy9 = verify(t1_param1%str1, "bc", .FALSE.,KIND=4)
  integer :: vrfyarry1(4) = verify(strparamarry1, " abc",KIND=4)
  integer :: vrfyarry2(2) = verify((/"abc   ", " abc  "/), "abc",KIND=4)

  integer :: lentrim1 = len_trim("xxyybc    ",KIND=4)
  integer :: lentrim2 = len_trim("zxxyybcwwvvzz",KIND=4)
  integer :: lentrim3 = len_trim(strparam3,KIND=4)
  integer :: lentrim4 = len_trim(t1_param1%str1,KIND=4)
  integer :: lentrimarry1(4) = len_trim(strparamarry1,KIND=4)
  integer :: lentrimarry2(2) = len_trim((/"abc   ", " abc  "/),KIND=4)

  integer, parameter :: iparam2 = 66
  character :: char1 = char(65, KIND=1)
  character :: char2 = char(iparam2, KIND=1)
  character :: chararry1(6) = char((/48,49,50,51,52,53/),KIND=1)

!  print *,"! idx1"
!  print *,idx1;
  result(1) = idx1
!  print *,"! idx2"
!  print *,idx2;
  result(2) = idx2
!  print *,"! idx3"
!  print *,idx3;
  result(3) = idx3
!  print *,"! idx4"
!  print *,idx4;
  result(4) = idx4
!  print *,"! idx5"
!  print *,idx5;
  result(5) = idx5
!  print *,"! idx6"
!  print *,idx6;
  result(6) = idx6
!  print *,"! idx7"
!  print *,idx7;
  result(7) = idx7
!  print *,"! idx8"
!  print *,idx8;
  result(8) = idx8
!  print *,"! idx9"
!  print *,idx9;
  result(9) = idx9
!  print *,"! idxarry1"
!  print *,idxarry1;
  result(10:13) = idxarry1
!  print *,"! idxarry2"
!  print *,idxarry2;
  result(14:15) = idxarry2
!  print *,"! lentrim1"
!  print *,lentrim1;
  result(16) = lentrim1
!  print *,"! lentrim2"
!  print *,lentrim2;
  result(17) = lentrim2
!  print *,"! lentrim3"
!  print *,lentrim3;
  result(18) = lentrim3
!  print *,"! lentrim4"
!  print *,lentrim4;
  result(19) = lentrim4
!  print *,"! lentrimarry1"
!  print *,lentrimarry1;
  result(20:23) = lentrimarry1
!  print *,"! lentrimarry2"
!  print *,lentrimarry2;
  result(24:25) = lentrimarry2

!  print *,"! scan1"
!  print *,scan1;
  result(26) = scan1
!  print *,"! scan2"
!  print *,scan2;
  result(27) = scan2
!  print *,"! scan3"
!  print *,scan3;
  result(28) = scan3
!  print *,"! scan4"
!  print *,scan4;
  result(29) = scan4
!  print *,"! scan5"
!  print *,scan5;
  result(30) = scan5
!  print *,"! scan6"
!  print *,scan6;
  result(31) = scan6
!  print *,"! scan7"
!  print *,scan7;
  result(32) = scan7
!  print *,"! scan8"
!  print *,scan8;
  result(33) = scan8
!  print *,"! scan9"
!  print *,scan9;
  result(34) = scan9
!  print *,"! scanarry1"
!  print *,scanarry1;
  result(35:38) = scanarry1
!  print *,"! scanarry2"
!  print *,scanarry2;
  result(39:40) = scanarry2

!  print *,"! vrfy1"
!  print *,vrfy1;
  result(41) = vrfy1
!  print *,"! vrfy2"
!  print *,vrfy2;
  result(42) = vrfy2
!  print *,"! vrfy3"
!  print *,vrfy3;
  result(43) = vrfy3
!  print *,"! vrfy4"
!  print *,vrfy4;
  result(44) = vrfy4
!  print *,"! vrfy5"
!  print *,vrfy5;
  result(45) = vrfy5
!  print *,"! vrfy6"
!  print *,vrfy6;
  result(46) = vrfy6
!  print *,"! vrfy7"
!  print *,vrfy7;
  result(47) = vrfy7
!  print *,"! vrfy8"
!  print *,vrfy8;
  result(48) = vrfy8
!  print *,"! vrfy9"
!  print *,vrfy9;
  result(49) = vrfy9
!  print *,"! vrfyarry1"
!  print *,vrfyarry1;
  result(50:53) = vrfyarry1
!  print *,"! vrfyarry2"
!  print *,vrfyarry2;
  result(54:55) = vrfyarry2

!   print *, char1
  result(56) = ichar(char1)
!   print *, char2
  result(57) =ichar(char2)
!   print *, chararry1(1)
!   print *, chararry1(2)
!   print *, chararry1(3)
!   print *, chararry1(4)
!   print *, chararry1(5)
!   print *, chararry1(6)
  result(58:63) =ichar(chararry1)

!check kind
result(64) = kind(index("xxyybcwwvv", "bc", KIND=2));
!print *, result(64)
result(65) = kind(index("xxyybcwwvv", "bc", KIND=4));
!print *, result(65)
result(66) = kind(index("xxyybcwwvv", "bc", KIND=8));
!print *, result(66)
result(67) = kind(scan("xxyybcwwvv", "bc", KIND=2))
!print *, result(67)
result(68) = kind(scan("xxyybcwwvv", "bc", KIND=4))
!print *, result(68)
result(69) = kind(scan("xxyybcwwvv", "bc", KIND=8))
!print *, result(69)
result(70) = kind(verify("bbbbc", "bc", KIND=2))
!print *, result(70)
result(71) = kind(verify("bbbbc", "bc", KIND=4))
!print *, result(71)
result(72) = kind(verify("bbbbc", "bc", KIND=8))
!print *, result(72)
result(73) = kind(len_trim("xxyybc    ",KIND=2))
!print *, result(73)
result(74) = kind(len_trim("xxyybc    ",KIND=4))
!print *, result(74)
result(75) = kind(len_trim("xxyybc    ",KIND=8))
!print *, result(75)
result(76) = kind(char(65, KIND=1))
!print *, result(76)

  call check(result,expect,NTEST)

end program

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
	SUBROUTINE SN_GHGT ( isnfln, nparms, dattim, rlat, rlon, tmpk,
     +			     mxpt, maxdst, tblflg, mstflg, clmflg,
     +			     stid, dist, pres, hght, npt, iret )
C************************************************************************
C* SN_GHGT								*
C*									*
C* This subroutine finds the pressure and height of temperature, TMPK,	*
C* from sounding data.  If TBLFLG is .TRUE., a sounding table is used;	*
C* if MSTFLG is .TRUE., a moist-adiabatic cloud height is computed;	*
C* if CLMFLG is .TRUE., a climatology table is used to compute the	*
C* cloud heights.  If cloud heights cannot be found by stardard		*
C* sounding search, the moist-adiabatical method is employed. If no	*
C* cloud height can be found, it will go back 12h of the sounding data	*
C* and find the cloud height.						*
C*									*
C* SN_GHGT  ( ISNFLN, NPARMS, DATTIM, RLAT, RLON, TMPK, MXPT, MAXDST, 	*
C*	      TBLFLG, MSTFLG, CLMFLG, STID, DIST, PRES, HGHT, NPT,	*
C*	      IRET )							*
C*									*
C* Input parameters:							*
C*	ISNFLN		INTEGER		Sounding file no		*
C*	NPARMS		INTEGER		Number of parameters		*
C*	DATTIM		CHAR*		GEMPAK data/time		*
C*	RLAT		REAL		Latitude of the input point	*
C*	RLON		REAL		Longitude of the input point	*
C*	TMPK		REAL		Temperature in K		*
C*	MXPT		INTEGER		Maximum intersection points 	*
C*	MAXDST		INTEGER		Maximum search distance	(meters)*
C*	TBLFLG		LOGICAL		Sounding table flag		*
C*	MSTFLG		LOGICAL		Moist-adiabatic cloud flag	*
C*	CLMFLG		LOGICAL		Climatology table flag		*
C*									*
C* Output parameters:							*
C*	STID		CHAR*		Station closest the the point	*
C*	DIST		REAL		Dist between the pt and the stn	*
C*	PRES (NPT)	REAL		Pressure in mb			*
C*	HGHT (NPT)	REAL		Height in meters		*
C*	NPT		INTEGER		No of intersections		*
C*	IRET		INTEGER		Return code			*
C*					+10 = climatology data is used	*
C*					 +8 = moist-adiabatic height	*
C*					 +3 = too many intersections	*
C*					  0 = normal return		*
C*					 -4 = file not open		*
C*					-30 = cannot compute moist-	*
C*					      adiabatic cloud height	*
C**									*
C* Log:									*
C* T. Lee/GSC		 8/99	Created					*
C* T. Lee/GSC		 9/99	Checked no of intersections		*
C* T. Lee/GSC		11/99	Added max search distance and MAXSTN	*
C* T. Lee/GSC		12/99	Added sounding table flag		*
C* T. Lee/GSC		 1/00	Added moist-adiabatic cloud height	*
C* T. Lee/GSC		 4/00	Modulized moist-adiabatic cloud height	*
C* T. Lee/GSC		 5/00	Added climatology table			*
C* J. Wu/GSC            07/00   Moved INCLUDE 'ERMISS.FNC' before the   *
C*                              DATA statement                          *
C* T. Lee/GSC		 4/01	Used station number for station w/o ID	*
C* T. Lee/SAIC		 4/02	Fixed uninitialized memory read		*
C* T. Lee/SAIC		 4/03	Used T-12h data instead of climatology	*
C* T. Lee/SAIC		 5/03	Initialized variables			*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	PARAMETER	( JPRE=1, JTMP=2, JDWP=3, JHGT=6 )
	PARAMETER	( MAXSTN=100, NPMS=7, P600=600., P500=500. )
	CHARACTER*(*)	dattim, stid
	REAL		pres (*), hght (*)
C*
	CHARACTER	stn (LLSTFL)*8, stntbl*20, stnam*8, savstn*8
	CHARACTER	tbfil*72, parms(MMPARM)*4, cstn*8
	INTEGER		istn (LLSTFL)
	REAL		dd (LLSTFL), data (LLMXDT),
     +			datout (LLMXLV, MMPARM)
	REAL		datlev (MMPARM), outdat (MMPARM),
     +			adata  (MMPARM), bdata  (MMPARM)
	LOGICAL		proces, found, intflg (7), angflg (7)
	LOGICAL		tblflg, mstflg, clmflg, savflg
	INCLUDE		'ERMISS.FNC'
	DATA		intflg / 7*.true. /, angflg / 7*.false./
C*
C------------------------------------------------------------------------
C
C*	Initialize the output.
C
	stid = ' '
	dist = RMISSD
	DO i = 1, mxpt
	    pres ( i ) = RMISSD
	    hght ( i ) = RMISSD
	END DO
C
C*	Cannot compute moist-adiabatic cloud height with table data.
C
	IF  ( ( tblflg .or. clmflg ) .and. mstflg )  THEN
	    iret = -30
	    RETURN
	  ELSE
	    IF ( tblflg .or. clmflg )  THEN
		nstn = 1
		dist = 0.
		IF  ( tblflg )  THEN
		    stid = 'STANDARD'
		  ELSE
		    stid = 'CLIMATOLOGY'
		END IF
	    END IF
	END IF
C
C*	Check file.
C
	IF  ( tblflg .or. clmflg )  THEN
	    IF  ( isnfln .lt. 1 )  THEN
		iret = -4
		RETURN
	    END IF
	  ELSE
	    CALL SN_CHKF ( isnfln, iret )
	    IF  ( iret .ne. 0 )  RETURN
C
C*	    Compute the distance between the input point and the
C*	    sounding stations.
C
	    stntbl = 'snstns.tbl'
	    CALL LC_DIST  ( stntbl, rlat, rlon, MAXSTN, maxdst, stn,
     +			    istn, dd, nstn, iret )
C
C*	    Set the time.
C
	    CALL SN_STIM  ( isnfln, dattim, iret )
	END IF
C
C*	Convert temperature to Celsius.
C
 	tmpc = PR_TMKC ( tmpk )
C
C*	Initialize buffer for the nearest station.
C
	savstn = ' '
	savflg = .false.
	savdst = RMISSD
	tsave  = RMISSD
C
C*	Read the data.
C
	found  = .false.
	proces = .false.
	k    = 1
	npt  = 0
C
	DO WHILE  ( ( .not. found ) .and. ( k .le. nstn ) )
C
C*	    Read data from standard sounding table.
C
	    IF  ( tblflg )  THEN
		CALL SN_RTBL ( isnfln, nparms, data, nlev, iret )
		IF  ( iret .ne. 0 )  RETURN
		istnm = 999999
		selv  = 0.
		CALL PC_SSTN  ( stid, istnm, rlat, rlon, selv,
     +				ispri, ihhmm, nlev, iret )
C
C*	      Get cloud height from climatology.
C
	      ELSE IF  ( clmflg )  THEN
		CALL SN_GCLM  ( isnfln, nparms, dattim, rlat, tmpk,
     +				data, h, p, iret )
		IF  ( iret .eq. 0 ) THEN
		    npt = 1
		    pres ( npt ) = p
		    hght ( npt ) = h
		END IF
		RETURN
	      ELSE
C
C*	      Read data from a sounding file.
C
		stid   = stn  ( k )
		IF  ( stid .eq. ' ' )  THEN
		    CALL ST_INCH ( istn (k), cstn, ier )
		    stid = cstn
		END IF
		dist   = dd   ( k )
C
C*		Set station.
C
		CALL SN_SSTN  ( isnfln, stid, stnam, istnm, slat, slon,
     +				selv, iret )
		IF  ( iret .eq. 0 )  THEN
      		    CALL SN_RDAT  ( isnfln, nlev, data, ihhmm, iret )
      		    CALL PC_SSTN  ( stid, istnm, slat, slon, selv,
     +				    ispri, ihhmm, nlev, iret )
		END IF
	    END IF
C
C*	    Save the station information.
C
	    IF  ( iret .eq. 0 )  THEN
		proces = .true.
		IF  ( .not. savflg )  THEN
		    CALL PC_FTOP ( data, nparms, levtop, outdat, ier )
		    IF  ( .not. ERMISS ( outdat (JTMP)  ) )  THEN
			tsave  = outdat (JTMP)
			savflg = .true.
			savstn = stid
			savdst = dist
		    END IF
		END IF
	    END IF
C
C*	    Moist-adiabatic cloud height methodology: Find the most
C*	    unstable level below 600 mb. Lift the parcel to 500 mb.
C*	    If the lifted index is less than 0, (air parcel is warmer
C*	    than the environment), continue lifting the parcel
C*	    moist-adiabatically until the parcel temperature is cooled
C*	    off to the cloud temperature, TMPK.
C
	    IF  ( proces .and. mstflg )  THEN
		CALL SN_MGHT  ( data, nparms, tmpk, pres, hght, npt,
     +				iret )
		IF  ( iret .ge. 0 )  THEN
		    found = .true.
		    RETURN
		END IF
C
	      ELSE
C
		IF  ( proces .and. ( nlev .gt. 1 ) )  THEN
C
C*		    Store data below 20 mb only into the buffer.
C
		    DO i = 1, NPMS
			outdat ( i ) = RMISSD
		    END DO
C
		    jtop = 0
C
		    DO  num = 1, nlev
			CALL PC_GLEV  ( num, data, nparms,
     +					datlev, ier )
			CALL PC_COMP ( 5, datlev, outdat, ier )
			pp = outdat ( JPRE )
			tt = outdat ( JTMP )
			hh = outdat ( JHGT )
			IF  ( ERMISS ( pp ) .or. ERMISS ( tt ) .or.
     +			      ERMISS ( hh ) .or.
     +			    ( ( pp .lt. 20.) .and. ( .not. tblflg ) ) )
     +			    THEN
C
			  ELSE
			    jtop = jtop + 1
			    DO np = 1, NPMS
				datout ( jtop, np ) = outdat ( np )
			    END DO
			END IF
		    END DO
C
		    jbot = jtop
C
C*		    Search cloud temperature top down from the sounding
C*		    data or table.
C
		    tmax = RMISSD
		    DO WHILE ( jbot .ge. 2 )
			tmax = AMAX1 ( tmax, datout ( jbot, JTMP ) )
			jbot = jbot - 1
			sign =  ( tmpc - datout ( jtop, JTMP ) ) *
     +				( tmpc - datout ( jbot, JTMP ) )
C
			IF  ( sign .le. 0. )  THEN
			    found = .true.
C
C*			    Check if SIGN .eq. 0.  Else, interpolation
C*			    with respect to temperature.
C
			    IF  ( sign .eq. 0. )  THEN
				IF ( tmpc .eq. datout(jtop, JTMP) ) THEN
				    IF  ( npt .ge. mxpt )  THEN
					iret = +3
					RETURN
				      ELSE
					npt = npt + 1
					pres (npt) = datout (jtop, JPRE)
					hght (npt) = datout (jtop, JHGT)
				    END IF
				END IF
C
			      ELSE
				jb = jbot
				jt = jb + 1
				DO i = 1, NPMS
				    adata ( i ) = datout ( jt, i )
				    bdata ( i ) = datout ( jb, i )
				END DO
				CALL PC_INTT  ( tmpc, adata, bdata,
     +						NPMS, intflg, angflg,
     +						JTMP, outdat, iret )
				IF  ( iret .eq. 0 )  THEN
				    IF  ( npt .ge. mxpt )  THEN
					iret = +3
					RETURN
				      ELSE
					npt = npt + 1
					pres ( npt ) = outdat ( JPRE )
					hght ( npt ) = outdat ( JHGT )
				    END IF
				END IF
			    END IF
C
C*			    Reassign top level data.
C
			    jtop = jbot
			    DO i = 1, NPMS
		       		datout ( jtop, i ) = datout ( jbot, i )
			    END DO
C
			END IF
C
		    END DO
C
C*		    Check the lowest level data.
C
		    tmax = AMAX1 ( tmax, datout ( 1, JTMP ) )
		    IF  ( proces .and. ( tmpc .eq. datout( 1, JTMP ) ) )
     +			  THEN
			found = .true.
			IF  ( npt .ge. mxpt )  THEN
			    iret = +3
			    RETURN
		          ELSE
			    npt   = npt + 1
			    pres ( npt ) = datout ( jbot, JPRE )
			    hght ( npt ) = datout ( jbot, JHGT )
			END IF
		    END IF
C
C*		    Per AWC's request, if cloud height is not found,
C*		    use moist adiabatic method to find cloud height.
C*		    If found, send a return code for status message.
C
		    IF  ( ( .not. found ) .and. ( .not. tblflg ) )  THEN
			CALL SN_MGHT  ( data, nparms, tmpk, pres,
     +					hght, npt, iret )
			IF  ( iret .ge. 0 )  THEN
			    found = .true.
			    iret  = +8
			    RETURN
			END IF
		    END IF
C
		END IF
	    END IF
C
	    IF  ( .not. found )  k = k + 1
C
	END DO
C
C*	If no station is found, climatological data table will be used for
C*	standard atmospheric search.  For nearest sounding search an error
C*	code will be returned.
C
	IF  ( .not. found )  THEN
	    IF   ( tblflg )  THEN
		tbfil = 'climate.tbl'
		CALL SN_OPNT ( tbfil, lun, parms, nparms, iret )
		IF  ( ( iret .ne. 0 ) .or. ( nparms .lt. 3 ) )  THEN
		    RETURN
		  ELSE
		    CALL SN_GCLM  ( lun, nparms, dattim, rlat, tmpk,
     +				    data, h, p, iret )
		END IF
C
C*		Write status message for climatology.
C
		IF  ( .not. savflg )  THEN
		    stid = 'CLIMATOLOGY'
		    dist = 0.
		  ELSE
		    stid = savstn
		    dist = savdst
		END IF
C
		IF  ( iret .eq. 0 )  THEN
		    npt = 1
		    pres ( npt ) = p
		    hght ( npt ) = h
		    iret = +10
		  ELSE
		    npt = 0
		END IF
		CALL FL_CLOS ( lun, ier )
	      ELSE
		stid = savstn
		dist = savdst
		npt = 0
C
C*		Check if cloud temperature is warmer than the entire sounding
C*		data.
C
		IF  ( .not. ERMISS (tmax) .and. .not. ERMISS ( dist )
     +		    .and.  (tmpc .gt. tmax) .and. ( nstn .gt. 0 ) ) THEN
		    iret = + 6
		  ELSE
		    iret = -33
		END IF
	    END IF
	END IF
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
C$Procedure FGROUP_1 ( First simple group of a string )

      SUBROUTINE FGROUP_1 ( STRING, LGR, RGR, BEG, END, DEPTH )
      IMPLICIT NONE

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

      CHARACTER*(*)         STRING
      CHARACTER*(*)         LGR
      CHARACTER*(*)         RGR
      INTEGER               BEG
      INTEGER               END
      INTEGER               DEPTH

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     STRING     I   Any character string
C     LGR        I   A string that marks the start of a group
C     RGR        I   A string that marks the end   of a group
C     BEG        O   Index of the left  character of first full group
C     END        O   Index of the right character of first full group
C     DEPTH      O   Nesting depth of the first full group
C
C$ Detailed_Input
C
C     STRING     Any character string.
C
C     LGR        A character or string of characters that marks the
C                start of a group. Typically this will be a left
C                parenthesis, left square bracket, or left brace.
C
C     RGR        A character or string of characters that marks the
C                end of a group. Usually this will be a right
C                parenthesis, right square bracket, or right brace.
C
C$ Detailed_Output
C
C     BEG        Index of the left character of the first simple group.
C
C     END        Index of the last character of the first simple group.
C
C     DEPTH      Integer giving the nesting depth of the simple group
C                delimited by BEG and END.
C
C$ Files
C
C     None.
C
C$ Exceptions
C
C     1) If the input string is not balanced (unpaired groupers),
C        prior to the detection of the first simple group,
C        the error 'SPICE(UNBALANCEDGROUP)' will be signalled.
C        The values of BEG, END and DEPTH will be set to zero.
C
C     2) If no left or right groupers occur in the string, BEG and END
C        will point to the first and last characters of the string
C        respectively. (The last character has index equal to the
C        declared length of the string.)  DEPTH will be set to zero.
C
C     3) If either grouper is an introductory substring of the other
C        the error 'SPICE(NONDISTINCTPAIR)' will be signalled.
C        The values of BEG, END and DEPTH will be set to zero.
C
C
C$ Particulars
C
C     One symbol begins a group, another distinct symbol ends it.
C
C     Associated with each character of the string is a depth number.
C     The depth of a character is determined via the following
C     proceedure.
C
C        Initially assume that each character has depth zero.
C
C        Read the string from left to right.
C
C        If a left grouper is encountered
C
C           the depth is incremented by one and all characters of the
C           left grouper are assigned this depth.
C
C        If a right grouper is encountered
C
C           the depth is decremented by one and all characters of the
C           right grouper are assigned this depth.
C
C        Otherwise
C
C           The character is given a depth equal to the last depth.
C
C     For example, suppose that the groupers are '(' and ')'.  Then
C     the depths of each character of the string below is given
C     immediately beneath it.
C
C              1 + 2 * ( (3+3)*4 +(2+(3*4 + 3) ) ) + 27
C              00000000                          000000
C                      11    11111             11
C                        2222     222        22
C                                    33333333
C
C     The group of depth zero is simply the entire string. A group of
C     depth n begins at left grouper of depth n and ends on the
C     first right grouper of depth n-1. Taking the above example
C     again the groups of depth 2 are underlined.
C
C              1 + 2 * ( (3+3)*4 +(2+(3*4 + 3) ) ) + 27
C                        -----    --------------
C
C     A simple group, is a group that contains no groups of greater
C     depth, i.e. a substring that begins with a left grouper,
C     ends with a right grouper, and contains neither type of
C     grouper between them.  The simple groups of our example
C     are underlined below.
C
C              1 + 2 * ( (3+3)*4 +(2+(3*4 + 3) ) ) + 27
C                        -----       ---------
C
C     There are three routines available for locating the
C     groups of a string.  They are:
C
C        FGROUP --- Find the first simple group of a string.
C        DGROUP --- Find the deepest group of a string.
C        NGROUP --- Find the kth group of depth n.
C
C     Related routines are
C
C        SGROUP --- Resolve super groupers
C        BGROUP --- Determine how balanced a group is.
C
C     This routine finds the endpoints of the first simple group of
C     a string and its depth.
C
C$ Examples
C
C     The tables below lists a collection of sample input strings.
C     The substring STRING(BEG:END) is underlined.  The depth that
C     would be return for the nesting level is listed under the column
C     DEPTH.
C
C     E X A M P L E     I N P U T S    A N D    O U T P U T S
C
C     This first example illustrates how the routine works when
C     standard grouping symbols are used.
C
C           LGR     = '('
C           RGR     = ')'
C
C           STRING                                    DEPTH
C           ==============================            =====
C
C          '( ( 2 + 4 ) * 7 ) + 19        '             2
C             ---------
C
C          '.NOT. ( ( A .OR B ) .AND. C ) '             2
C                   -----------
C
C          '( ( (( ((( X + Y )            '             7
C                    ---------
C
C          'THIS HAS NO DELIMITERS        '             0
C           ------------------------------
C
C          'THE MAN ( JIM ) WENT HOME     '             1
C                   -------
C
C          '( 12*[3 + 5] )                '             1
C           --------------
C
C
C
C     The following example illustrate the need for care when choosing
C     left and right groupers of more than one character.
C
C
C           LGR    = '\beg'
C           RGR    = '\end'
C
C           STRING                                    DEPTH
C           ==============================            =====
C
C          'Data: \begin 1 2  3 4 \enddata'              1
C                 --------------------
C
C          '\begin \beg time \endit \end  '              2
C                  --------------
C
C
C
C     Finally the next example shows that you need to be careful
C     about passing trailing blanks attached to the right and left
C     groupers.
C
C
C           LGR    = ' [ '
C           RGR    = '] '
C
C           STRING                                    DEPTH
C           ==============================            =====
C
C          'A + [B+7]*32                  '             0
C           ------------------------------
C
C          'A + [ B + 7 ]*32              '       results in error
C                                                 SPICE(UNBALANCEDGROUP)
C
C          'A + [ B + 7] *32              '             1
C              ----------
C
C
C     ===============================================================
C
C
C     A    C O D E     E X A M P L E
C
C
C     The following loop shows how one might use this routine to
C     parse an arithmetic expression.
C
C        See if the string is BALANCED
C
C        DO WHILE ( BGROUP( STRING, LGR, RGR, SGR ) .NE. 0 )
C
C           Attempt to balance the string.
C
C        END DO
C
C        CALL FGROUP ( STRING, '(', ')', BEG, END, LEVEL )
C
C        DO WHILE ( LEVEL .GT. 0 )
C
C           Simplyfy the expression from BEG to END
C
C           Insert the simplified expression in place of the previous
C           one in STRING.
C
C           CALL FGROUP ( STRING, '(', ')', BEG, END, LEVEL )
C
C        END DO
C
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     W.L. Taber     (JPL)
C     I.M. Underwood (JPL)
C
C$ Version
C
C     Beta Version 1.0.0, 24-OCT-1988 (WLT) (IMU)
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               RETURN

C
C     Local variables
C

      INTEGER               LENGTH
      INTEGER               LGLEN
      INTEGER               RGLEN
      INTEGER               LGOFF
      INTEGER               RGOFF

      INTEGER               NEST
      INTEGER               I

      INTEGER               R
      INTEGER               L
      INTEGER               MINLEN

C%&END_DECLARATIONS

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'FGROUP_1' )
      END IF

C
C     Get the lengths of the left and right groupers as well as
C     the input string.
C

      LENGTH = LEN ( STRING )
      LGLEN  = LEN ( LGR    )
      RGLEN  = LEN ( RGR    )

      LGOFF  = LGLEN  - 1
      RGOFF  = RGLEN  - 1

      MINLEN = MIN ( LGLEN, RGLEN )

C
C     Check for dumb errors
C
      IF (      ( RGR    .EQ. LGR(1:MINLEN) )
     .     .OR. ( LGR    .EQ. RGR(1:MINLEN) ) ) THEN

         CALL SIGERR ( 'SPICE(NONDISTINCTPAIR)' )
         CALL CHKOUT ( 'FGROUP_1'               )
         RETURN

      END IF

C
C     Now we have the expected case.  Initialize the items to be used
C     for the search loop.
C
      NEST   =  0
      BEG    =  1
      END    =  LENGTH
      I      =  1

      DO WHILE ( I .LE. LENGTH )

C
C        Figure out where the ends of the groupers might appear in
C        the input string.
C
         R = I + RGOFF
         L = I + LGOFF

         R = MIN ( R, LENGTH )
         L = MIN ( L, LENGTH )

C
C        If this is a right grouper, we are done one way or the other.
C
         IF ( STRING(I:R) .EQ. RGR ) THEN

            IF ( NEST .GT. 0 ) THEN
               END   = R
               DEPTH = NEST
            ELSE
               BEG   = 0
               END   = 0
               DEPTH = 0

               CALL SIGERR ( 'SPICE(UNBALANCEDGROUP)' )

            END IF

            CALL CHKOUT ( 'FGROUP_1' )
            RETURN

         ELSE IF ( STRING(I:L) .EQ. LGR ) THEN

            NEST   =  NEST + 1
            BEG    =  I
            I      =  L + 1

         ELSE

            I     = I + 1

         END IF

      END DO

C
C     If the routine makes it this far, we had better have a string
C     with no groupers, i.e. current nesting zero.  Otherwise we
C     have an error.
C
      IF ( NEST .NE. 0 ) THEN

         BEG   = 0
         END   = 0
         DEPTH = 0

         CALL SIGERR ( 'SPICE(UNBALANCEDGROUP)' )

      ELSE

         DEPTH = 0
         BEG   = 1
         END   = LENGTH

      END IF

      CALL CHKOUT ( 'FGROUP_1' )
      RETURN

      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text: The Davis Jaycees was the local branch of the  http://en.wikipedia.org/wiki/United_States_Junior_Chamber United States Junior Chamber. The branch was shut down due to one of the most shocking incidents in Davis history. On November 19, 1983, police raided the Jaycees clubhouse, which was located between Davis and Woodland. Undercover officers discovered illegal sex acts, illegal gambling, illegal drug and alcohol sales, and the auctioning of a night with a prostitute. A total of 22 arrests were made.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c-----------------------------------------------------------------------
      subroutine setquad
C
C     Set up the quadrant points and associated element-element links
C     for plotting r-s-t surfaces in 3D and for establishing general
C     zipper planes.        pff 11-20-90  21:30 and  1-8-91
C
C     (Note: this algorithm is inherently CONFORMING, as constant r-s-t
C            planes have little significance in a non-conforming
C            discretization where the continuum assumption is not honored.)
C
#     include "basics.inc"
      include 'basicsp.inc'
C
C     Algorithm:  Set up indices, e.g:
C                                                     ^ s
C         These points define planes                  |
C         which propagate from one         +----------|---------+
C         element to the next, as          |          |         |
C         shown near pt. 3:                |         4+         |
C                                          |          |         |
C                                          |          |         |
C                                          +--1+------------+2--+--> r
C                                          |          |         |
C                                          |          |         |
C                                          | - - - - 3+ - - - - |
C                                          |          |         |
C                                          +--------------------+
C
C     Set up the connectivity list by performing the K^2 loop which
C     searches for coincident nodes, thus identifying an element-element
C     interaction.
C
CXXX
CXXX   recall, that normals will always be defined positive if they point
CXXX   away from the element center...,,, thus, once the topology of the
CXXX   normal (in or out) is established for one element, it will carry through
CXXX   to the rest of the elements.
CXXX
CXXX   2)  we can have a pseudo XYZ plane for now which takes as the default
CXXX   normal to be x y or z at the given point.
CXXX
CXXX   3)  when determining which  element contains the originating point, just
CXXX       use the first found point, and choose the closest quadpoint with an
CXXX       appropriately oriented normal vector (in some distance/normal norm).
CXXX
C
      DIMENSION XYZ(3,8)
      DIMENSION INEIGH(6,6)
      DATA      INEIGH/  1, -1, -1,  1,  1, -1
     $                , -1,  1,  1, -1, -1,  1
     $                , -1,  1,  1, -1, -1,  1
     $                ,  1, -1, -1,  1,  1, -1
     $                ,  1, -1, -1,  1,  1, -1
     $                , -1,  1,  1, -1, -1,  1 /
C
      IF (.NOT.IF3D) return
      CALL GENCEN
      NFACES=6
      EPSM = 1.0E-10
      NEL5=120*NELM
      CALL RZERO(XYZQ,NEL5)

c     if (nel.gt.20000) then
c        call prsi('Will not set quad for nel > 20K. NEL=$',nel)
c        return
c     endif

C     Define quadpoints:
      DO 100 IE=1,NEL
         CALL MAPXC(XYZ,IE)
C
C     Plane 1
         XYZQ(1,1,1,IE)=.125*(3.0*(XYZ(1,1)+XYZ(1,5))+XYZ(1,2)+XYZ(1,6))
         XYZQ(2,1,1,IE)=.125*(3.0*(XYZ(2,1)+XYZ(2,5))+XYZ(2,2)+XYZ(2,6))
         XYZQ(3,1,1,IE)=.125*(3.0*(XYZ(3,1)+XYZ(3,5))+XYZ(3,2)+XYZ(3,6))
C
         XYZQ(1,2,1,IE)=.125*(3.0*(XYZ(1,3)+XYZ(1,7))+XYZ(1,4)+XYZ(1,8))
         XYZQ(2,2,1,IE)=.125*(3.0*(XYZ(2,3)+XYZ(2,7))+XYZ(2,4)+XYZ(2,8))
         XYZQ(3,2,1,IE)=.125*(3.0*(XYZ(3,3)+XYZ(3,7))+XYZ(3,4)+XYZ(3,8))
C
         XYZQ(1,3,1,IE)=.125*(3.0*(XYZ(1,1)+XYZ(1,3))+XYZ(1,2)+XYZ(1,4))
         XYZQ(2,3,1,IE)=.125*(3.0*(XYZ(2,1)+XYZ(2,3))+XYZ(2,2)+XYZ(2,4))
         XYZQ(3,3,1,IE)=.125*(3.0*(XYZ(3,1)+XYZ(3,3))+XYZ(3,2)+XYZ(3,4))
C
         XYZQ(1,4,1,IE)=.125*(3.0*(XYZ(1,5)+XYZ(1,7))+XYZ(1,6)+XYZ(1,8))
         XYZQ(2,4,1,IE)=.125*(3.0*(XYZ(2,5)+XYZ(2,7))+XYZ(2,6)+XYZ(2,8))
         XYZQ(3,4,1,IE)=.125*(3.0*(XYZ(3,5)+XYZ(3,7))+XYZ(3,6)+XYZ(3,8))
C
C     Plane 2
         XYZQ(1,1,2,IE)=.125*(3.0*(XYZ(1,2)+XYZ(1,6))+XYZ(1,1)+XYZ(1,5))
         XYZQ(2,1,2,IE)=.125*(3.0*(XYZ(2,2)+XYZ(2,6))+XYZ(2,1)+XYZ(2,5))
         XYZQ(3,1,2,IE)=.125*(3.0*(XYZ(3,2)+XYZ(3,6))+XYZ(3,1)+XYZ(3,5))
C
         XYZQ(1,2,2,IE)=.125*(3.0*(XYZ(1,4)+XYZ(1,8))+XYZ(1,3)+XYZ(1,7))
         XYZQ(2,2,2,IE)=.125*(3.0*(XYZ(2,4)+XYZ(2,8))+XYZ(2,3)+XYZ(2,7))
         XYZQ(3,2,2,IE)=.125*(3.0*(XYZ(3,4)+XYZ(3,8))+XYZ(3,3)+XYZ(3,7))
C
         XYZQ(1,3,2,IE)=.125*(3.0*(XYZ(1,2)+XYZ(1,4))+XYZ(1,1)+XYZ(1,3))
         XYZQ(2,3,2,IE)=.125*(3.0*(XYZ(2,2)+XYZ(2,4))+XYZ(2,1)+XYZ(2,3))
         XYZQ(3,3,2,IE)=.125*(3.0*(XYZ(3,2)+XYZ(3,4))+XYZ(3,1)+XYZ(3,3))
C
         XYZQ(1,4,2,IE)=.125*(3.0*(XYZ(1,6)+XYZ(1,8))+XYZ(1,5)+XYZ(1,7))
         XYZQ(2,4,2,IE)=.125*(3.0*(XYZ(2,6)+XYZ(2,8))+XYZ(2,5)+XYZ(2,7))
         XYZQ(3,4,2,IE)=.125*(3.0*(XYZ(3,6)+XYZ(3,8))+XYZ(3,5)+XYZ(3,7))
C
C     Plane 5
         XYZQ(1,1,5,IE)=.125*(3.0*(XYZ(1,1)+XYZ(1,3))+XYZ(1,5)+XYZ(1,7))
         XYZQ(2,1,5,IE)=.125*(3.0*(XYZ(2,1)+XYZ(2,3))+XYZ(2,5)+XYZ(2,7))
         XYZQ(3,1,5,IE)=.125*(3.0*(XYZ(3,1)+XYZ(3,3))+XYZ(3,5)+XYZ(3,7))
C
         XYZQ(1,2,5,IE)=.125*(3.0*(XYZ(1,2)+XYZ(1,4))+XYZ(1,6)+XYZ(1,8))
         XYZQ(2,2,5,IE)=.125*(3.0*(XYZ(2,2)+XYZ(2,4))+XYZ(2,6)+XYZ(2,8))
         XYZQ(3,2,5,IE)=.125*(3.0*(XYZ(3,2)+XYZ(3,4))+XYZ(3,6)+XYZ(3,8))
C
         XYZQ(1,3,5,IE)=.125*(3.0*(XYZ(1,1)+XYZ(1,2))+XYZ(1,5)+XYZ(1,6))
         XYZQ(2,3,5,IE)=.125*(3.0*(XYZ(2,1)+XYZ(2,2))+XYZ(2,5)+XYZ(2,6))
         XYZQ(3,3,5,IE)=.125*(3.0*(XYZ(3,1)+XYZ(3,2))+XYZ(3,5)+XYZ(3,6))
C
         XYZQ(1,4,5,IE)=.125*(3.0*(XYZ(1,3)+XYZ(1,4))+XYZ(1,7)+XYZ(1,8))
         XYZQ(2,4,5,IE)=.125*(3.0*(XYZ(2,3)+XYZ(2,4))+XYZ(2,7)+XYZ(2,8))
         XYZQ(3,4,5,IE)=.125*(3.0*(XYZ(3,3)+XYZ(3,4))+XYZ(3,7)+XYZ(3,8))
C
C     Plane 6
         XYZQ(1,1,6,IE)=.125*(3.0*(XYZ(1,5)+XYZ(1,7))+XYZ(1,1)+XYZ(1,3))
         XYZQ(2,1,6,IE)=.125*(3.0*(XYZ(2,5)+XYZ(2,7))+XYZ(2,1)+XYZ(2,3))
         XYZQ(3,1,6,IE)=.125*(3.0*(XYZ(3,5)+XYZ(3,7))+XYZ(3,1)+XYZ(3,3))
C
         XYZQ(1,2,6,IE)=.125*(3.0*(XYZ(1,6)+XYZ(1,8))+XYZ(1,2)+XYZ(1,4))
         XYZQ(2,2,6,IE)=.125*(3.0*(XYZ(2,6)+XYZ(2,8))+XYZ(2,2)+XYZ(2,4))
         XYZQ(3,2,6,IE)=.125*(3.0*(XYZ(3,6)+XYZ(3,8))+XYZ(3,2)+XYZ(3,4))
C
         XYZQ(1,3,6,IE)=.125*(3.0*(XYZ(1,5)+XYZ(1,6))+XYZ(1,1)+XYZ(1,2))
         XYZQ(2,3,6,IE)=.125*(3.0*(XYZ(2,5)+XYZ(2,6))+XYZ(2,1)+XYZ(2,2))
         XYZQ(3,3,6,IE)=.125*(3.0*(XYZ(3,5)+XYZ(3,6))+XYZ(3,1)+XYZ(3,2))
C
         XYZQ(1,4,6,IE)=.125*(3.0*(XYZ(1,7)+XYZ(1,8))+XYZ(1,3)+XYZ(1,4))
         XYZQ(2,4,6,IE)=.125*(3.0*(XYZ(2,7)+XYZ(2,8))+XYZ(2,3)+XYZ(2,4))
         XYZQ(3,4,6,IE)=.125*(3.0*(XYZ(3,7)+XYZ(3,8))+XYZ(3,3)+XYZ(3,4))
C
C     Plane 3
         XYZQ(1,1,3,IE)=.125*(3.0*(XYZ(1,1)+XYZ(1,5))+XYZ(1,3)+XYZ(1,7))
         XYZQ(2,1,3,IE)=.125*(3.0*(XYZ(2,1)+XYZ(2,5))+XYZ(2,3)+XYZ(2,7))
         XYZQ(3,1,3,IE)=.125*(3.0*(XYZ(3,1)+XYZ(3,5))+XYZ(3,3)+XYZ(3,7))
C
         XYZQ(1,2,3,IE)=.125*(3.0*(XYZ(1,2)+XYZ(1,6))+XYZ(1,4)+XYZ(1,8))
         XYZQ(2,2,3,IE)=.125*(3.0*(XYZ(2,2)+XYZ(2,6))+XYZ(2,4)+XYZ(2,8))
         XYZQ(3,2,3,IE)=.125*(3.0*(XYZ(3,2)+XYZ(3,6))+XYZ(3,4)+XYZ(3,8))
C
         XYZQ(1,3,3,IE)=.125*(3.0*(XYZ(1,1)+XYZ(1,2))+XYZ(1,3)+XYZ(1,4))
         XYZQ(2,3,3,IE)=.125*(3.0*(XYZ(2,1)+XYZ(2,2))+XYZ(2,3)+XYZ(2,4))
         XYZQ(3,3,3,IE)=.125*(3.0*(XYZ(3,1)+XYZ(3,2))+XYZ(3,3)+XYZ(3,4))
C
         XYZQ(1,4,3,IE)=.125*(3.0*(XYZ(1,5)+XYZ(1,6))+XYZ(1,7)+XYZ(1,8))
         XYZQ(2,4,3,IE)=.125*(3.0*(XYZ(2,5)+XYZ(2,6))+XYZ(2,7)+XYZ(2,8))
         XYZQ(3,4,3,IE)=.125*(3.0*(XYZ(3,5)+XYZ(3,6))+XYZ(3,7)+XYZ(3,8))
C
C     Plane 4
         XYZQ(1,1,4,IE)=.125*(3.0*(XYZ(1,3)+XYZ(1,7))+XYZ(1,1)+XYZ(1,5))
         XYZQ(2,1,4,IE)=.125*(3.0*(XYZ(2,3)+XYZ(2,7))+XYZ(2,1)+XYZ(2,5))
         XYZQ(3,1,4,IE)=.125*(3.0*(XYZ(3,3)+XYZ(3,7))+XYZ(3,1)+XYZ(3,5))
C
         XYZQ(1,2,4,IE)=.125*(3.0*(XYZ(1,4)+XYZ(1,8))+XYZ(1,2)+XYZ(1,6))
         XYZQ(2,2,4,IE)=.125*(3.0*(XYZ(2,4)+XYZ(2,8))+XYZ(2,2)+XYZ(2,6))
         XYZQ(3,2,4,IE)=.125*(3.0*(XYZ(3,4)+XYZ(3,8))+XYZ(3,2)+XYZ(3,6))
C
         XYZQ(1,3,4,IE)=.125*(3.0*(XYZ(1,3)+XYZ(1,4))+XYZ(1,1)+XYZ(1,2))
         XYZQ(2,3,4,IE)=.125*(3.0*(XYZ(2,3)+XYZ(2,4))+XYZ(2,1)+XYZ(2,2))
         XYZQ(3,3,4,IE)=.125*(3.0*(XYZ(3,3)+XYZ(3,4))+XYZ(3,1)+XYZ(3,2))
C
         XYZQ(1,4,4,IE)=.125*(3.0*(XYZ(1,7)+XYZ(1,8))+XYZ(1,5)+XYZ(1,6))
         XYZQ(2,4,4,IE)=.125*(3.0*(XYZ(2,7)+XYZ(2,8))+XYZ(2,5)+XYZ(2,6))
         XYZQ(3,4,4,IE)=.125*(3.0*(XYZ(3,7)+XYZ(3,8))+XYZ(3,5)+XYZ(3,6))
C
  100 CONTINUE
C
C     Find Neighbor-Neighbor connections:
C
      DO 500 IE=1,NEL
      DO 500 IFACE=1,NFACES
      DO 500 IQPT=1,4
C
         NEIGHI=INT(XYZQ(4,IQPT,IFACE,IE))
         IF (NEIGHI.EQ.0) THEN
            IQP1=IQPT+1
            IQP1=MOD1(IQP1,4)
            EPSI=(XYZQ(1,IQPT,IFACE,IE)-XYZQ(1,IQP1,IFACE,IE))**2
     $          +(XYZQ(2,IQPT,IFACE,IE)-XYZQ(2,IQP1,IFACE,IE))**2
     $          +(XYZQ(3,IQPT,IFACE,IE)-XYZQ(3,IQP1,IFACE,IE))**2
C
c           DO 204 JE=1,NEL
c   pff: 4/13/93
            DO 204 JE=ie,Nel
            IF (JE.EQ.IE) GOTO 202
            DIST=(XCEN(IE)-XCEN(JE))**2+(YCEN(IE)-YCEN(JE))**2
     $          +(ZCEN(IE)-ZCEN(JE))**2
            DST2=(RCEN(IE)+RCEN(JE))**2
c           write(6,101) ie,je,dst2,dist,rcen(ie),rcen(je)
  101       format(2x,'dist:',2i3,4f12.7)
            IF (DST2.LT.DIST) GOTO 202
            DO 200 JFACE=1,NFACES
            DO 200 JQPT=1,4
c   pff: 4/13/93
            NEIGHJ=INT(XYZQ(4,IQPT,IFACE,IE))
            IF (NEIGHJ.ne.0) GOTO 200
c   pff: 4/13/93
               JQP1=JQPT+1
               JQP1=MOD1(JQP1,4)
               EPSJ=(XYZQ(1,JQPT,JFACE,JE)-XYZQ(1,JQP1,JFACE,JE))**2
     $             +(XYZQ(2,JQPT,JFACE,JE)-XYZQ(2,JQP1,JFACE,JE))**2
     $             +(XYZQ(3,JQPT,JFACE,JE)-XYZQ(3,JQP1,JFACE,JE))**2
               EPS=EPSM*(EPSJ+EPSI)
               EPS=MAX(EPS,EPSM)
               DIST=(XYZQ(1,IQPT,IFACE,IE)-XYZQ(1,JQPT,JFACE,JE))**2
     $             +(XYZQ(2,IQPT,IFACE,IE)-XYZQ(2,JQPT,JFACE,JE))**2
     $             +(XYZQ(3,IQPT,IFACE,IE)-XYZQ(3,JQPT,JFACE,JE))**2
               IF (DIST.LE.EPS) THEN
C                 We found our match, make a note of it
C                 and move to the next point.
                  NEIGHI=24*(JE-1)+4*(JFACE-1)+JQPT
                  NEIGHJ=24*(IE-1)+4*(IFACE-1)+IQPT
                  INRM=INEIGH(IFACE,JFACE)
                  XYZQ(4,NEIGHJ,1,1)=FLOAT(INRM*NEIGHI)
                  XYZQ(4,NEIGHI,1,1)=FLOAT(INRM*NEIGHJ)
                  GOTO 300
               endif
  200       CONTINUE
  202       CONTINUE
  204       CONTINUE
  300       CONTINUE
         endif
  500 CONTINUE
c
c     Initialize interp routine ---
c
      CALL INTERP(VAL1,xp(1),yp(1),zp(1),idum,wkv1,ierr)
C
C     Diagnostic (1-8-91)
C
c     CALL DRQUAD
C
      return
      end
c-----------------------------------------------------------------------
      subroutine drquad
#     include "basics.inc"
C
      NFACES=6
      DO 600 IE=1,NEL
      DO 600 IFACE=1,NFACES
      DO 600 IQPT=1,4
         ICLR=2*IFACE
c        ICLR=MOD1(ICLR,14)
c        ICLR=5
         CALL COLOR(ICLR)
         CALL XDIAM3(XYZQ(1,IQPT,IFACE,IE))
  600 CONTINUE
C
      return
      end
c-----------------------------------------------------------------------
      subroutine xdiam3(xyz)
c     draws a diamond around XX,YY
#     include "basics.inc"
C
      DIMENSION XYZ(3)
      XX=XISO(XYZ(1),XYZ(2),XYZ(3))
      YY=YISO(XYZ(1),XYZ(2),XYZ(3))
      CALL MOVEC(XX+.0025*XFAC ,YY+.0025*YFAC)
      CALL DRAWC(XX-.0025*XFAC ,YY+.0025*YFAC)
      CALL DRAWC(XX-.0025*XFAC ,YY-.0025*YFAC)
      CALL DRAWC(XX+.0025*XFAC ,YY-.0025*YFAC)
      CALL DRAWC(XX+.0025*XFAC ,YY+.0025*YFAC)
      return
      end
c-----------------------------------------------------------------------
      subroutine mapxc(xyz,ie)
C     Map the element corner points (x,y,z) to standard vector.
#     include "basics.inc"
      DIMENSION INDX(8)
      DIMENSION XYZ(3,8)
      INDX(1)=1
      INDX(2)=2
      INDX(3)=4
      INDX(4)=3
      INDX(5)=5
      INDX(6)=6
      INDX(7)=8
      INDX(8)=7
      NDIM2 = 2**NDIM
C
      DO 50 IX=1,NDIM2
         I=INDX(IX)
         XYZ(1,IX)=X(IE,I)
         XYZ(2,IX)=Y(IE,I)
         XYZ(3,IX)=Z(IE,I)
   50 CONTINUE
      return
      end
c-----------------------------------------------------------------------
      subroutine decod(ix,iy,iz,ie,Jxyz,nx,ny,nz)
C
      ixyz = abs(jxyz)
      NXY  = NX*NY
      NXYZ = NXY*NZ
C
      IE  = (IXYZ + NXYZ - 1)/NXYZ
      IXY = IXYZ-(IE-1)*NXYZ
      IZ  = (IXY + NXY - 1)/NXY
      IXX = IXY-(IZ-1)*NXY
      IY  = (IXX + NX - 1)/NX
      IX  = IXX-(IY-1)*NX
C
C     Check
C
      ICHECK = (IE-1)*NXYZ + (IZ-1)*NXY + (IY-1)*NX + IX
      IF (ICHECK.NE.IXYZ) THEN
         WRITE(6,*) 'hey, what gives here?  icheck,ixyz:',icheck,ixyz
      endif
      return
      end
c-----------------------------------------------------------------------
      subroutine setsrfp(srfqal,nqual,ifld)
      CHARACTER*1 SRFQAL(1)
      CHARACTER*1 CB1,CB
C
#     include "basics.inc"
C
C     Select a list of planes according to the specified boundary
C     qualifications given in SRFQAL.  If the number of qualifications,
C     NQUAL, is negative, select all boundaries NOT meeting the
C     qualifications.
C
C     Note:  Modification 10/2/92 -
C            Set INRM such that the "+" direction associated with the
C            plane has the same direction as the closest fundamental
C            unit vector (e_1,e_2, or e_3).    pff.
C
C
c     NLSTP=0                      pff 12/21/96
      NFACES=2*NDIM
C
      if (nqual.gt.0) then
c
         do iq=1,nqual
            cb1=srfqal(iq)
            if (cb1.eq.'1'.or.cb1.eq.'2'.or.cb1.eq.'3'
     $      .or.cb1.eq.'4'.or.cb1.eq.'5'.or.cb1.eq.'6') then
               read(cb1,1) iface
    1          format(i1)
               do ie=1,nel
                  nlstp=nlstp+1
                  ijkpln=eface1(iface)
                  ix=1
                  if (mod(ijkpln,2).eq.0) ix=nx
                  list(nlstp)=6*nx*(ie-1)+nx*(ijkpln-1)+ix
                  call signpl(list(nlstp))
               enddo
            else
               do ie=1,nel
               do iface=1,nfaces
                  cb=cbc(iface,ie,ifld)
                  if (cb.eq.cb1) then
                     nlstp=nlstp+1
                     ijkpln=eface1(iface)
                     ix=1
                     if (mod(ijkpln,2).eq.0) ix=nx
                     list(nlstp)=6*nx*(ie-1)+nx*(ijkpln-1)+ix
                     call signpl(list(nlstp))
                  endif
               enddo
               enddo
            endif
         enddo
C
      ELSE
C
c        open(unit=55,file='bdry')
         NQ=-NQUAL
         DO 300 IE=1,NEL
         DO 300 IFACE=1,NFACES
            CB=CBC(IFACE,IE,IFLD)
            DO 200 IQ=1,NQ
               CB1=SRFQAL(IQ)
C              if it's one of these, don't plot.
               IF (CB1.EQ.CB) GOTO 300
  200       CONTINUE
C
C           Since the b.c. isn't one of the one's specified, we plot
C           this surface....
            NLSTP=NLSTP+1
            IJKPLN=EFACE1(IFACE)
c           call outvtx(ie,iface,55)
            IX=1
            IF (MOD(IJKPLN,2).EQ.0) IX=NX
            LIST(NLSTP)=6*NX*(IE-1)+NX*(IJKPLN-1)+IX
            CALL SIGNPL(LIST(NLSTP))
  300    CONTINUE
c        close(unit=55)
C
      endif
      return
      end
c-----------------------------------------------------------------------
      subroutine sortl
#     include "basics.inc"
      include 'basicsp.inc'
      COMMON /CTMPS1/ IND(NELM,6),WK(NELM,6)
C
      NXY =NX*NY
      NXYZ=NZ*NXY
c
c     March over all mirror configurations, all planes
c
      i = 0
      do im  =1,nmirror
         do ipl =1,nlstp
            i = i+1
            LISTA=ABS(LIST(ipl))
            CALL DECOD(IPLANE,IPLN,IE,IDUM,LISTA,NX,6,NELM)
c
c           Z-sort based on centerpoint of selected plane
c
            IX=(NX+1)/2
            IY=(NY+1)/2
            IZ=(NZ+1)/2
            IF (IPLN.LE.2) THEN
               IX=IPLANE
            ELSEIF (IPLN.LE.4) THEN
               IY=IPLANE
            ELSE
               IZ=IPLANE
            endif
c
            IEOFF=NXYZ*(IE-1)+NXY*(IZ-1)+NX*(IY-1)+IX
            if (im.eq.1) then
               ZDEPTH(I,1)=ZISO(XP(IEOFF),YP(IEOFF),ZP(IEOFF))
            else
               call mirror(xmi,ymi,zmi,xp(ieoff),yp(ieoff),zp(ieoff),im)
               ZDEPTH(I,1)=ZISO(xmi,ymi,zmi)
            endif
            lmir(i) = im
         enddo
         call icopy(lisw(1+(im-1)*nlstp),list,nlstp)
      enddo
      CALL SORT(ZDEPTH,   IND,NLSTP*nmirror)
      CALL ISWAP(LISW ,WK,IND,NLSTP*nmirror)
      CALL ISWAP(lmir ,WK,IND,NLSTP*nmirror)
C
      return
      end
c-----------------------------------------------------------------------
      subroutine signpl(lplane)
C
C     Find the closest component of (1,1,1) to the normal associated
C     with this plane, and assign the sign of the dot product of that
C     component with the normal to Lplane.
C     11/92 pff.
C
#     include "basics.inc"
      include 'basicsp.inc'
      DIMENSION VEC0(3),VEC1(3),VEC2(3),VEC3(3)
C
      CALL DECOD(IPLANE,IPLN,IE,IDUM,Lplane,NX,6,NELM)
C
      I=(1+NX)/2
      J=(1+NY)/2
      K=(1+NZ)/2
C
      Nxyz=Nx*Ny*Nz
      Nxy =Nx*Ny
C
      IF (IPLN.LE.2) THEN
         I =IPLANE
         ieoff1 = nxyz*(ie-1)+nxy*(k-1)+nx*(j  )+i
         ieoff2 = nxyz*(ie-1)+nxy*(k  )+nx*(j-1)+i
      ELSEIF (IPLN.LE.4) THEN
         J =IPLANE
         ieoff1 = nxyz*(ie-1)+nxy*(k-1)+nx*(j-1)+i+1
         ieoff2 = nxyz*(ie-1)+nxy*(k  )+nx*(j-1)+i
      ELSE
         K =IPLANE
         ieoff1 = nxyz*(ie-1)+nxy*(k-1)+nx*(j-1)+i+1
         ieoff2 = nxyz*(ie-1)+nxy*(k-1)+nx*(j  )+i
      endif
      ieoff  = nxyz*(ie-1)+nxy*(k-1)+nx*(j-1)+i
C
      VEC0(1)=Xp(ieoff)
      VEC0(2)=Yp(ieoff)
      VEC0(3)=Zp(ieoff)
C
      VEC1(1)=Xp(ieoff1)
      VEC1(2)=Yp(ieoff1)
      VEC1(3)=Zp(ieoff1)
C
      VEC2(1)=Xp(ieoff2)
      VEC2(2)=Yp(ieoff2)
      VEC2(3)=Zp(ieoff2)
C
      CALL SUB3  (VEC1,VEC1,VEC0,3)
      CALL SUB3  (VEC2,VEC2,VEC0,3)
      CALL CROSS (VEC0,VEC1,VEC2)
      CALL NORM3D(VEC0)
C
      Vmax=0.0
      DO 10 I=1,3
         IF (ABS(VEC0(i)).GT.Vmax) THEN
             Vmax=ABS(VEC0(i))
             Imax=i
         endif
   10 CONTINUE
      IF (VEC0(Imax).LT.0.0) Lplane=-abs(Lplane)
c     write(6,*) ie,vec0(imax),lplane,' signpl'
      return
      end
c-----------------------------------------------------------------------
      subroutine gencen
#     include "basics.inc"
C
C     Generate the element centers
C
      IF (IF3D) THEN
         DO 100 IEL=1,NEL
            XCEN(IEL)=( X(IEL,1)+X(IEL,2)+X(IEL,3)+X(IEL,4)
     $               +  X(IEL,5)+X(IEL,6)+X(IEL,7)+X(IEL,8))/8.
            YCEN(IEL)=( Y(IEL,1)+Y(IEL,2)+Y(IEL,3)+Y(IEL,4)
     $               +  Y(IEL,5)+Y(IEL,6)+Y(IEL,7)+Y(IEL,8))/8.
            ZCEN(IEL)=( Z(IEL,1)+Z(IEL,2)+Z(IEL,3)+Z(IEL,4)
     $               +  Z(IEL,5)+Z(IEL,6)+Z(IEL,7)+Z(IEL,8))/8.
  100    CONTINUE
      ELSE
         DO 200 IEL=1,NEL
            XCEN(IEL)=(X(IEL,1)+X(IEL,2)+X(IEL,3)+X(IEL,4))/4.
            YCEN(IEL)=(Y(IEL,1)+Y(IEL,2)+Y(IEL,3)+Y(IEL,4))/4.
            ZCEN(IEL)=0.0
  200    CONTINUE
      endif
C
C     Compute the maximum radius from the center
C
      CALL RZERO(RCEN,NEL)
      IF (IF3D) THEN
         DO 300 IEL=1,NEL
         DO 300 J=1,8
            RAD=(X(IEL,J)-XCEN(IEL))**2 + (Y(IEL,J)-YCEN(IEL))**2
     $         +(Z(IEL,J)-ZCEN(IEL))**2
            RCEN(IEL)=MAX(RCEN(IEL),RAD)
  300    CONTINUE
      ELSE
         DO 400 IEL=1,NEL
         DO 400 J=1,4
            RAD=(X(IEL,J)-XCEN(IEL))**2 + (Y(IEL,J)-YCEN(IEL))**2
            RCEN(IEL)=MAX(RCEN(IEL),RAD)
  400    CONTINUE
      endif
      CALL VSQRT(RCEN,NEL)
C
      return
      end
c-----------------------------------------------------------------------
      subroutine clip_plane
C     Set up list of candidate RST planes
#     include "basics.inc"
      include 'basicsp.inc'
      COMMON /PFFLG/  ILGRNG,ISCND,IENTRP,INEWTX
      DIMENSION VEC(3),VEC1(3),VEC2(3),VEC3(3)
      CHARACTER*1 cplane,dplane
      LOGICAL IFANY,IFTMP
c
      IND(I,J,K,IEL)=I+NX*(J-1)+NX*NY*(K-1)+NX*NY*NZ*(IEL-1)
C
      CALL PRS('Clip in X, Y, or Z ?$')
      CALL RES(cplane,1)
      call capit(cplane,1)
      write(6,*)'cplane : ',cplane
c
      CALL PRS('Enter location of clipping plane:$')
      CALL RER(clipping_plane)
c
      CALL PRS
     $('Enter "<" or ">" if you wish to remove below or above plane.$')
      CALL RES(dplane,1)
c
C     3-D
C
      inew = 0
      DO 6000 I=1,NLSTP
C
C           New plotting algorithms for color fill and fishnet - pff 7-1-90
C           (3D will come later - 10-10-90)
C
         LISTA=ABS(LIST(I))
         CALL DECOD(IPLANE,IPLN,IE,IDUM,LISTA,NX,6,NELM)
         I1=1
         I2=NX
         J1=1
         J2=NY
         K1=1
         K2=NZ
         IF (IPLN.LE.2) THEN
            I1=IPLANE
            I2=IPLANE
         ELSEIF (IPLN.LE.4) THEN
            J1=IPLANE
            J2=IPLANE
         ELSE
            K1=IPLANE
            K2=IPLANE
         endif
c
         DO 5001 IZ=K1,K2
         DO 5001 IY=J1,J2
         DO 5001 IX=I1,I2
            IPOINT=IND(IX,IY,IZ,IE)
c
            if (cplane.eq.'X') qq = xp(ipoint)
            if (cplane.eq.'Y') qq = yp(ipoint)
            if (cplane.eq.'Z') qq = zp(ipoint)
c
c           See if we delete this plane from the list
            if (dplane.eq.'<'.and.qq.lt.clipping_plane) goto 5002
            if (dplane.eq.'>'.and.qq.gt.clipping_plane) goto 5002
c
 5001    CONTINUE
c        If we get here, then all points are included
         inew = inew+1
         list(inew) = list(i)
c
 5002    CONTINUE
 6000 CONTINUE
      nlstp = inew
      return
      end
c-----------------------------------------------------------------------
      subroutine mirror(xm,ym,zm,xi,yi,zi,m)
c
c     Based on mirror specification number, m,  compute the reflected
c     image of xi,yi,zi
c
c     For now, we support only one mirror -- that is the y=0 plane
c
      xm =  xi
      ym = -yi
      zm =  zi
c
      return
      end
c-----------------------------------------------------------------------
      subroutine vmirror(tm,xm,ym,zm,n,m)
      real tm(n,n),xm(n,n),ym(n,n),zm(n,n)
c
c     Based on mirror specification number, m,  compute the reflected
c     image of xi,yi,zi
c
c     For now, we support only one mirror -- that is the y=0 plane
c
      if (m.eq.2) then
         do i=1,n*n
            ym(i,1) = -ym(i,1)
         enddo
      endif
c
      return
      end
c-----------------------------------------------------------------------
      subroutine outvtx(ie,ifce,io)
c
#     include "basics.inc"
c
      integer fv3(4,6)
      save    fv3
      data    fv3 /  1,2,5,6
     $            ,  2,3,6,7
     $            ,  3,4,7,8
     $            ,  1,4,5,8
     $            ,  1,2,3,4
     $            ,  5,6,7,8  /
c
c     Output coordinates of this face
c
      xa = 0
      ya = 0
      za = 0
c
      do iv=1,4
         xa = xa+0.25*x(ie,fv3(iv,ifce))
         ya = ya+0.25*y(ie,fv3(iv,ifce))
         za = za+0.25*z(ie,fv3(iv,ifce))
      enddo
      write(io,1) xa,ya,za
    1 format(3f14.4)
c
      return
      end
c-----------------------------------------------------------------------
      subroutine assign_bc
c
c     prompt user for new bc to be applied to currently selected planes
c
#     include "basics.inc"
      include 'basicsp.inc'
c
      return
      end
c-----------------------------------------------------------------------
      subroutine gaujordf(a,m,n,indr,indc,ipiv,ierr)
C
C     Gauss-Jordan matrix inversion with full pivoting
c
c     Num. Rec. p. 30, 2nd Ed., Fortran
c
C
C     A is an nxn matrix, assumed contiguous in memory
C     W is a  work array of dimension n
C
      real    a(m,n)
      integer indr(m),indc(n),ipiv(n)
c
      real rmult
      integer wdsize
c
      ierr = 0
c
      eps = 1.e-12
      one = 1.
      wdsize = 8
      if (one+eps.eq.1.) wdsize = 4
c
      call izero(ipiv,n)
c
c
      do j=1,n
         amx=0.
c
c        Pivot search
         do i=1,m
            if (ipiv(i).ne.1) then
               do k=1,n
                  if (ipiv(k).eq.0) then
                    if (abs(a(i,k)).ge.amx) then
                       amx = abs(a(i,k))
                       ir  = i
                       jc  = k
                    endif
                 elseif (ipiv(k).gt.1) then
                    ierr = -ipiv(k)
                    return
                 endif
              enddo
           endif
        enddo
        ipiv(jc) = ipiv(jc) + 1
c
c       Swap rows
        if (ir.ne.jc) then
           do k=1,n
              tmp     = a(ir,k)
              a(ir,k) = a(jc,k)
              a(jc,k) = tmp
           enddo
        endif
        indr(j)=ir
        indc(j)=jc
        if (abs(a(jc,jc)).lt.eps) then
           write(6,*) 'small Gauss Jordan Piv:',jc,a(jc,jc)
           ierr = jc
           return
        endif
        piv = 1./a(jc,jc)
        a(jc,jc)=1.
        do k=1,n
           a(jc,k) = a(jc,k)*piv
        enddo
c
        do i=1,m
           if (i.ne.jc) then
              rmult   = a(i,jc)
              a(i,jc) = 0.
              do k=1,n
                 a(i,k) = a(i,k) - rmult*a(jc,k)
              enddo
           endif
        enddo
      enddo
c
c     Unscramble matrix
      do k=n,1,-1
         if (indr(k).ne.indc(k)) then
            do i=1,m
               tmp=a(i,indr(k))
               a(i,indr(k))=a(i,indc(k))
               a(i,indc(k))=tmp
            enddo
         endif
      enddo
c
      return
      end
c-----------------------------------------------------------------------
      subroutine shiftcg
c
c     prompt user for new cg and translate geometry to this location
c
#     include "basics.inc"
      include 'basicsp.inc'
c
      real cgx,cgy,cgz
      real ngx,ngy,ngz
      character*80 fnm
c
      ntot= nx*ny*nz*nel
      mtot= nelm*8
      if (if3d) then
         call vol_int(cgx,vol,xp)
         call vol_int(cgy,vol,yp)
         call vol_int(cgz,vol,zp)
         cgx = cgx/vol
         cgy = cgy/vol
         cgz = cgz/vol
         call prsrrr('Old CG: $',cgx,cgy,cgz)
c
c        call prs   ('Input new CG: $')
c        call rerrr (ngx,ngy,ngz)
c
         call prs ('Input file containing new matrix:$')
         call blank(fnm,80)
         call res (fnm,80)
         open (unit=88,file=fnm,status='old',err=999)
         read (88,*,err=9999,end=9999) ngx,ngy,ngz
         close (unit=88)
c
         dx = ngx-cgx
         dy = ngy-cgy
         dz = ngz-cgz
         call cadd(xp,dx,ntot)
         call cadd(yp,dy,ntot)
         call cadd(zp,dz,ntot)
         call cadd(x ,dx,mtot)
         call cadd(y ,dy,mtot)
         call cadd(z ,dz,mtot)
      else
         call vol_int(cgx,vol,xp)
         call vol_int(cgy,vol,yp)
         cgx = cgx/vol
         cgy = cgy/vol
         call prsrr ('Old CG: $',cgx,cgy)
         call prs   ('Input new CG: $')
         call rerr  (ngx,ngy)
         dx = ngx-cgx
         dy = ngy-cgy
         call cadd(xp,dx,ntot)
         call cadd(yp,dy,ntot)
         call cadd(x ,dx,mtot)
         call cadd(y ,dy,mtot)
      endif
c
      call gencen
      call coef
      return
c
  999 continue
      call prs('Problem opening file.$')
      return
c
 9999 continue
      call prs('Problem reading file.$')
      close (unit=88)
      return
c
      end
c-----------------------------------------------------------------------
      subroutine shiftmom
c
c     prompt user for new cg and translate geometry to this location
c
#     include "basics.inc"
      include 'basicsp.inc'
c
      real ixx,iyy,izz,ixy,ixz,iyz
      real m(3,3),lam(3),wrk(9)
      real n(3,3),r(3,3)
      integer indr(3),indc(3),ipiv(3)
      character*80 fnm
c
      ntot= nx*ny*nz*nel
      if (if3d) then
         call vol_int(cgx,vol,xp)
         call vol_int(cgy,vol,yp)
         call vol_int(cgz,vol,zp)
         cgx = cgx/vol
         cgy = cgy/vol
         cgz = cgz/vol
         cgxm=-cgx
         cgym=-cgy
         cgzm=-cgz
c
         call cadd2(wkv1,xp,cgxm,ntot)
         call cadd2(wkv2,yp,cgym,ntot)
         call cadd2(wkv3,zp,cgzm,ntot)
         call vol_int2(ixx,vol,wkv1,wkv1)
         call vol_int2(iyy,vol,wkv2,wkv2)
         call vol_int2(izz,vol,wkv3,wkv3)
         call vol_int2(ixy,vol,wkv1,wkv2)
         call vol_int2(ixz,vol,wkv1,wkv3)
         call vol_int2(iyz,vol,wkv2,wkv3)
c
         ixx = ixx/vol
         iyy = iyy/vol
         izz = izz/vol
         ixy = ixy/vol
         ixz = ixz/vol
         iyz = iyz/vol
c
         a = iyy+izz
         b = izz+ixx
         c = ixx+iyy
c
         f = iyz
         g = ixz
         h = ixy
c
         m(1,1) =  a
         m(1,2) = -h
         m(1,3) = -g
c
         m(2,1) = -h
         m(2,2) =  b
         m(2,3) = -f
c
         m(3,1) = -g
         m(3,2) = -f
         m(3,3) =  c
c

         call ssyev('V','U',3,M,3,lam,wrk,9,info)
c
         if (info.ne.0) then
            call prsi('Returning. ssyev info not 0:$',info)
            return
         endif
c
         call prsrrr('Vec 1:$',m(1,1),m(2,1),m(3,1))
         call prsrrr('Vec 2:$',m(1,2),m(2,2),m(3,2))
         call prsrrr('Vec 3:$',m(1,3),m(2,3),m(3,3))
c
         call outmat(m,3,3,'Mom')
c
c        call prs   ('Input new matrix $')
c        call prs   ('Vec 1:$')
c        call rerrr (n(1,1),n(2,1),n(3,1))
c        call prs   ('Vec 2:$')
c        call rerrr (n(1,2),n(2,2),n(3,2))
c        call prs   ('Vec 3:$')
c        call rerrr (n(1,3),n(2,3),n(3,3))
c
         call prs ('Input file containing new matrix:$')
         call blank(fnm,80)
         call res (fnm,80)
         open (unit=88,file=fnm,status='old',err=999)
         read (88,*,err=9999,end=9999) ((n(i,j),j=1,3),i=1,3)
         close (unit=88)
c
c
c        X" = R X'  ==>    R  =  X" (X')^-1
c
         do k=1,3
            call norm3d(m(1,k))
            call norm3d(n(1,k))
         enddo
c
         call outmat   (m,3,3,'mb4')
c        call gaujord  (m,3,3,wrk,info)
         call gaujordf (m,3,3,indr,indc,ipiv,info)
         call outmat   (m,3,3,'maf')
         call outmat   (n,3,3,'nb4')
         call mxm      (n,3,m,3,r,3)
         call outmat   (r,3,3,'Rot')
c
         if (info.ne.0) then
            call prsi('Returning. gaujord info not 0:$',info)
            return
         endif
c
c        Rotation matrix set, apply to every point on the planet
c
         do i=1,ntot
            xn = r(1,1)*wkv1(i)+r(1,2)*wkv2(i)+r(1,3)*wkv3(i)+cgx
            yn = r(2,1)*wkv1(i)+r(2,2)*wkv2(i)+r(2,3)*wkv3(i)+cgy
            zn = r(3,1)*wkv1(i)+r(3,2)*wkv2(i)+r(3,3)*wkv3(i)+cgz
            xp(i) = xn
            yp(i) = yn
            zp(i) = zn
         enddo
c
         do ie=1,nel
         do ic=1,8
            xo = x(ie,ic)-cgx
            yo = y(ie,ic)-cgy
            zo = z(ie,ic)-cgz
            xn = r(1,1)*xo+r(1,2)*yo+r(1,3)*zo
            yn = r(2,1)*xo+r(2,2)*yo+r(2,3)*zo
            zn = r(3,1)*xo+r(3,2)*yo+r(3,3)*zo
            x(ie,ic) = xn + cgx
            y(ie,ic) = yn + cgy
            z(ie,ic) = zn + cgz
         enddo
         enddo
      endif
c
      call gencen
      call coef
      return
c
  999 continue
      call prs('Problem opening file.$')
      return
c
 9999 continue
      call prs('Problem reading file.$')
      close (unit=88)
      return
c
      end
c-----------------------------------------------------------------------
      subroutine mapz(f,regdir)
#     include "basics.inc"
c
      real f(1)
      real gll_pts(nxm),gll_wts(nxm),i(nxm*nxm),it(nxm*nxm)
      real w(nxm*nxm*nxm)
      character*1 regdir
c
      integer e
c
c     This routine maps data from GLL points to uniform points in
c     the z direction only.   This is really a special purpose routine...
c
      call legend(gll_pts,gll_wts,nx)
c
c     Use Fornberg's routine to compute interpolation coefficients
c
      nx1 = nx-1
      l   = 1
      mx  = nx
      do k=1,mx
         z0 = (2.*(k-1))/(mx-1) - 1.
         call fd_weights_full(z0,gll_pts,nx1,0,it(l))
         l  = l+nx
      enddo
      call transpose(i,mx,it,nx)
c
      if (regdir.eq.'z'.or.regdir.eq.'Z') then
         nxy  = nx*ny
         nxyz = nx*ny*nz
         l = 1
         do e=1,nel
            call mxm (f(l),nxy,it,nz,w,nz)
            call copy(f(l),w,nxyz)
            l = l+nxyz
         enddo
      elseif (regdir.eq.'x'.or.regdir.eq.'X') then
         nyz  = ny*nz
         nxyz = nx*ny*nz
         l = 1
         do e=1,nel
            call mxm (i,nx,f(l),nx,w,nyz)
            call copy(f(l),w,nxyz)
            l = l+nxyz
         enddo
      else
         call prs('Sorry, direction is not currently supported$')
         write(6,*) 'Direction: ',regdir
      endif
      return
      end
c-----------------------------------------------------------------------
      subroutine setsrfp_periodic(e,f,e1,f1,ifld)

#     include "basics.inc"

      integer e,f,e1,f1,e0,f0
c
c     Select a pair of planes w / periodic bcs
c
      nlstp=0
      nfaces=2*ndim
c
      nlstp=nlstp+1
      ijkpln=eface1(f)
      ix=1
      if (mod(ijkpln,2).eq.0) ix=nx
      list(nlstp)=6*nx*(e-1)+nx*(ijkpln-1)+ix
      call signpl(list(nlstp))
c
      e1 = bc(1,f,e,ifld)
      f1 = bc(2,f,e,ifld)
c
      e0 = bc(1,f1,e1,ifld)
      f0 = bc(2,f1,e1,ifld)
c
      if (e0.ne.e .or. f0.ne.f) then
         call prsiii('Element RECIPROCITY FAIL:$',e,e0,e1)
         call prsiii('Face    RECIPROCITY FAIL:$',f,f0,f1)
      else
         nlstp=nlstp+1
         ijkpln=eface1(f1)
         ix=1
         if (mod(ijkpln,2).eq.0) ix=nx
         list(nlstp)=6*nx*(e1-1)+nx*(ijkpln-1)+ix
         call signpl(list(nlstp))
         call prsiv('Reciprocity OK:$',e,f,e1,f1)
      endif
c
      return
      end
c-----------------------------------------------------------------------
      subroutine periodic_check

#     include "basics.inc"

      integer e,f,e1,f1,e0,f0
c
c     Select a pair of planes w / periodic bcs
c
      ifld = 2
      if (ifflow) ifld = 1
c
      nfaces = 2*ndim
      do e=1,nel
      do f=1,nfaces
         write(6,*) e,f,cbc(f,e,ifld),' CBC'
         if (cbc(f,e,ifld).eq.'P  ') then
            call setsrfp_periodic(e,f,e1,f1,ifld)
            call tem
            call prsiv('Plotted e,f:$',e,f,e1,f1)
            call prs  ('Continue (y/n) ?$')
            call res  (ans,1)
            if (ans.eq.'n'.or.ans.eq.'N') return
         endif
      enddo
      enddo
c
      return
      end
c-----------------------------------------------------------------------

max_stars_repo_name: RonRahaman/nekRS

max_issues_repo_name: neams-th-coe/nekRS

max_forks_repo_name: neams-th-coe/nekRS


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

C $Log: geom3d.f,v $
C Revision 1.4  2009/03/25 12:36:44  gdsjaar
C Add copyright and license notice to all files.
C Permission to assert copyright has been granted; blot is now open source, BSD
C
C Revision 1.3  2009/01/22 21:34:21  gdsjaar
C There were several inline dbnums common blocks. Replaced with the
C include so they all have the same size with the added variable types.
C
C Added minor support for nodeset and sideset variables.
C
C It can print the count and the names, but that is all
C at this time.
C
C Revision 1.2  2004/10/18 16:30:00  gdsjaar
C Add capability to handle tet elements.
C
C There are still a few minor problems, but the display is substantially correct.
C
C Revision 1.1  1994/04/07 20:01:28  gdsjaar
C Initial checkin of ACCESS/graphics/blotII2
C
c Revision 1.2  1990/12/14  08:51:00  gdsjaar
c Added RCS Id and Log to all files
c
C=======================================================================
      SUBROUTINE GEOM3D (LENF, NLNKF, LINKF, IF2EL,
     &   LENL, LINSET, LINDEF, LREF)
C=======================================================================

C   --*** GEOM3D *** (MESH) Group 3D lines connecting nodes
C   --   Written by Amy Gilkey - revised 03/29/88
C   --              Sam Key, 04/85
C   --
C   --GEOM3D finds all lines making up the surface faces.
C   --The lines are classified and sorted into the following parts:
C   --   -1) Edge
C   --    0) Element block boundary
C   --    n) Interior within element block 'n'
C   --
C   --Parameters:
C   --   LENF - IN - the cumulative face counts by element block
C   --   NLNKF - IN - the number of nodes per face
C   --   LINKF - IN - the connectivity for all faces
C   --   IF2EL - IN - the element number of each face
C   --   LENL - OUT - the cumulative line counts by element block
C   --   LINSET - OUT - the sorted line set, for surface faces only
C   --   LINDEF - SCRATCH - the line set scratch array;
C   --      size = 6 * maximum number of lines
C   --   LREF - SCRATCH - the last line set index for node, length = NUMNPF
C   --
C   --Common Variables:
C   --   Uses NELBLK of /DBNUMS/
C   --   Uses NUMNPF, LLNSET of /D3NUMS/

      common /debugc/ cdebug
      common /debugn/ idebug
      character*8 cdebug

      include 'dbnums.blk'
      COMMON /D3NUMS/ IS3DIM, NNPSUR, NUMNPF, LLNSET
      LOGICAL IS3DIM

      INTEGER LENF(0:NELBLK)
      INTEGER NLNKF(NELBLK)
      INTEGER LINKF(*)
      INTEGER IF2EL(*)
      INTEGER LENL(-2:NELBLK), LINSET(LLNSET,*)
      INTEGER LINDEF(0:5,*)
      INTEGER LREF(*)

      CALL INIINT (NUMNPF, 0, LREF)

C   --Define line set for surface faces only

      NHEAP = 0
      DO 170 IELB = 1, NELBLK
         IXL0 = IDBLNK (IELB, 0, LENF, NLNKF) - 1

C      --Handle normal multi-sided elements

         IF (NLNKF(IELB) .GT. 2) THEN
            DO 130 IFAC = LENF(IELB-1)+1, LENF(IELB)
               IEL = IF2EL(IFAC)

C            --Extract element lines

               N2 = LINKF(IXL0+NLNKF(IELB))
               LREFN2 = LREF(N2)
               LSTREF = LREFN2

               DO 120 ILINK = 1, NLNKF(IELB)
                  N1 = N2
                  N2 = LINKF(IXL0+ILINK)
                  if (n1 .eq. n2) goto 120
                  NMIN = MIN (N1, N2)
                  NMAX = MAX (N1, N2)

C               --Search for line in existing lines

                  LREFN1 = LREFN2
                  IF (ILINK .LT. NLNKF(IELB)) THEN
                     LREFN2 = LREF(N2)
                  ELSE
                     LREFN2 = LSTREF
                  END IF

                  IL = MIN (LREFN1, LREFN2)
  100             CONTINUE
                  IF (IL .GT. 0) THEN
                     IF (LINDEF(1,IL) .EQ. NMIN) THEN
                        IF (LINDEF(2,IL) .EQ. NMAX) GOTO 110
                        IL = LINDEF(4,IL)
                     ELSE IF (LINDEF(1,IL) .EQ. NMAX) THEN
                        IL = LINDEF(4,IL)
                     ELSE
                        IL = LINDEF(5,IL)
                     END IF
                     GOTO 100
                  END IF
  110             CONTINUE

                  IF (IL .GT. 0) THEN
                     IOLELB = IABS (LINDEF(0,IL))

                     IF (IELB .EQ. IOLELB) THEN
C                     --If both faces of the same element (i.e., an edge),
C                     --change from element block n set to edge set,
C                     --else change to element block n set
                        IF (IEL .EQ. LINDEF(3,IL)) THEN
                           LINDEF(0,IL) = -1
                        ELSE
                           LINDEF(0,IL) = IELB
                        END IF

                     ELSE
C                     --If not same element block, change line from
C                     --element block n set to element block boundary set
                        LINDEF(0,IL) = 0

                     END IF
                  ELSE

C                  --No match, add to edge set and point to face

                     NHEAP = NHEAP + 1
                     LINDEF(0,NHEAP) = -IELB
                     LINDEF(1,NHEAP) = NMIN
                     LINDEF(2,NHEAP) = NMAX
                     LINDEF(3,NHEAP) = IEL
                     LINDEF(4,NHEAP) = LREF(NMIN)
                     LINDEF(5,NHEAP) = LREF(NMAX)

                     LREF(N1) = NHEAP
                     LREF(N2) = NHEAP
                  END IF
  120          CONTINUE

               IXL0 = IXL0 + NLNKF(IELB)
  130       CONTINUE

C      --Handle 2-node truss elements

         ELSE IF (NLNKF(IELB) .EQ. 2) THEN
            DO 150 IFAC = LENF(IELB-1)+1, LENF(IELB)
               IEL = IF2EL(IFAC)

C            --Extract element lines for 2-node elements

               N2 = LINKF(IXL0+NLNKF(IELB))
               DO 140 ILINK = 1, NLNKF(IELB)
                  N1 = N2
                  N2 = LINKF(IXL0+ILINK)
                  NMIN = MIN (N1, N2)
                  NMAX = MAX (N1, N2)

C               --Add to edge set as element block n set

                  NHEAP = NHEAP + 1
                  LINDEF(0,NHEAP) = IELB
                  LINDEF(1,NHEAP) = NMIN
                  LINDEF(2,NHEAP) = NMAX
                  LINDEF(3,NHEAP) = IEL
                  LINDEF(4,NHEAP) = 0
                  LINDEF(5,NHEAP) = 0
  140          CONTINUE

               IXL0 = IXL0 + NLNKF(IELB)
  150       CONTINUE

         ELSE
            DO 160 IFAC = LENF(IELB-1)+1, LENF(IELB)
               IXL0 = IXL0 + NLNKF(IELB)
  160       CONTINUE
         END IF
  170 CONTINUE

C   --Count lines in each line set

      DO 180 IELB = -2, NELBLK
         LENL(IELB) = 0
  180 CONTINUE

      DO 190 IL = 1, NHEAP
         IF (LINDEF(0,IL) .LT. -1) LINDEF(0,IL) = -1
         IELB = LINDEF(0,IL)
         LENL(IELB) = LENL(IELB) + 1
  190 CONTINUE

C   --Get starting indices for each line set

      LENL(-2) = 0

      NHEAP = 0
      DO 200 IELB = -1, NELBLK
         L = LENL(IELB)
         LENL(IELB) = NHEAP
         NHEAP = NHEAP + L
  200 CONTINUE

C   --Move lines into appropriate line set

      DO 210 IL = 1, NHEAP
         IELB = LINDEF(0,IL)
         IX = LENL(IELB) + 1
         LENL(IELB) = IX
C...         CALL CPYINT (2, LINDEF(1,IL), LINSET(1,IX))
         linset(1,ix) = lindef(1,il)
         linset(2,ix) = lindef(2,il)
  210 CONTINUE

      RETURN
      END

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


! { dg-do compile }
! { dg-require-effective-target vect_double }

subroutine foo (x,nnd)
  dimension x(nnd)
  integer i

  do i=1,nnd
    x(i) = 1.d0 + (1.d0*i)/nnd
  end do

end subroutine foo

! { dg-final { scan-tree-dump-times "vectorized 1 loops" 1 "vect" { target { vect_unpack  && vect_intfloat_cvt } } } }