text:
! simple code for random number sequences
! from numerical recipes..
function ranfb(ix)
  real :: y,ranfb
  integer :: ix,ix1,ix2,ia1,ia2,ic,j1,j2,j3

  ia1=273
  ia2=25301
  ic=226908345

  if (ix.lt.0) ix=-ix
  if (ix.gt.2**30) ix=MOD(ix,2**30)

  ix1=ix/(2**15)
  ix2=ix-ix1*(2**15)
  j1=MOD(ix1*ia2,2**15)*2**15
  j2=MOD(ix2*ia1,2**15)*2**15
  j3=ix2*ia2
  ix=MOD(j1+j2,2**30)
  ix=MOD(ix+j3,2**30)
  ix=MOD(ix+ic,2**30)
  y=FLOAT(ix)

  ranfb=y/2**30

end function

subroutine setRND
  include 'globals.h'
  integer i
  real ranfb
  do i=1,nrndtable
    rndtable(i)=ranfb(rnseed)
  enddo
end subroutine

function pRND(i)
  include 'globals.h'
  real :: pRND
  integer :: i
  pRND=rndtable(mod(i,nrndtable)+1)
end function

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
!
! The Laboratory of Algorithms
!
! The MIT License
!
! Copyright 2011-2015 Andrey Pudov.
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the 'Software'), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
!

module MFileReader

    use MArrays

    implicit none
    private

    integer, parameter :: DEFAULT_SIZE    = 16
    integer, parameter :: LINE_LENGTH_MAX = 4096

    type, public :: TFileReader
    contains
        procedure, nopass :: readListOfIntegers
        procedure, nopass :: readAdjacencyList
        procedure, nopass :: readAdjacencyWeightedList
    end type
contains
    subroutine readListOfIntegers(name, array)
        character(len=*), intent(in)                       :: name
        integer, dimension(:), allocatable, intent(in out) :: array
        integer, dimension(:), allocatable                 :: temporary_array

        type(TArrays) arrays

        integer unit
        integer status
        integer value
        integer length
        integer index

        allocate(array(DEFAULT_SIZE))
        unit   = 0
        status = 0
        length = 0

        open(unit = unit, file = name)

        do while (status == 0)
            read(unit, '(I)', iostat = status) value

            if (status  == 0) then
                length = length + 1
                call arrays%increaseWhenRequired(array, length)
                array(length) = value
            end if
        end do

        call arrays%resize(array, length)

        close(unit)
    end subroutine

    subroutine readAdjacencyList(name, array)
        character(len=*), intent(in)                         :: name
        integer, dimension(:,:), allocatable, intent(in out) :: array

        type(TArrays) arrays

        character(len=LINE_LENGTH_MAX) line

        integer unit
        integer status
        integer length
        integer index

        integer value
        integer begin
        integer end
        integer first
        integer beginning

        allocate(array(DEFAULT_SIZE, DEFAULT_SIZE))
        unit   = 0
        status = 0
        length = 0

        open(unit = unit, file = name)

        do while (status == 0)
            read(unit, '(A)', iostat = status) line

            if (status  == 0) then
                beginning = .true.

                index = 1
                do while (index <= len(line))
                    begin = index
                    ! accept only ascii digit characters
                    do while ((iachar(line(index:index)) >= 48) .and. (iachar(line(index:index)) <= 57))
                        end = index
                        index = index + 1
                    end do

                    if (begin <= end) then
                        read(line(begin:end), '(I)', iostat = status) value

                        if (beginning == .true.) then
                            beginning = .false.
                            first     = value
                        else
                            call arrays%increaseWhenRequired(array, max(first, value), max(first, value))
                            array(first, value) = 1
                        end if
                    end if

                    index = index + 1
                end do
            end if
        end do

        call arrays%resize(array, first, first)

        close(unit)
    end subroutine

    subroutine readAdjacencyWeightedList(name, array)
        character(len=*), intent(in)                         :: name
        integer, dimension(:,:), allocatable, intent(in out) :: array

        type(TArrays) arrays

        character(len=LINE_LENGTH_MAX) line

        integer unit
        integer status
        integer length
        integer index

        integer vertex
        integer weight

        integer begin
        integer end
        integer first
        integer beginning

        allocate(array(DEFAULT_SIZE, DEFAULT_SIZE))
        unit   = 0
        status = 0
        length = 0

        open(unit = unit, file = name)

        do while (status == 0)
            read(unit, '(A)', iostat = status) line

            if (status  == 0) then
                beginning = .true.
                print *, line

                index = 1
                do while (index <= len(line))
                    begin = index
                    ! accept only ascii digit characters
                    do while ((iachar(line(index:index)) >= 48) .and. (iachar(line(index:index)) <= 57))
                        end = index
                        index = index + 1
                    end do

                    if (begin <= end) then
                        if (beginning == .true.) then
                            read(line(begin:end), '(I)', iostat = status) vertex

                            beginning = .false.
                            first     = vertex
                        else
                            read(line(begin:end), '(I,I)', iostat = status) vertex, weight

                            call arrays%increaseWhenRequired(array, max(first, vertex), max(first, vertex))
                            array(first, vertex) = weight
                        end if
                    end if

                    index = index + 1
                end do
            end if
        end do

        call arrays%resize(array, first, first)

        close(unit)
    end subroutine
end module

ext: f

max_stars_repo_name: andreypudov/Algorithms

max_issues_repo_name: jlokimlin/Algorithms

max_forks_repo_name: jlokimlin/Algorithms


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 p
  stop  42
end program

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
!##############################################################################
!# ****************************************************************************
!# <name> poisson </name>
!# ****************************************************************************
!#
!# <purpose>
!# This program is a simple test program for discretising the equation
!#
!#              - Laplace(u) = f
!#
!# for a scalar function u.
!#
!# There are a couple of examples provided how to solve this problem.
!# Each example has its own characteristics how to solve the problem.
!#
!# There are examples for 1D, 2D and 3D triangulations. The
!# poissonXd_method0_simple modules show the very basic steps how to solve the
!# Poisson problem. On top of that, the poissonXd_method1_XXXX modules
!# extend the poissonXd_method0_simple to more general situations (mmultigrid,
!# fictitious boundary,...), still in the style of poissonXd_method0_simple.
!# The poissonXd_method2_XXXX routines then decompose the different
!# tasks of the problem into different subroutines, thus bringing more
!# structure into the code and highlighting more and more features of the
!# kernel.
!#
!# </purpose>
!##############################################################################

program poisson

  use fsystem
  use genoutput
  use paramlist

  use poisson1d_method0_simple
  use poisson1d_method1_mg
  use poisson2d_method0_agmg
  use poisson2d_method0_block
  use poisson2d_method0_cmsort
  use poisson2d_method0_neumann
  use poisson2d_method0_ravthomas
  use poisson2d_method0_simple
  use poisson2d_method0_smart
  use poisson2d_method1_em30
  use poisson2d_method1_fbc
  use poisson2d_method1_hadapt
  use poisson2d_method1_l2prj
  use poisson2d_method1_mg
  use poisson2d_method1_ncc
  use poisson2d_method1_prolmat
  use poisson2d_method1_robin
  use poisson2d_method2
  use poisson2d_method2_cmsort
  use poisson2d_method2_collect
  use poisson2d_method2_mg
  use poisson3d_method0_agmg
  use poisson3d_method0_simple
  use poisson3d_method1_em30
  use poisson3d_method1_mg

  implicit none

  ! local parameter list
  type(t_parlist) :: rparamlist

  ! local variables
  character(len=SYS_STRLEN) :: slogdir,slogfile,smaster
  integer :: i

  ! The very first thing in every application:
  ! Initialise system-wide settings:

  call sys_init()

  ! Initialise the output system.
  !
  ! Normally, we write all the output to the screen and to a file
  ! "./log/output.txt".
  ! In the case that environment variables "$logdir"/"$resultsfile" exists,
  ! we write all the output to that file. This can be used e.g. in
  ! regression tests to compare results to reference results.
  if (sys_getenv_string("LOGDIR",slogdir) .and. &
      sys_getenv_string("RESULTFILE",slogfile)) then
    call output_init (trim(slogdir)//"/"//trim(slogfile))
  else
    call output_init ("./log/output.txt")
  end if

  ! The very second thing in every program:
  ! Initialise the FEAT 2.0 storage management:
  call storage_init(999, 100)

  ! Initialise the parameter list
  call parlst_init(rparamlist)

  ! Read the master file
  call sys_getcommandLineArg(1,smaster,sdefault='./data/master.dat')
  call parlst_readfromfile(rparamlist, smaster)

  ! Call the problem to solve. Poisson 1D method 1 - simple:
  call parlst_getvalue_int(rparamlist, "", "POISSON1D_METHOD0_SIMPLE", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-1D-Problem with method 0 - simple")
    call output_line ("-----------------------------------------------------")
    call poisson1d_0_simple
  end if

  ! Call the problem to solve. Poisson 1D method 1 - multigrid:
  call parlst_getvalue_int(rparamlist, "", "POISSON1D_METHOD1_MG", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-1D-Problem with method 1 - multigrid")
    call output_line ("--------------------------------------------------------")
    call poisson1d_1_mg
  end if

  ! Call the problem to solve. Poisson 2D method 1 - simple:
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD0_SIMPLE", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 0 - simple")
    call output_line ("-----------------------------------------------------")
    call poisson2d_0_simple
  end if

  ! Call the problem to solve. Poisson 2D method 1 - smart:
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD0_SMART", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 0 - smart")
    call output_line ("----------------------------------------------------")
    call poisson2d_0_smart
  end if

  ! Call the problem to solve. Poisson 2D method 1 - pure Neumann problem:
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD0_NEUMANN", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 0 - Neumann")
    call output_line ("------------------------------------------------------")
    call poisson2d_0_neumann
  end if

  ! Call the problem to solve. Poisson 2D method 1 - CM-sorting:
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD0_CMSORT", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 0 - CM-sorting")
    call output_line ("---------------------------------------------------------")
    call poisson2d_0_cmsort
  end if

  ! Call the problem to solve. Poisson 2D method 0 - block variant:
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD0_BLOCK", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 0 - block")
    call output_line ("-----------------------------------------------------")
    call poisson2d_0_block
  end if

  ! Call the problem to solve. Poisson 2D method 0 - Raviart-Thomas variant:
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD0_RAVTHOMAS", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 0 - RavThomas")
    call output_line ("--------------------------------------------------------")
    call poisson2d_0_ravthomas
  end if

  ! Call the problem to solve. Poisson 2D method 1 - algebraic multigrid:
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD0_AGMG", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 0 - AGMG")
    call output_line ("---------------------------------------------------")
    call poisson2d_0_agmg
  end if

  ! Call the problem to solve. Poisson 2D method 1 - nonconstant coefficients:
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD1_NCC", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 1 - ncc")
    call output_line ("--------------------------------------------------")
    call poisson2d_1_ncc
  end if

  ! Call the problem to solve. Poisson 2D method 1 - multigrid:
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD1_MG", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 1 - multigrid")
    call output_line ("--------------------------------------------------------")
    call poisson2d_1_mg
  end if

  ! Call the problem to solve. Poisson 1: Support for nonconforming elements
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD1_EM30", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 1 - EM30")
    call output_line ("---------------------------------------------------")
    call poisson2d_1_em30
  end if

  ! Call the problem to solve. Poisson 1: Support for nonconforming elements
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD1_ROBIN", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 1 - Robin BC")
    call output_line ("-------------------------------------------------------")
    call poisson2d_1_robin
  end if

  ! Call the problem to solve. Poisson 1: Fictitious boundary support
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD1_FBC", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 1 - FBC")
    call output_line ("--------------------------------------------------")
    call poisson2d_1_fbc
  end if

  ! Call the problem to solve. Poisson 1: h-adaptivity
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD1_HADAPT", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 1 - hadapt")
    call output_line ("-----------------------------------------------------")
    call poisson2d_1_hadapt
  end if

  ! Call the problem to solve. Poisson 2D method 1 - L2-projection:
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD1_L2PRJ", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 1 - L2-projection")
    call output_line ("------------------------------------------------------------")
    call poisson2d_1_l2prj
  end if

  ! Call the problem to solve. Poisson 2D method 1 - Prolongation matrix:
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD1_PROLMAT", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 1 - Prol.-Matrix")
    call output_line ("-----------------------------------------------------------")
    call poisson2d_1_prolmat
  end if

  ! Call the problem to solve. Poisson 2D method 2:
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD2", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 2")
    call output_line ("--------------------------------------------")
    call poisson2d_2
  end if

  ! Call the problem to solve. Poisson 3: Sorting with Cuthill McKee
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD2_CMSORT", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 2 - CM-sorting")
    call output_line ("---------------------------------------------------------")
    call poisson2d_2_cmsort
  end if

  ! Call the problem to solve. Poisson 5:
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD2_MG", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 2 - multigrid")
    call output_line ("--------------------------------------------------------")
    call poisson2d_2_mg
  end if

  ! Call the problem to solve. Poisson 3: Collection support
  call parlst_getvalue_int(rparamlist, "", "POISSON2D_METHOD2_COLLECT", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-2D-Problem with method 2 - collection")
    call output_line ("---------------------------------------------------------")
    call poisson2d_2_collect
  end if

  ! Call the problem to solve. Poisson3D-1:
  call parlst_getvalue_int(rparamlist, "", "POISSON3D_METHOD0_SIMPLE", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-3D-Problem with method 0 - simple")
    call output_line ("-----------------------------------------------------")
    call poisson3d_0_simple
  end if

  ! Call the problem to solve. Poisson3D-1:
  call parlst_getvalue_int(rparamlist, "", "POISSON3D_METHOD0_AGMG", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-3D-Problem with method 0 - AGMG")
    call output_line ("---------------------------------------------------")
    call poisson3d_0_agmg
  end if

  ! Call the problem to solve. Poisson3D-1:
  call parlst_getvalue_int(rparamlist, "", "POISSON3D_METHOD1_MG", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-3D-Problem with method 1 - multigrid")
    call output_line ("--------------------------------------------------------")
    call poisson3d_1_mg
  end if

  ! Call the problem to solve. Poisson3D-7:
  call parlst_getvalue_int(rparamlist, "", "POISSON3D_METHOD1_EM30", i, 0)
  if (i.ne. 0) then
    call output_lbrk ()
    call output_line ("Calculating Poisson-3D-Problem with method 1 - EM30")
    call output_line ("---------------------------------------------------")
    call poisson3d_1_em30
  end if

  ! Release parameter list
  call parlst_done(rparamlist)

  ! Print out heap statistics - just to check if everything
  ! is cleaned up.
  ! This should display "Handles in use=0" and "Memory in use=0"!
  call output_lbrk ()
  call storage_info(.true.)

  ! Clean up the storage management, finish
  call storage_done()

end program

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
module constants

 !Include all modifiable parameters for use below:
use parameters

! Contains all the non-modifiable parameters as well as all
! quantities which never change throughout a simulation
! for the suite of casl f90 codes.

 !Grid dimensions +/- 1 & 2:
integer,parameter:: nxp1=nx+1,nxm1=nx-1,nxm2=nx-2
integer,parameter:: nyp1=ny+1,nym1=ny-1

 !Fine grid used normally in contour -> grid conversion:
integer,parameter:: mgf=4,nxf=mgf*nx,nyf=mgf*ny
 !mgf:  fine grid/coarse grid ratio (4 is required by subroutine
 !      coarsen in contours.f90)
 !Fine grid dimensions +/- 1 & 2:
integer,parameter:: nxfm1=nxf-1,nxfm2=nxf-2
integer,parameter:: nyfp1=nyf+1,nyfm1=nyf-1

 !Ultra-fine grid used in contouring:
integer,parameter:: mgu=16,nxu=mgu*nx,nyu=mgu*ny
 !mgu:  ultra-fine grid/coarse grid ratio (16 is the default)
 !Ultra-fine grid dimensions +/- 1 & 2:
integer,parameter:: nxum1=nxu-1,nxum2=nxu-2
integer,parameter:: nyup1=nyu+1,nyum1=nyu-1

 !For reading & writing direct access data:
integer,parameter:: ngridp=nx*nyp1,nbytes=4*(ngridp+1)

 !Maximum number of contour levels (used in surgery and congen):
integer,parameter:: nlevm=2000
 !nlevm: up to 2*nlevm contour levels are allowed

 !Maximum number of contour nodes:
integer,parameter:: npm=625*nx*ny
 !Maximum number of contours:
integer,parameter:: nm=npm/20+npm/200
 !Maximum number of nodes on any single contour:
integer,parameter:: nprm=npm/10
 !Maximum number of nodes in any contour level:
integer,parameter:: nplm=npm/2

 !Generic double precision numerical constants:
double precision,parameter:: zero=0.0,one=1.0 ,two=2.0,three=3.0
double precision,parameter:: four=4.0,five=5.0,six=6.0,twelve=12.0
double precision,parameter:: f12=one/two,    f13=one/three,f23=two/three
double precision,parameter:: f14=one/four,   f32=three/two,f43=four/three
double precision,parameter:: f53=five/three, f56=five/six, f74=7.0/four
double precision,parameter:: f112=one/twelve,f1112=11.0/twelve,f16=one/six
double precision,parameter:: f18=one/8.0,f116=one/16.0
double precision,parameter:: twopi=two*pi
double precision,parameter:: small=1.e-12,small3=small*small*small
double precision,parameter:: oms=one-small

 !Domain lengths and inverses:
double precision,parameter:: hlx=f12*ellx,hlxi=one/hlx,xmin=-hlx,xmax=hlx
double precision,parameter:: ymax=ymin+elly,ycen=(ymin+ymax)/two
double precision,parameter:: hly=oms*f12*elly,ybeg=ycen-hly

 !Fractional thickness of the upper layer (note: h1 + h2 = 1):
double precision,parameter:: h2=one-h1,h1h2=h1*h2
double precision,parameter:: h1inv=one/h1,h2inv=one/h2

 !Constants used to define relative vorticity and interface displacements:
double precision,parameter:: alphac=one-alpha,kdbarsq=kdbar*kdbar
double precision,parameter:: h1kdbarsq=h1*kdbarsq
double precision,parameter:: h1h2kdbarsq=h1h2*kdbarsq
double precision,parameter:: h1h2ackdbarsq=h1h2*alphac*kdbarsq
double precision,parameter:: massfac=one/(h1+alpha*h2)

 !Define Rossby deformation wavenumbers (for each mode):
double precision,parameter:: gamma1=f12-sqrt(f14-alphac*h1*h2)
double precision,parameter:: gamma2=f12+sqrt(f14-alphac*h1*h2)
double precision,parameter:: kd1sq=gamma1*kdbarsq, kd1=sqrt(kd1sq)
double precision,parameter:: kd2sq=gamma2*kdbarsq, kd2=sqrt(kd2sq)

 !Define layer -> mode transformation coefficients:
double precision,parameter:: vec11=h1/(one-gamma1)          !mode 1, layer 1
double precision,parameter:: vec12=(h2-gamma1)/(one-gamma1) !mode 1, layer 2
double precision,parameter:: vec21=h1/(h2-gamma2)           !mode 2, layer 1
double precision,parameter:: vec22=one                      !mode 2, layer 2
double precision,parameter:: vec2sum=vec21+vec22

 !Define mode -> layer transformation coefficients:
double precision,parameter:: determinant=vec11*vec22-vec12*vec21
double precision,parameter:: vect11=vec22/determinant       !layer 1, mode 1
double precision,parameter:: vect12=-vec12/determinant      !layer 1, mode 2
double precision,parameter:: vect21=-vec21/determinant      !layer 2, mode 1
double precision,parameter:: vect22=vec11/determinant       !layer 2, mode 2

 !Initialise coefficients for decomposing energy & enstrophy into modes:
double precision,parameter:: coeff1=h1*vect11**2+alpha*h2*vect21**2
double precision,parameter:: coeff2=h1*vect12**2+alpha*h2*vect22**2

 !Maximum time step:
double precision,parameter:: dtmax=0.1*min(tgsave,tcsave)

 !Set maximum Rossby wave frequency (used in adapt in evolution.f90):
double precision,parameter:: srwfm=beta/((twopi/ellx)**2+kd1sq)

 !Basic constants:
double precision,parameter:: domarea=ellx*elly,aspect=ellx/elly
double precision,parameter:: glx=ellx/nx,glxi=nx/ellx
double precision,parameter:: gly=elly/ny,glyi=ny/elly
double precision,parameter:: garea=glx*gly,dsumi=one/(nx*ny)

 !Logical control variables:
logical,parameter:: sponge=(thdmax .gt. zero)
logical,parameter:: heating=((rtherm1 .gt. zero) .or. (rtherm2 .gt. zero) .or. sponge)
logical,parameter:: friction=(rekman .gt. zero)
logical,parameter:: damping=(heating .or. friction)
logical,parameter:: barot=(alpha .eq. one)
logical,parameter:: stoch=(eirate .gt. zero)

end module

max_stars_repo_name: julianmak/pydra

max_issues_repo_name: julianmak/pydra

max_forks_repo_name: julianmak/pydra


! glsphere.f90
!
! OpenGL example that rotates the camera around GLU spheres.
!
! Author:  Philipp Engel
! GitHub:  https://github.com/interkosmos/fortran-sdl2/
! Licence: ISC
module camera
    implicit none
    private

    type, public :: camera_class
        real(kind=8) :: x, y, z
    contains
        procedure :: init   => camera_init
        procedure :: update => camera_update
    end type camera_class
contains
    subroutine camera_init(this, x, y, z)
        class(camera_class), intent(inout) :: this
        real(kind=8),        intent(in)    :: x
        real(kind=8),        intent(in)    :: y
        real(kind=8),        intent(in)    :: z

        this%x = x
        this%y = y
        this%z = z
    end subroutine camera_init

    subroutine camera_update(this, angle)
        real(kind=8), parameter            :: PI = acos(-1.0_8)
        class(camera_class), intent(inout) :: this
        real,                intent(in)    :: angle
        real(kind=8)                       :: dir

        dir = angle * PI / 180
        this%x = cos(dir)
        this%z = sin(dir)
    end subroutine camera_update
end module camera

module sphere
    use :: sdl2
    use :: glu
    implicit none
    private

    type, public :: sphere_class
        real         :: x
        real         :: y
        real         :: z
        real(kind=8) :: radius
        integer      :: slices
        integer      :: stacks
        real         :: ambient(4)
        real         :: diffuse(4)
        real         :: specular(4)
        real         :: emission(4)
        real         :: shininess
        type(c_ptr)  :: quadric
    contains
        procedure :: init    => sphere_init
        procedure :: destroy => sphere_destroy
        procedure :: render  => sphere_render
    end type sphere_class
contains
    subroutine sphere_init(this, x, y, z, radius, diffuse)
        class(sphere_class), intent(inout) :: this
        real,                intent(in)    :: x
        real,                intent(in)    :: y
        real,                intent(in)    :: z
        real(kind=8),        intent(in)    :: radius
        real,                intent(in)    :: diffuse(4)

        this%x = x
        this%y = y
        this%z = z

        this%quadric = glunewquadric()
        this%radius  = radius
        this%slices  = 32
        this%stacks  = 32

        this%diffuse   = diffuse
        this%ambient   = [ 0.0, 0.0, 0.0, 1.0 ]
        this%specular  = [ 1.0, 1.0, 1.0, 1.0 ]
        this%emission  = [ 0.0, 0.0, 0.0, 1.0 ]
        this%shininess = 100
    end subroutine sphere_init

    subroutine sphere_destroy(this)
        class(sphere_class), intent(inout) :: this

        call gludeletequadric(this%quadric)
    end subroutine sphere_destroy

    subroutine sphere_render(this)
        class(sphere_class), intent(inout) :: this

        call glpushmatrix()
            call gltranslatef(this%x, this%y, this%z)
            call glmaterialfv(GL_FRONT, GL_AMBIENT,  this%ambient)
            call glmaterialfv(GL_FRONT, GL_DIFFUSE,  this%diffuse)
            call glmaterialfv(GL_FRONT, GL_SPECULAR, this%specular)
            call glmaterialf(GL_FRONT, GL_SHININESS, this%shininess)
            call glmaterialfv(GL_FRONT, GL_EMISSION, this%emission)
            call glusphere(this%quadric, this%radius, this%slices, this%stacks)
        call glpopmatrix()
    end subroutine sphere_render
end module sphere

program main
    use, intrinsic :: iso_c_binding
    use, intrinsic :: iso_fortran_env, only: stderr => error_unit, stdout => output_unit
    use :: sdl2
    use :: glu
    use :: camera
    use :: sphere
    implicit none

    character(len=*), parameter :: WINDOW_TITLE  = 'Fortran SDL 2.0'
    integer,          parameter :: SCREEN_WIDTH  = 1024
    integer,          parameter :: SCREEN_HEIGHT = 768

    character(len=32)        :: title
    type(c_ptr)              :: context
    type(c_ptr)              :: window
    type(sdl_event)          :: event
    type(camera_class)       :: eye
    type(sphere_class)       :: spheres(2)
    integer(kind=1), pointer :: keys(:)
    integer                  :: fc, fps, i, t1

    ! Initialise SDL.
    if (sdl_init(SDL_INIT_EVERYTHING) < 0) then
        write (stderr, '(2a)') 'SDL Error: ', sdl_get_error()
        stop
    end if

    ! Create the SDL window.
    window = sdl_create_window(WINDOW_TITLE // c_null_char, &
                               SDL_WINDOWPOS_UNDEFINED, &
                               SDL_WINDOWPOS_UNDEFINED, &
                               SCREEN_WIDTH, &
                               SCREEN_HEIGHT, &
                               ior(SDL_WINDOW_OPENGL, SDL_WINDOW_SHOWN))

    if (.not. c_associated(window)) then
        write (stderr, '(2a)') 'SDL Error: ', sdl_get_error()
        stop
    end if

    ! Create OpenGL context.
    context = sdl_gl_create_context(window)

    ! Enable VSync.
    if (sdl_gl_set_swap_interval(1) < 0) then
        write (stderr, '(2a)') 'SDL Error: ', sdl_get_error()
    end if

    ! Initialise OpenGL.
    call gl_init(SCREEN_WIDTH, SCREEN_HEIGHT)

    ! Initialise camera.
    call eye%init(1.0_8, 1.0_8, 1.0_8)
    call eye%update(90.0)

    ! Initialise spheres.
    call spheres(1)%init(0.0, 0.0, 0.0, 0.2_8, [ 0.1, 0.5, 0.8, 1.0 ])
    call spheres(2)%init(0.2, 0.2, 0.3, 0.1_8, [ 0.8, 0.0, 0.2, 1.0 ])

    ! Main loop.
    loop: do
        t1 = sdl_get_ticks()

        ! Event handling.
        if (sdl_poll_event(event) > 0) then
            select case (event%type)
                case (SDL_QUITEVENT)
                    exit loop
            end select
        end if

        ! Check keyboard input.
        keys(0:) => sdl_get_keyboard_state()
        if (keys(int(SDL_SCANCODE_ESCAPE, kind=1)) == 1) exit loop

        ! Render the scene.
        call render(eye, spheres)
        call sdl_gl_swap_window(window)

        ! Show frames per second.
        fps = current_fps(t1)
        fc  = fc + 1

        if (fc >= fps) then
            fc = 0
            write (title, '(2a, i2, 2a)') WINDOW_TITLE, ' (', fps, ' FPS)', c_null_char
            call sdl_set_window_title(window, title)
        end if
    end do loop

    do i = 1, size(spheres)
        call spheres(i)%destroy()
    end do

    call sdl_gl_delete_context(context)
    call sdl_destroy_window(window)
    call sdl_quit()
contains
    function current_fps(t1)
        !! Calculates current frames per seconds.
        integer, intent(in) :: t1           !! First time value.
        integer             :: current_fps  !! Result.
        integer             :: dt           !! Time delta.
        integer, save       :: fc           !! Frame counter.
        integer, save       :: ft           !! Frame time.
        integer, save       :: fps          !! Frames per second.

        dt = sdl_get_ticks() - t1

        if (ft >= 1000) then
            fps = fc
            fc  = 0
            ft  = 0
        else
            fc = fc + 1
            ft = ft + dt
        end if

        current_fps = fps
    end function current_fps

    subroutine gl_init(screen_width, screen_height)
        !! Initialises OpenGL.
        integer, intent(in) :: screen_width, screen_height
        real(kind=8)        :: aspect

        ! Set drawing region.
        call glviewport(0, 0, screen_width, screen_height)

        ! Set projection matrix.
        call glmatrixmode(GL_PROJECTION)
        call glloadidentity()
        aspect = real(screen_width, kind=8) / real(screen_height, kind=8)
        call gluperspective(45.0_8, aspect, 0.1_8, 50.0_8)

        ! Set model view matrix.
        call glmatrixmode(GL_MODELVIEW)
        call glenable(GL_DEPTH_TEST)
        call glcullface(GL_BACK)
        call glenable(GL_CULL_FACE)

        ! Set lighting.
        call gllightmodelfv(GL_LIGHT_MODEL_AMBIENT, [ 0.2, 0.2, 0.2 ])
        call gllightmodeli(GL_LIGHT_MODEL_LOCAL_VIEWER, 0)
        call gllightmodeli(GL_LIGHT_MODEL_TWO_SIDE, 1)
        call glenable(GL_LIGHT0)
        call glenable(GL_LIGHTING)

        call glloadidentity()
        call glclearcolor(0.0, 0.0, 0.0, 1.0)
    end subroutine gl_init

    subroutine render(eye, spheres)
        !! The display routine, called every frame.
        type(camera_class), intent(inout) :: eye
        type(sphere_class), intent(inout) :: spheres(:)
        integer                           :: i
        real, save                        :: angle

        angle = modulo(angle + 0.25, 360.0)
        call eye%update(angle)

        call glclear(ior(GL_COLOR_BUFFER_BIT, GL_DEPTH_BUFFER_BIT)) ! Clear the screen and depth buffer.
        call glmatrixmode(GL_MODELVIEW)                             ! Alter model view matrix.
        call glloadidentity()                                       ! Reset current model view matrix.
        call glulookat(eye%x, eye%y, eye%z, &
                       0.0_8, 0.0_8, 0.0_8, &
                       0.0_8, 1.0_8, 0.0_8)

        do i = 1, size(spheres)
            call spheres(i)%render()
        end do
    end subroutine render
end program main


text:
! @@name:	ordered.2f
! @@type:	F-fixed
! @@compilable:	no
! @@linkable:	no
! @@expect:	failure
      SUBROUTINE WORK(I)
      INTEGER I
      END SUBROUTINE WORK

      SUBROUTINE ORDERED_WRONG(N)
      INTEGER N

        INTEGER I
!$OMP   DO ORDERED
        DO I = 1, N
! incorrect because an iteration may not execute more than one ordered region
!$OMP     ORDERED
            CALL WORK(I)
!$OMP     END ORDERED

!$OMP     ORDERED     ! { error "Insert compiler error message text here" }
            CALL WORK(I+1)
!$OMP     END ORDERED
        END DO
      END SUBROUTINE ORDERED_WRONG

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
+incdir+../rtl
+incdir+../tb

../rtl/key_gen_sbox.sv
../rtl/aes_key_expander.sv
../tb/aes_key_expander_tb.sv


! test_pot_gb.f90
! Pair potential, Gay-Berne
MODULE test_pot_module

  !------------------------------------------------------------------------------------------------!
  ! This software was written in 2016/17                                                           !
  ! by Michael P. Allen <m.p.allen@warwick.ac.uk>/<m.p.allen@bristol.ac.uk>                        !
  ! and Dominic J. Tildesley <d.tildesley7@gmail.com> ("the authors"),                             !
  ! to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"),     !
  ! published by Oxford University Press ("the publishers").                                       !
  !                                                                                !
  ! LICENCE                                                                                !
  ! Creative Commons CC0 Public Domain Dedication.                                                 !
  ! To the extent possible under law, the authors have dedicated all copyright and related         !
  ! and neighboring rights to this software to the PUBLIC domain worldwide.                        !
  ! This software is distributed without any warranty.                                             !
  ! You should have received a copy of the CC0 Public Domain Dedication along with this software.  !
  ! If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.                               !
  !                                                                                !
  ! DISCLAIMER                                                                                !
  ! The authors and publishers make no warranties about the software, and disclaim liability       !
  ! for all uses of the software, to the fullest extent permitted by applicable law.               !
  ! The authors and publishers do not recommend use of this software for any purpose.              !
  ! It is made freely available, solely to clarify points made in the text. When using or citing   !
  ! the software, you should not imply endorsement by the authors or publishers.                   !
  !------------------------------------------------------------------------------------------------!

  USE, INTRINSIC :: iso_fortran_env, ONLY : error_unit

  IMPLICIT NONE
  PRIVATE

  ! Public routine
  PUBLIC :: force

  ! Public data
  INTEGER, PARAMETER, PUBLIC :: n = 2 ! Pair potential

CONTAINS

  SUBROUTINE force  ( r, e, pot, f, t )
    USE maths_module, ONLY : cross_product
    IMPLICIT NONE
    REAL, DIMENSION(:,:),           INTENT(in)  :: r, e
    REAL,                           INTENT(out) :: pot
    REAL, DIMENSION(:,:), OPTIONAL, INTENT(out) :: f, t

    ! Parameters of the Gay-Berne potential
    !
    ! The key parameters are
    !    mu, nu  ................ the exponents
    !    kappa and kappa' ....... the anisotropies
    !    kappa is the ratio of intermolecular separations
    !    sigma_e / sigma_s  i.e. end-to-end / side-by-side
    !    kappa' is the ratio of well depths
    !    epsilon_s / epsilon_e  i.e. side-by-side / end-to-end
    ! The derived parameters are chi and chi'
    !    chi = (kappa**2 - 1) / (kappa**2+1)
    !    chi' = (z - 1) / (z + 1)
    !    where z = (kappa') ** ( 1 / mu )
    !
    ! For convenience kappa' is spelt xappa, chi' is spelt xhi
    ! We choose units such that sigma_s = 1.0 and epsilon_0 = 1.0
    ! Two of the following three varieties should be commented out

    ! Original Gay-Berne-deMiguel potential [J. Chem. Phys, 74, 3316; Mol. Phys. 74, 405 (1991)]
    INTEGER, PARAMETER :: mu = 2, nu = 1
    REAL,    PARAMETER :: kappa = 3.0, xappa = 5.0

!!$    ! Luckhurst-Phippen potential [Liq. Cryst., 8, 451 (1990)]
!!$    INTEGER, PARAMETER :: mu = 1, nu = 2
!!$    REAL,    PARAMETER :: kappa = 3.0, xappa = 5.0

!!$    ! Berardi-Zannoni potential [J. Chem. Soc. Faraday Trans., 89, 4069 (1993)]
!!$    INTEGER, PARAMETER :: mu = 1, nu = 3
!!$    REAL,    PARAMETER :: kappa = 3.0, xappa = 5.0

    REAL, PARAMETER :: chi = (kappa**2 - 1.0) / (kappa**2+1.0)
    REAL, PARAMETER :: xhi = (xappa**(1.0/mu) - 1.0) / (xappa**(1.0/mu) + 1.0)

    REAL, DIMENSION(3) :: rij, sij, fij, gi, gj
    REAL               :: rij_mag, rij_sq, ci, cj, cij, cp, cm, prefac
    REAL               :: cpchi, cmchi, sigma
    REAL               :: cpxhi, cmxhi, eps1, eps2, epsilon
    REAL               :: dsig_dci, dsig_dcj, dsig_dcij
    REAL               :: deps_dci, deps_dcj, deps_dcij
    REAL               :: dpot_dci, dpot_dcj, dpot_dcij, dpot_drij
    REAL               :: rho, rho6, rho12, rhoterm, drhoterm, cutterm, dcutterm
    REAL,    PARAMETER :: r_cut = 4.0   ! Normally would use a larger value of r_cut
    REAL,    PARAMETER :: tol = 1.e-6
    INTEGER, PARAMETER :: i = 1, j = 2 ! Notation to match appendix

    ! Routine to demonstrate the calculation of forces and torques from the
    ! Gay-Berne potential, including the spherical cutoff contribution.
    ! Written for ease of comparison with the text, rather than efficiency!

    ! Check dimensions to be sure
    IF ( ANY ( SHAPE(r) /= [3,n] ) ) THEN
       WRITE ( unit=error_unit, fmt='(a,4i15)' ) 'r shape error', SHAPE(r), 3, n
       STOP 'Error in test_pot_gb'
    END IF
    IF ( ANY ( SHAPE(e) /= [3,n] ) ) THEN
       WRITE ( unit=error_unit, fmt='(a,4i15)' ) 'e shape error', SHAPE(e), 3, n
       STOP 'Error in test_pot_gb'
    END IF

    ! Check unit vectors
    IF ( ABS(SUM(e(:,i)**2)-1.0) > tol .OR. ABS(SUM(e(:,j)**2)-1.0) > tol ) THEN
       WRITE ( unit=error_unit, fmt='(a)'      ) 'Warning, non-unit vectors'
       WRITE ( unit=error_unit, fmt='(4f10.5)' ) e(:,i), SUM(e(:,i)**2)
       WRITE ( unit=error_unit, fmt='(4f10.5)' ) e(:,j), SUM(e(:,j)**2)
    END IF

    rij     = r(:,i) - r(:,j)
    rij_sq  = SUM ( rij**2 )
    rij_mag = SQRT(rij_sq)  ! Magnitude of separation
    sij     = rij / rij_mag ! Unit vector along rij

    ! Orientation-dependent terms
    ci  = DOT_PRODUCT ( sij, e(:,i) )
    cj  = DOT_PRODUCT ( sij, e(:,j) )
    cij = DOT_PRODUCT ( e(:,i), e(:,j) )
    cp  = ci + cj
    cm  = ci - cj

    ! Sigma formula
    cpchi = cp/(1.0+chi*cij)
    cmchi = cm/(1.0-chi*cij)
    sigma = 1.0/SQRT(1.0-0.5*chi*(cp*cpchi+cm*cmchi))

    ! Epsilon formula
    eps1    = 1.0/SQRT(1.0-(chi*cij)**2) ! Depends on chi, not xhi
    cpxhi   = cp/(1.0+xhi*cij)
    cmxhi   = cm/(1.0-xhi*cij)
    eps2    = 1.0-0.5*xhi*(cp*cpxhi+cm*cmxhi) ! Depends on xhi
    epsilon = (eps1**nu) * (eps2**mu)

    ! Potential at rij
    rho      = rij_mag - sigma + 1.0
    rho6     = 1.0 / rho**6
    rho12    = rho6**2
    rhoterm  = 4.0*(rho12 - rho6)                 ! Needed for forces and torques
    drhoterm = -24.0 * (2.0 * rho12 - rho6) / rho ! Needed for forces and torques
    pot      = epsilon*rhoterm

    ! Potential at r_cut
    rho      = r_cut - sigma + 1.0
    rho6     = 1.0 / rho**6
    rho12    = rho6**2
    cutterm  = 4.0*(rho12 - rho6)                 ! Needed for cutoff forces and torques
    dcutterm = -24.0 * (2.0 * rho12 - rho6) / rho ! Needed for cutoff forces and torques
    pot      = pot - epsilon * cutterm

    IF ( .NOT. PRESENT(f) ) RETURN

    IF ( .NOT. PRESENT(t) ) THEN
       WRITE ( unit=error_unit, fmt='(a)' ) 'Both f and t expected'
       STOP 'Error in test_pot_gb'
    END IF
    IF ( ANY ( SHAPE(f) /= [3,n] ) ) THEN
       WRITE ( unit=error_unit, fmt='(a,4i15)' ) 'f shape error', SHAPE(f), 3, n
       STOP 'Error in test_pot_gb'
    END IF
    IF ( ANY ( SHAPE(t) /= [3,n] ) ) THEN
       WRITE ( unit=error_unit, fmt='(a,4i15)' ) 't shape error', SHAPE(t), 3, n
       STOP 'Error in test_pot_gb'
    END IF

    ! Derivatives of sigma
    prefac    = 0.5*chi*sigma**3
    dsig_dci  = prefac*(cpchi+cmchi)
    dsig_dcj  = prefac*(cpchi-cmchi)
    prefac    = prefac*(0.5*chi)
    dsig_dcij = -prefac*(cpchi**2-cmchi**2)

    ! Derivatives of epsilon
    prefac    = -mu*xhi*(eps1**nu)*eps2**(mu-1)
    deps_dci  = prefac*(cpxhi+cmxhi)
    deps_dcj  = prefac*(cpxhi-cmxhi)
    prefac    = prefac*(0.5*xhi)
    deps_dcij = -prefac*(cpxhi**2-cmxhi**2)                           ! From derivative of eps2
    deps_dcij = deps_dcij + nu*(chi**2)*(eps1**(nu+2))*(eps2**mu)*cij ! From derivative of eps1

    ! Derivatives of potential
    dpot_drij  = epsilon * drhoterm
    dpot_dci   = rhoterm * deps_dci  - epsilon * drhoterm * dsig_dci
    dpot_dcj   = rhoterm * deps_dcj  - epsilon * drhoterm * dsig_dcj
    dpot_dcij  = rhoterm * deps_dcij - epsilon * drhoterm * dsig_dcij

    ! Standard formula for forces and torque gradients
    fij = -dpot_drij*sij - dpot_dci*(e(:,i)-ci*sij)/rij_mag - dpot_dcj*(e(:,j)-cj*sij)/rij_mag
    gi  = dpot_dci*sij  + dpot_dcij*e(:,j)
    gj  = dpot_dcj*sij  + dpot_dcij*e(:,i)

    ! Derivatives of potential at cutoff
    dpot_drij  = epsilon * dcutterm
    dpot_dci   = cutterm * deps_dci  - epsilon * dcutterm * dsig_dci
    dpot_dcj   = cutterm * deps_dcj  - epsilon * dcutterm * dsig_dcj
    dpot_dcij  = cutterm * deps_dcij - epsilon * dcutterm * dsig_dcij

    ! Standard formula for forces and torque gradients (without dpot_drij term)
    fij(:) = fij(:) + dpot_dci*(e(:,i)-ci*sij)/rij_mag + dpot_dcj*(e(:,j)-cj*sij)/rij_mag
    gi(:)  = gi(:) - ( dpot_dci*sij  + dpot_dcij*e(:,j) )
    gj(:)  = gj(:) - ( dpot_dcj*sij  + dpot_dcij*e(:,i) )

    ! Forces
    f(:,i) = fij
    f(:,j) = -fij

    ! Torques
    t(:,i) = -cross_product ( e(:,i), gi ) ! Torque on i due to j
    t(:,j) = -cross_product ( e(:,j), gj ) ! Torque on j due to i

  END SUBROUTINE force

END MODULE test_pot_module


text: Describe Users/LisaStevens here.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*
*     ------------------------------------------------------------------
*	G T R A C A
*     ------------------------------------------------------------------
*
      SUBROUTINE GTRACA(I,L,IFIRST,NFOUND,BUF,
     &           IBUF,LBUF,JLEI,LU)
*
* 1s inactive, 2 electrons in (1s 1p), Singlet S
*
* The following is the list in the order : l  j i   r l coef
*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER ( NCOUP =2)
      DIMENSION IRACAH(5,NCOUP)
      DIMENSION  RACAH(NCOUP)
      DATA  IRACAH/
     &0, 2,2, 1,1,
     &1, 1,1, 2,2/
*234567
*. Output
      DIMENSION BUF(LBUF),IBUF(4,LBUF)
      PARAMETER (SQRT2 = 1.414D0     )
*
      DATA RACAH/
     & 2.0D0,2.0D0/
*
*. Position so IFIRST -1 elements have been skipped
*

      IF(IFIRST.NE.1) THEN
       ISKIP = 0
       DO 100 IELMNT = 1, NCOUP
        JELMNT = IELMNT
        IF(IRACAH(3,IELMNT) .EQ. I .AND.
     &     IRACAH(1,IELMNT).EQ.L ) THEN
*
          IF(JLEI.EQ.0) THEN
            ISKIP = ISKIP +1
          ELSE IF (JLEI.EQ.1.AND.
     &    IRACAH(2,IELMNT).LT.IRACAH(3,IELMNT))THEN
            ISKIP = ISKIP + 1
          ELSE IF (JLEI.EQ.2.AND.
     &    IRACAH(2,IELMNT).EQ.IRACAH(3,IELMNT))THEN
            ISKIP = ISKIP + 1
          ELSE IF (JLEI.EQ.3.AND.
     &    IRACAH(2,IELMNT).LE.IRACAH(3,IELMNT))THEN
            ISKIP = ISKIP + 1
          END IF
          IF(ISKIP.EQ.IFIRST-1) GOTO 101
*
        END IF
  100  CONTINUE
  101  CONTINUE
      ELSE IF ( IFIRST .EQ. 1 ) THEN
         JELMNT = 0
      END IF
*
* Obtain the next elements, atmost LBUF
*
       NFOUND = 0
       DO 200 IELMNT = JELMNT+1, NCOUP
        IF(IRACAH(3,IELMNT) .EQ. I .AND.
     &     IRACAH(1,IELMNT).EQ.L ) THEN
          IF(JLEI.EQ.0) THEN
           NFOUND = NFOUND + 1
           IBUF(1,NFOUND) = IRACAH(2,IELMNT)
           IBUF(2,NFOUND) = IRACAH(3,IELMNT)
           IBUF(3,NFOUND) = IRACAH(4,IELMNT)
           IBUF(4,NFOUND) = IRACAH(5,IELMNT)
            BUF(  NFOUND) =  RACAH(  IELMNT)
          ELSE IF (JLEI.EQ.1.AND.
     &    IRACAH(2,IELMNT).LT.IRACAH(3,IELMNT))THEN
           NFOUND = NFOUND + 1
           IBUF(1,NFOUND) = IRACAH(2,IELMNT)
           IBUF(2,NFOUND) = IRACAH(3,IELMNT)
           IBUF(3,NFOUND) = IRACAH(4,IELMNT)
           IBUF(4,NFOUND) = IRACAH(5,IELMNT)
            BUF(  NFOUND) =  RACAH(  IELMNT)
          ELSE IF (JLEI.EQ.2.AND.
     &    IRACAH(2,IELMNT).EQ.IRACAH(3,IELMNT))THEN
           NFOUND = NFOUND + 1
           IBUF(1,NFOUND) = IRACAH(2,IELMNT)
           IBUF(2,NFOUND) = IRACAH(3,IELMNT)
           IBUF(3,NFOUND) = IRACAH(4,IELMNT)
           IBUF(4,NFOUND) = IRACAH(5,IELMNT)
            BUF(  NFOUND) =  RACAH(  IELMNT)
          ELSE IF (JLEI.EQ.3.AND.
     &    IRACAH(2,IELMNT).LE.IRACAH(3,IELMNT))THEN
           NFOUND = NFOUND + 1
           IBUF(1,NFOUND) = IRACAH(2,IELMNT)
           IBUF(2,NFOUND) = IRACAH(3,IELMNT)
           IBUF(3,NFOUND) = IRACAH(4,IELMNT)
           IBUF(4,NFOUND) = IRACAH(5,IELMNT)
            BUF(  NFOUND) =  RACAH(  IELMNT)
          END IF
          IF(NFOUND.EQ.LBUF   ) GOTO 201
        END IF
  200  CONTINUE
  201  CONTINUE
*
       NTEST = 0
       IF( NTEST.NE.0) THEN
        WRITE(6,*) ' GTRAC1 to your service '
        WRITE(6,*) ' ======================='
        WRITE(6,*) 'Number of elements obtained',NFOUND
        WRITE(6,*) ' IBUF and BUF '
        DO 300 IELMNT = 1, NFOUND
          WRITE(6,'(E12.7,3X,4I4)')
     &    BUF(IELMNT),(IBUF(K,IELMNT),K=1,4)
  300  CONTINUE
*
      END IF
      RETURN
      END


text:
program direct_formaated_write
    use, intrinsic :: iso_fortran_env, only : error_unit
    use :: format_mod
    implicit none
    integer(kind=I4), parameter :: nr_recs = 11
    integer :: unit_nr
    integer(kind=I4) :: idx, istat
    real(kind=DP) :: x
    character(len=1024) :: msg, file_name

    if (command_argument_count() /= 1) then
        write (unit=error_unit, fmt='(A)') 'error: file name expected'
        stop 1
    end if
    call get_command_argument(1, file_name)
    open (newunit=unit_nr, file=trim(file_name), form='formatted', &
          access='direct', recl=rec_width, action='write', status='new', &
          iostat=istat, iomsg=msg)
    if (istat /= 0) then
        write (unit=error_unit, fmt='(2A)') 'error: ', trim(msg)
        stop 2
    end if
    x = 0.5_DP
    do idx = nr_recs, 1, -1
        write (unit=unit_nr, fmt=fmt_str, rec=idx) idx, x
        x = x + 1.3_DP
    end do
    close (unit=unit_nr)

end program direct_formaated_write

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC


text:
subroutine handle_special_cases (ncidt, ncido, nvars, ntime, unlimdimid)
  use shr_kind_mod, only: r8 => shr_kind_r8

  implicit none

  include 'netcdf.inc'
!
! Input arguments
!
  integer ncidt, ncido
  integer nvars, ntime
  integer unlimdimid
!
! Local workspace
!
  character*(nf_max_name) :: name

  integer i, n    ! indices
  integer vt, vo
  integer natts
  integer xtype
  integer nvdims
  integer dimlen
  integer totsiz
  integer vardids(nf_max_var_dims) ! variable dimension ids
  integer start(nf_max_var_dims)
  integer count(nf_max_var_dims)
  integer ret

  character*(nf_max_name) :: attname
!
! Externals
!
  logical, external :: is_special_case
!
! Copy all the special case variables from template to output
! Probably should have separate "define" and "put" loops to optimize redef/enddef
!
  start(:) = 1

  do vt=1,nvars
    call wrap_inq_var (ncidt, vt, name, xtype, nvdims, vardids, natts)

    if (is_special_case (name, ncidt)) then
      totsiz = 1                       ! Init to size that can be multiplied
      do n=1,nvdims
        if (vardids(n) == unlimdimid) then
          count(n) = ntime
          totsiz = totsiz * ntime
        else
          call wrap_inq_dimlen (ncidt, vardids(n), dimlen)
          count(n) = dimlen
          totsiz = totsiz * dimlen
        end if
      end do

      call wrap_def_var (ncido, name, xtype, nvdims, vardids, vo)
!
! Copy attributes from input to output, then the variable itself
!
      do n=1,natts
        call wrap_inq_attname (ncidt, vt, n, attname)
        call wrap_copy_att (ncidt, vt, attname, ncido, vo)
      end do

      if (nf_enddef (ncido) /= NF_NOERR) stop 999

      call cpvar (ncidt, ncido, vt, vo, name, &
                  totsiz, xtype, start, count )

      ret = nf_redef (ncido)

    end if
  end do
end subroutine handle_special_cases

max_stars_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_issues_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_forks_repo_licenses: zlib-acknowledgement FTL RSA-MD


text:
** Copyright (c) 1989, 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.

*--- Large loop with many induction pointers. Can cause block breaking
*--- and register usage problems.
	program ka42
	common /res/iexp(80), ires(80)
	common /com/ a(10), b(10), c(10), d(10), e(10)
	common /com/ f(10), g(10), h(10), p(10), q(10)
	common /com/ r(10), s(10), t(10), u(10), v(10)
	common /com/ w(10), x(10), y(10), z(10), za(10)
	common /com/ ab(10), bb(10), cb(10), db(10), eb(10)
	common /com/ fb(10), gb(10), hb(10), pb(10), qb(10)
	common /com/ rb(10), sb(10), tb(10), ub(10), vb(10)
	common /com/ wb(10), xb(10), yb(10), zb(10), zab(10)
	common /com/ ac(10), bc(10), cc(10), dc(10), ec(10)
	common /com/ fc(10), gc(10), hc(10), pc(10), qc(10)
	common /com/ rc(10), sc(10), tc(10), uc(10), vc(10)
	common /com/ wc(10), xc(10), yc(10), zc(10), zac(10)
	common /com/ ad(10), bd(10), cd(10), dd(10), ed(10)
	common /com/ fd(10), gd(10), hd(10), pd(10), qd(10)
	common /com/ rd(10), sd(10), td(10), ud(10), vd(10)
	common /com/ wd(10), xd(10), yd(10), zd(10), zad(10)

	do 10 i = 2,10
	    a(i) = a(i-1) + b(i)
	    c(i) = d(i) + a(i) * e(i)
	    f(i) = f(i-1) + g(i)
	    h(i) = p(i) + f(i) * q(i)
	    r(i) = r(i-1) + s(i)
	    t(i) = u(i) + r(i) * v(i)
	    w(i) = w(i-1) + x(i)
	    y(i) = z(i) + w(i) * za(i)
	    ab(i) = ab(i-1) + bb(i)
	    cb(i) = db(i) + ab(i) * eb(i)
	    fb(i) = fb(i-1) + gb(i)
	    hb(i) = pb(i) + fb(i) * qb(i)
	    rb(i) = rb(i-1) + sb(i)
	    tb(i) = ub(i) + rb(i) * vb(i)
	    wb(i) = wb(i-1) + xb(i)
	    yb(i) = zb(i) + wb(i) * zab(i)
	    ac(i) = ac(i-1) + bc(i)
	    cc(i) = dc(i) + ac(i) * ec(i)
	    fc(i) = fc(i-1) + gc(i)
	    hc(i) = pc(i) + fc(i) * qc(i)
	    rc(i) = rc(i-1) + sc(i)
	    tc(i) = uc(i) + rc(i) * vc(i)
	    wc(i) = wc(i-1) + xc(i)
	    yc(i) = zc(i) + wc(i) * zac(i)
	    ad(i) = ad(i-1) + bd(i)
	    cd(i) = dd(i) + ad(i) * ed(i)
	    fd(i) = fd(i-1) + gd(i)
	    hd(i) = pd(i) + fd(i) * qd(i)
	    rd(i) = rd(i-1) + sd(i)
	    td(i) = ud(i) + rd(i) * vd(i)
	    wd(i) = wd(i-1) + xd(i)
	    yd(i) = zd(i) + wd(i) * zad(i)
10	continue

	call dosums

	call check (ires, iexp, 80)
	end
	subroutine dosums
	common /com/ a(10), b(10), c(10), d(10), e(10)
	common /com/ f(10), g(10), h(10), p(10), q(10)
	common /com/ r(10), s(10), t(10), u(10), v(10)
	common /com/ w(10), x(10), y(10), z(10), za(10)
	common /com/ ab(10), bb(10), cb(10), db(10), eb(10)
	common /com/ fb(10), gb(10), hb(10), pb(10), qb(10)
	common /com/ rb(10), sb(10), tb(10), ub(10), vb(10)
	common /com/ wb(10), xb(10), yb(10), zb(10), zab(10)
	common /com/ ac(10), bc(10), cc(10), dc(10), ec(10)
	common /com/ fc(10), gc(10), hc(10), pc(10), qc(10)
	common /com/ rc(10), sc(10), tc(10), uc(10), vc(10)
	common /com/ wc(10), xc(10), yc(10), zc(10), zac(10)
	common /com/ ad(10), bd(10), cd(10), dd(10), ed(10)
	common /com/ fd(10), gd(10), hd(10), pd(10), qd(10)
	common /com/ rd(10), sd(10), td(10), ud(10), vd(10)
	common /com/ wd(10), xd(10), yd(10), zd(10), zad(10)
	call sum (a)		! test 1
	call sum (b)		! test 2
	call sum (c)		! test 3
	call sum (d)		! test 4
	call sum (e)		! test 5
	call sum (f)		! test 6
	call sum (g)		! test 7
	call sum (h)		! test 8
	call sum (p)		! test 9
	call sum (q)		! test 10
	call sum (r)		! test 11
	call sum (s)		! test 12
	call sum (t)		! test 13
	call sum (u)		! test 14
	call sum (v)		! test 15
	call sum (w)		! test 16
	call sum (x)		! test 17
	call sum (y)		! test 18
	call sum (z)		! test 19
	call sum (za)		! test 20
	call sum (ab)		! test 21
	call sum (bb)		! test 22
	call sum (cb)		! test 23
	call sum (db)		! test 24
	call sum (eb)		! test 25
	call sum (fb)		! test 26
	call sum (gb)		! test 27
	call sum (hb)		! test 28
	call sum (pb)		! test 29
	call sum (qb)		! test 30
	call sum (rb)		! test 31
	call sum (sb)		! test 32
	call sum (tb)		! test 33
	call sum (ub)		! test 34
	call sum (vb)		! test 35
	call sum (wb)		! test 36
	call sum (xb)		! test 37
	call sum (yb)		! test 38
	call sum (zb)		! test 39
	call sum (ac)		! test 40
	call sum (bc)		! test 41
	call sum (cc)		! test 42
	call sum (dc)		! test 43
	call sum (ec)		! test 44
	call sum (fc)		! test 45
	call sum (gc)		! test 46
	call sum (hc)		! test 47
	call sum (pc)		! test 48
	call sum (qc)		! test 49
	call sum (rc)		! test 50
	call sum (sc)		! test 51
	call sum (tc)		! test 52
	call sum (uc)		! test 53
	call sum (vc)		! test 54
	call sum (wc)		! test 55
	call sum (xc)		! test 56
	call sum (yc)		! test 57
	call sum (zc)		! test 58
	call sum (ad)		! test 59
	call sum (bd)		! test 60
	call sum (cd)		! test 61
	call sum (dd)		! test 62
	call sum (ed)		! test 63
	call sum (fd)		! test 64
	call sum (gd)		! test 65
	call sum (hd)		! test 66
	call sum (pd)		! test 67
	call sum (qd)		! test 68
	call sum (rd)		! test 69
	call sum (sd)		! test 70
	call sum (td)		! test 71
	call sum (ud)		! test 72
	call sum (vd)		! test 73
	call sum (wd)		! test 74
	call sum (xd)		! test 75
	call sum (yd)		! test 76
	call sum (zd)		! test 77
	call sum (zab)		! test 78
	call sum (zac)		! test 79
	call sum (zad)		! test 80
	end
	subroutine sum (a)
	real a(10)
	common /res/iexp(80), ires(80)
	integer icnt
	save icnt
	data icnt/0/
	fsum = 0.0
	do 10 i = 1, 10
10	fsum = fsum + a(i)
	icnt = icnt + 1
	ires(icnt) = fsum
	end
	blockdata
	common /res/iexp(80), ires(80)
	common /com/ a(10), b(10), c(10), d(10), e(10)
	common /com/ f(10), g(10), h(10), p(10), q(10)
	common /com/ r(10), s(10), t(10), u(10), v(10)
	common /com/ w(10), x(10), y(10), z(10), za(10)
	common /com/ ab(10), bb(10), cb(10), db(10), eb(10)
	common /com/ fb(10), gb(10), hb(10), pb(10), qb(10)
	common /com/ rb(10), sb(10), tb(10), ub(10), vb(10)
	common /com/ wb(10), xb(10), yb(10), zb(10), zab(10)
	common /com/ ac(10), bc(10), cc(10), dc(10), ec(10)
	common /com/ fc(10), gc(10), hc(10), pc(10), qc(10)
	common /com/ rc(10), sc(10), tc(10), uc(10), vc(10)
	common /com/ wc(10), xc(10), yc(10), zc(10), zac(10)
	common /com/ ad(10), bd(10), cd(10), dd(10), ed(10)
	common /com/ fd(10), gd(10), hd(10), pd(10), qd(10)
	common /com/ rd(10), sd(10), td(10), ud(10), vd(10)
	common /com/ wd(10), xd(10), yd(10), zd(10), zad(10)
	data a /1,2,3,4,5,6,7,8,9,10/
	data b /1,2,3,4,5,6,7,8,9,10/
	data c /1,2,3,4,5,6,7,8,9,10/
	data d /1,2,3,4,5,6,7,8,9,10/
	data e /1,2,3,4,5,6,7,8,9,10/
	data f /1,2,3,4,5,6,7,8,9,10/
	data g /1,2,3,4,5,6,7,8,9,10/
	data h /1,2,3,4,5,6,7,8,9,10/
	data p /1,2,3,4,5,6,7,8,9,10/
	data q /1,2,3,4,5,6,7,8,9,10/
	data r /1,2,3,4,5,6,7,8,9,10/
	data s /1,2,3,4,5,6,7,8,9,10/
	data t /1,2,3,4,5,6,7,8,9,10/
	data u /1,2,3,4,5,6,7,8,9,10/
	data v /1,2,3,4,5,6,7,8,9,10/
	data w /1,2,3,4,5,6,7,8,9,10/
	data x /1,2,3,4,5,6,7,8,9,10/
	data y /1,2,3,4,5,6,7,8,9,10/
	data z /1,2,3,4,5,6,7,8,9,10/
	data za/1,2,3,4,5,6,7,8,9,10/
	data ab/1,2,3,4,5,6,7,8,9,10/
	data bb/1,2,3,4,5,6,7,8,9,10/
	data cb/1,2,3,4,5,6,7,8,9,10/
	data db/1,2,3,4,5,6,7,8,9,10/
	data eb/1,2,3,4,5,6,7,8,9,10/
	data fb/1,2,3,4,5,6,7,8,9,10/
	data gb/1,2,3,4,5,6,7,8,9,10/
	data hb/1,2,3,4,5,6,7,8,9,10/
	data pb/1,2,3,4,5,6,7,8,9,10/
	data qb/1,2,3,4,5,6,7,8,9,10/
	data rb/1,2,3,4,5,6,7,8,9,10/
	data sb/1,2,3,4,5,6,7,8,9,10/
	data tb/1,2,3,4,5,6,7,8,9,10/
	data ub/1,2,3,4,5,6,7,8,9,10/
	data vb/1,2,3,4,5,6,7,8,9,10/
	data wb/1,2,3,4,5,6,7,8,9,10/
	data xb/1,2,3,4,5,6,7,8,9,10/
	data yb/1,2,3,4,5,6,7,8,9,10/
	data zb/1,2,3,4,5,6,7,8,9,10/
	data ac/1,2,3,4,5,6,7,8,9,10/
	data bc/1,2,3,4,5,6,7,8,9,10/
	data cc/1,2,3,4,5,6,7,8,9,10/
	data dc/1,2,3,4,5,6,7,8,9,10/
	data ec/1,2,3,4,5,6,7,8,9,10/
	data fc/1,2,3,4,5,6,7,8,9,10/
	data gc/1,2,3,4,5,6,7,8,9,10/
	data hc/1,2,3,4,5,6,7,8,9,10/
	data pc/1,2,3,4,5,6,7,8,9,10/
	data qc/1,2,3,4,5,6,7,8,9,10/
	data rc/1,2,3,4,5,6,7,8,9,10/
	data sc/1,2,3,4,5,6,7,8,9,10/
	data tc/1,2,3,4,5,6,7,8,9,10/
	data uc/1,2,3,4,5,6,7,8,9,10/
	data vc/1,2,3,4,5,6,7,8,9,10/
	data wc/1,2,3,4,5,6,7,8,9,10/
	data xc/1,2,3,4,5,6,7,8,9,10/
	data yc/1,2,3,4,5,6,7,8,9,10/
	data zc/1,2,3,4,5,6,7,8,9,10/
	data ad/1,2,3,4,5,6,7,8,9,10/
	data bd/1,2,3,4,5,6,7,8,9,10/
	data cd/1,2,3,4,5,6,7,8,9,10/
	data dd/1,2,3,4,5,6,7,8,9,10/
	data ed/1,2,3,4,5,6,7,8,9,10/
	data fd/1,2,3,4,5,6,7,8,9,10/
	data gd/1,2,3,4,5,6,7,8,9,10/
	data hd/1,2,3,4,5,6,7,8,9,10/
	data pd/1,2,3,4,5,6,7,8,9,10/
	data qd/1,2,3,4,5,6,7,8,9,10/
	data rd/1,2,3,4,5,6,7,8,9,10/
	data sd/1,2,3,4,5,6,7,8,9,10/
	data td/1,2,3,4,5,6,7,8,9,10/
	data ud/1,2,3,4,5,6,7,8,9,10/
	data vd/1,2,3,4,5,6,7,8,9,10/
	data wd/1,2,3,4,5,6,7,8,9,10/
	data xd/1,2,3,4,5,6,7,8,9,10/
	data yd/1,2,3,4,5,6,7,8,9,10/
	data zd/1,2,3,4,5,6,7,8,9,10/
	data zab/10*7/
	data zac/10*8/
	data zad/10*9/
	data iexp/
     + 220, 55, 1759, 55, 55, 220, 55, 1759, 55, 55,
     + 220, 55, 1759, 55, 55, 220, 55, 1759, 55, 55,
     + 220, 55, 1759, 55, 55, 220, 55, 1759, 55, 55,
     + 220, 55, 1759, 55, 55, 220, 55, 1588, 55, 220,
     + 55, 1759, 55, 55, 220, 55, 1759, 55, 55, 220,
     + 55, 1759, 55, 55, 220, 55, 1807, 55, 220, 55,
     + 1759, 55, 55, 220, 55, 1759, 55, 55, 220, 55,
     + 1759, 55, 55, 220, 55, 2026, 55, 70, 80, 90
     + /
	end


text:
! MIT License
! Copyright (c) 2020 Sam Miller
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

submodule(mod_field) field_2d_gpu_operators
!< Summary: This is the submodule to the field_2d_t type that contains
!<          all of the GPU arithmetic operators

contains

! --------------------------------------------------------------------
! Reduction operators (min/max/sum val and loc)
! --------------------------------------------------------------------
  module function field_maxval_gpu(f) result(res)
    class(field_2d_t), intent(in) :: f
    real(rk) :: res
  endfunction field_maxval_gpu

  module function field_maxloc_gpu(f) result(res)
    class(field_2d_t), intent(in) :: f
    integer(ik), dimension(2) :: res
  endfunction field_maxloc_gpu

  module function field_minval_gpu(f) result(res)
    class(field_2d_t), intent(in) :: f
    real(rk) :: res
  endfunction field_minval_gpu

  module function field_minloc_gpu(f) result(res)
    class(field_2d_t), intent(in) :: f
    integer(ik), dimension(2) :: res
  endfunction field_minloc_gpu

  module function field_sum_gpu(f) result(res)
    class(field_2d_t), intent(in) :: f
    real(rk) :: res
  endfunction field_sum_gpu

! --------------------------------------------------------------------
! Arithmetic operators
! --------------------------------------------------------------------
  module subroutine field_add_field_gpu(lhs, f, res)
    !< Implementation of the field_2d_t + field_2d_t operation
    class(field_2d_t), intent(in) :: lhs !< left-hand-side of the operation
    class(field_2d_t), intent(in) :: f
    type(field_2d_t), intent(inout) :: res
  endsubroutine field_add_field_gpu

  module subroutine field_sub_field_gpu(lhs, f, res)
    !< Implementation of the field_2d_t + field_2d_t operation
    class(field_2d_t), intent(in) :: lhs !< left-hand-side of the operation
    class(field_2d_t), intent(in) :: f
    type(field_2d_t), intent(inout) :: res
  endsubroutine field_sub_field_gpu

  module subroutine field_mul_field_gpu(lhs, f, res)
    !< Implementation of the field_2d_t * field_2d_t operation
    class(field_2d_t), intent(in) :: lhs !< left-hand-side of the operation
    class(field_2d_t), intent(in) :: f
    type(field_2d_t), intent(inout) :: res
  endsubroutine field_mul_field_gpu

  module subroutine field_div_field_gpu(lhs, f, res)
    !< Implementation of the field_2d_t * field_2d_t operation
    class(field_2d_t), intent(in) :: lhs !< left-hand-side of the operation
    class(field_2d_t), intent(in) :: f
    type(field_2d_t), intent(inout) :: res
  endsubroutine field_div_field_gpu

  module subroutine field_add_real_1d_gpu(lhs, x, res)
    !< Implementation of the field_2d_t + real64 operation
    class(field_2d_t), intent(in) :: lhs !< left-hand-side of the operation
    real(rk), intent(in) :: x
    type(field_2d_t), intent(inout) :: res
  endsubroutine field_add_real_1d_gpu

  module subroutine field_add_real_2d_gpu(lhs, x, res)
    !< Implementation of the field_2d_t + real64 operation
    class(field_2d_t), intent(in) :: lhs !< left-hand-side of the operation
    real(rk), dimension(lhs%ilo:, lhs%jlo:), intent(in) :: x
    type(field_2d_t), intent(inout) :: res
  endsubroutine field_add_real_2d_gpu

  module subroutine field_sub_real_1d_gpu(lhs, x, res)
    !< Implementation of the field_2d_t + real64 operation
    class(field_2d_t), intent(in) :: lhs !< left-hand-side of the operation
    real(rk), intent(in) :: x
    type(field_2d_t), intent(inout) :: res
  endsubroutine field_sub_real_1d_gpu

  module subroutine field_sub_real_2d_gpu(lhs, x, res)
    !< Implementation of the field_2d_t + real64 operation
    class(field_2d_t), intent(in) :: lhs !< left-hand-side of the operation
    real(rk), dimension(lhs%ilo:, lhs%jlo:), intent(in) :: x
    type(field_2d_t), intent(inout) :: res
  endsubroutine field_sub_real_2d_gpu

  module subroutine real_1d_sub_field_gpu(x, rhs, res)
    !< Implementation of the field_2d_t + real64 operation
    class(field_2d_t), intent(in) :: rhs !< right-hand-side of the operation
    real(rk), intent(in) :: x
    type(field_2d_t), intent(inout) :: res
  endsubroutine real_1d_sub_field_gpu

  module subroutine real_2d_sub_field_gpu(x, rhs, res)
    !< Implementation of the field_2d_t + real64 operation
    class(field_2d_t), intent(in) :: rhs !< right-hand-side of the operation
    real(rk), dimension(rhs%ilo:, rhs%jlo:), intent(in) :: x
    type(field_2d_t), intent(inout) :: res
  endsubroutine real_2d_sub_field_gpu

  module subroutine field_div_real_1d_gpu(lhs, x, res)
    !< Implementation of the field_2d_t / real64 operation
    class(field_2d_t), intent(in) :: lhs !< left-hand-side of the operation
    real(rk), intent(in) :: x
    type(field_2d_t), intent(inout) :: res
  endsubroutine field_div_real_1d_gpu

  module subroutine field_div_real_2d_gpu(lhs, x, res)
    !< Implementation of the field_2d_t / real64 operation
    class(field_2d_t), intent(in) :: lhs !< left-hand-side of the operation
    real(rk), dimension(lhs%ilo:, lhs%jlo:), intent(in) :: x
    type(field_2d_t), intent(inout) :: res
  endsubroutine field_div_real_2d_gpu

  module subroutine real_1d_div_field_gpu(x, rhs, res)
    !< Implementation of the field_2d_t / real64 operation
    class(field_2d_t), intent(in) :: rhs !< right-hand-side of the operation
    real(rk), intent(in) :: x
    type(field_2d_t), intent(inout) :: res
  endsubroutine real_1d_div_field_gpu

  module subroutine real_2d_div_field_gpu(x, rhs, res)
    !< Implementation of the field_2d_t / real64 operation
    class(field_2d_t), intent(in) :: rhs !< right-hand-side of the operation
    real(rk), dimension(rhs%ilo:, rhs%jlo:), intent(in) :: x
    type(field_2d_t), intent(inout) :: res
  endsubroutine real_2d_div_field_gpu

  module subroutine field_mul_real_2d_gpu(lhs, x, res)
    !< Implementation of the field_2d_t * array operation
    class(field_2d_t), intent(in) :: lhs !< left-hand-side of the operation
    real(rk), dimension(lhs%ilo:, lhs%jlo:), intent(in) :: x
    type(field_2d_t), intent(inout) :: res
  endsubroutine field_mul_real_2d_gpu

  module subroutine field_mul_real_1d_gpu(lhs, x, res)
    !< Implementation of the field_2d_t * real64 operation
    class(field_2d_t), intent(in) :: lhs !< left-hand-side of the operation
    real(rk), intent(in) :: x
    type(field_2d_t), intent(inout) :: res
  endsubroutine field_mul_real_1d_gpu
endsubmodule

max_issues_repo_name: smillerc/cato

max_forks_repo_name: smillerc/cato


text:
      subroutine conv(omega,xsec,ne1,vicorr)
c     multiply xsec by theta(omega-efermi) and
c     convolute xsec(omega) with  xloss/((omega-omega0)**2+xloss**2)/pi
c     the result is xsec0(omega0)

      implicit double precision (a-h, o-z)
      include '../HEADERS/const.h'
      include '../HEADERS/dim.h'
      dimension  omega(nex)
      complex*16 xsec(nex), xsec0(nex), xsecdx

      complex*16 conv1
      external conv1

      do ie = 1,ne1
         xsec0(ie) = 0.0d0
         omega0 = omega(ie)
c        Add one more point to correct for the finite grid
c        at large energies. Use linear interpolation.
         dx = max( omega(ne1) - omega(ne1-1), 50*vicorr)
         xlast = omega(ne1)+dx
         dx = dx / ( omega(ne1) - omega(ne1-1))
         xsecdx = xsec(ne1)+ (xsec(ne1)-xsec(ne1-1)) * dx

c        first interval
         do i = 1, ne1-1
            xsec0(ie) = xsec0(ie) +
     1      conv1(omega(i),omega(i+1),xsec(i),xsec(i+1),omega0,vicorr)
         enddo
c        last interval
         xsec0(ie) = xsec0(ie) +
     1   conv1(omega(ne1),xlast,xsec(ne1),xsecdx,omega0,vicorr)
         xsec0(ie) = xsec0(ie) /real(pi)
      enddo
      do ie = 1, ne1
         xsec(ie) = xsec0(ie)
      enddo
      return
      end

      complex*16 function conv1(x1,x2,y1,y2,x0,xloss)
c     convolution of function 1/(omega-omega0-i*xloss)/pi
c     makes linear interpolation for function between x1,x2 and
c     takes advantage that the integral can be taken analytically.
      implicit double precision (a-h, o-z)
      complex*16  y1, y2, t, coni,dum, a, b
      parameter (coni = (0.d0,1.d0))

      d = (x2-x1) / 2.d0
      a = dble(y2-y1) / 2.d0
      b = dble(y2+y1) / 2.d0
      t = d / ( (x1+x2)/2 - x0 - coni*xloss )
      if (abs(t) .ge. 0.1) then
         dum = 2.d0*a + (b - a/t) * log((1+t)/(1-t))
      else
         dum = 2.d0*b*(t+t**3 / 3.d0) - 2.d0/3.d0 * a*t**2
      endif
      conv1 = dimag (dum)

      d = (x2-x1) / 2.d0
      a = dimag(y2-y1) / 2.d0
      b = dimag(y2+y1) / 2.d0
      t = d / ( (x1+x2)/2 - x0 - coni*xloss )
      if (abs(t) .ge. 0.1) then
         dum = 2.d0*a + (b - a/t) * log((1+t)/(1-t))
      else
         dum = 2.d0*b*(t+t**3 / 3.d0) - 2.d0/3.d0 * a*t**2
      endif
      conv1 = conv1 + coni* dimag( dum)

      return
      end


text:
*
*
c      subroutine newpos(k,nup,rnew1,rnew2,iruntemp)
      subroutine newpos(k,nup)
*
*
*       determine new positions of two interacting objects after encounter
*       ------------------------------------------------------------------
*       M. HENON 1971, Astrophysics and Space Science Vol. 14, 151
*       ----------------------------------------------------------
*
*
      include 'common.h'
*
*
      real*8 e1,e2,a1,a2,q1,dz,b1,c1,gmin1,gmax1,gmin2,gmax2,
     &       rmin1,rmin2,rmax1,rmax2,rmaxz,s,u1,u2,v1,drds,q0,
     &       rtot,rles,xcut,ycut,xrr,rnn1,rnn2,ytot,qmax,rnew1,
     &       rnew2
c     &       rtot,rles,xcut,ycut,xrr,vrp,vtp,rnn1,rnn2
*
      real*4 ran2
*
      integer k,im1,im2,n,i,ipoin1,ipoin2,irm1,irm2,ipp,nup,m,imi1,
     &        imi2,ipo,kmin1,kmax1,kmin2,kmax2,imk1,imk2,maxk,
     &        mmin,mmax
*
c      common /timep/ vrp(nmax),vtp(nmax)
*
*
      n = nt
      i = 0
      im1 = iname(k)
      im2 = iname(k+1)
      gmin1 = 0.0d0
      gmin2 = 0.0d0
      gmax1 = 0.0d0
      gmax2 = 0.0d0
      u1 = 0.d0
      u2 = 0.d0
      v1 = 0.d0
      q0 = 0.d0
      ipoin1 = 0
      ipoin2 = 0
      irm1 = 0
      irm2 = 0
      imi1 = 0
      imi2 = 0
      imk1 = 0
      imk2 = 0
c      xmin(im1) = 0.0d0
c      xmin(im2) = 0.0d0
c      xmax(im1) = 0.0d0
c      xmax(im2) = 0.0d0
      xmaxz(im1) = 0.0d0
      xmaxz(im2) = 0.0d0
      xgmin(im1) = 0.0d0
      xgmin(im2) = 0.0d0
      xgmax(im1) = 0.0d0
      xgmax(im2) = 0.0d0
      rnew1 = r(k)
      rnew2 = r(k+1)
      vrp(im1) = vr(im1)
      vtp(im1) = vt(im1)
      vrp(im2) = vr(im2)
      vtp(im2) = vt(im2)
*
*       find maximal distance for new position of objects in super-zone
*
      ipp = nzst(nup)
      maxk = ipp + 1

*
      if(ipp.eq.n) then
        rmaxz = 1.d12
        maxk = n
      else
        rmaxz = 0.01d0*r(ipp) + 0.99d0*r(ipp+1)
      endif
*
*       determine total energy and angular momentum for interacting objects
*
      e1 = u(k) + 0.5d0*(vr(im1)**2 + vt(im1)**2)
      e2 = u(k+1) + 0.5d0*(vr(im2)**2 + vt(im2)**2)
      a1 = r(k)*vt(im1)
      a2 = r(k+1)*vt(im2)
*
*       check for escapers  -   energy > 0
*       do not find a new possition for massles particles
*
      if(e1.gt.0.0d0.or.body(im1).eq.0.d0) then

*
c          xmin(im1) = rnew1
c          xmax(im1) = rnew1
          xmaxz(im1) = rmaxz
          ipoin1 = 1
*
       endif
*
*       do not find a new possition for massles particles
*
      if(e2.gt.0.0d0.or.body(im2).eq.0.d0) then
*
c          xmin(im2) = rnew2
c          xmax(im2) = rnew2
          xmaxz(im2) = rmaxz
          ipoin2 = 1
*
      endif
*
*
*       find rmin and rmax for interacting objects
*
*       first star
*
      if(ipoin1.eq.1) go to 30
*
cThis call attempts to find m for a simple case (the majority).  Otherwise it returns
czero, and then the standard procedure is followed.
      m = mmin(k,e1,a1,n)
      if (m.gt.0) goto 13
      i = -1
   10 i = i + 1
      m = k - i
*
      if(m.eq.0) then
        b1 = 0.0d0
        c1 = u(1)
        go to 15
      endif
*
      q1 = 2.0d0*(e1 - u(m)) - a1**2/r(m)**2
      if(q1.ge.0.0d0) go to 10
*
      if(m.eq.k) then
        imk1 = 1
        rmin1 = 0.9999d0*r(k)
        write(6,*) '< r(k)   k,im1,r(k),u(k),e1,vr,vt,body = ',
     &  k,im1,r(k),u(k),e1,vr(im1),vt(im1),body(im1)
	call flush(6)
        kmin1 = 1
        go to 17
      endif
*
      if(m.eq.n) then
        b1 = u(m)*r(m)
        c1 = 0.0d0
        go to 15
      endif
*
 13   continue
      dz = 1.0d0/r(m+1) - 1.0d0/r(m)
      b1 = (u(m+1) - u(m))/dz
      c1 = (u(m)/r(m+1) - u(m+1)/r(m))/dz
   15 rmin1 = (-b1 + sqrt(b1*b1 - 2.0d0*a1*a1*(c1 - e1)))/a1/a1
      rmin1 = 1.0d0/rmin1
*
      if(rmin1.gt.r(k)) then
*
        xrr = 1.d0 - rmin1/r(k)
        xrr = abs(xrr)
*
        if(xrr.gt.1.0d-6) then
          write(6,*) '    rmin1  >  r(k)  iseed = ',iseed
          write(6,*) 'm,k,rmin1,r(k),r(m),r(m-1),u(m),u(m-1),e1,a1 =',
     &                m,k,rmin1,r(k),r(m),r(m-1),u(m),u(m-1),e1,a1
          call flush(6)
*
          rmin1 = 0.9999d0*r(k)
          imi1 = 1
        else
          rmin1 = r(k)
          imi1 = 1
        endif
*
      endif
*
      gmin1 = 2.0d0*a1*a1/rmin1**3 + 2.0d0*b1/rmin1**2
c      print*,gmin1,rmin1
*
      if(gmin1.lt.0.0d0) then
*
        write(6,*) '    gmin1  <   0.0   iseed = ',iseed
        write(6,*) 'm,k,rm1,r(k),r(m),r(m-1),u(m),u(m-1),e1,a1,b1,gm1=',
     &            m,k,rmin1,r(k),r(m),r(m-1),u(m),u(m-1),e1,a1,b1,gmin1

        call flush(6)
*
        gmin1 = 1.d10
      endif
*
      if(m.le.1) then
        kmin1 = 1
      else
        kmin1 = m - 1
      endif
*
      if((m.eq.n).or.(imi1.eq.1)) kmin1 = 1
*
 17   continue
*
      m = mmax(k,e1,a1,n)
      if (m.gt.0) goto 23

      i = -1
   20 i = i + 1
      m = k + i
*
      if(m.gt.n) then
        b1 = u(m-1)*r(m-1)
        c1 = 0.0d0
        go to 25
      endif
*
      q1 = 2.0d0*(e1 - u(m)) - a1**2/r(m)**2
c      print*,i,m,q1,e1,u(m),a1,r(m)
      if(q1.ge.0.0d0) go to 20
*
      if(m.eq.k) then
        if(imk1.eq.1) then
          rnew1 = r(k)
          imk1 = 2
          write(6,*) '-- O --  k,im1,r(k),u(k),e1= ',k,im1,r(k),u(k),e1
	  call flush(6)
          rmin1 = rnew1
          rmax1 = rnew1
          go to 30
        endif
*
        imk1 = 1
        rmax1 = 1.0001d0*r(k)
        write(6,*) '> r(k)    k,im1,r(k),u(k),e1 = ',k,im1,r(k),u(k),e1
	call flush(6)
        kmax1 = maxk
        go to 30
      endif
*
      if(m.eq.1) then
*
        if(imi1.eq.1) then
*
        write(6,*) ' m = 1   rnew1 = r(k)  iseed = ',iseed
        write(6,*) 'm,k,r(k),r(m),u(m),e1,a1 =',
     &              m,k,r(k),r(m),u(m),e1,a1
        call flush(6)
*
          imi1 = 2
          rnew1 = r(k)
          go to 30
        endif
*
        write(6,*) ' m = 1   rmax1 = 1.0001*r(k)   iseed = ',iseed
        write(6,*) 'm,k,r(k),r(m),u(m),e1,a1 =',
     &              m,k,r(k),r(m),u(m),e1,a1
        call flush(6)
*
        rmax1 = 1.0001d0*r(k)
        gmax1 = -1.d10
        go to 27
      endif
*
 23   continue
      dz = 1.0d0/r(m) - 1.0d0/r(m-1)
      b1 = (u(m) - u(m-1))/dz
      c1 = (u(m-1)/r(m) - u(m)/r(m-1))/dz
   25 continue
c      print*,c1,u(m-1),r(m),u(m),r(m-1),dz,m,k
      rmax1 = 0.5d0*(-b1 + sqrt(b1*b1 - 2.0d0*a1*a1*(c1 - e1)))
     &        /(c1 - e1)
*
      if(rmax1.lt.r(k)) then
*
        if(imi1.eq.1) then
          write(6,*) 'circular orbit     iseed,k = ',iseed,k
	  call flush(6)
          imi1 = 2
          rnew1 = r(k)
          rmin1 = rnew1
          rmax1 = rnew1
          go to 30
        endif
*
        xrr = 1.d0 - rmax1/r(k)
        xrr = abs(xrr)
*
        if(xrr.gt.1.0d-6) then
          write(6,*) ' rmax1 < r(k)  iseed =',iseed
          write(6,*) 'm,k,r(k),r(m),r(m-1),u(m),u(m-1),e1,a1 =',
     &                m,k,r(k),r(m),r(m-1),u(m),u(m-1),e1,a1
          call flush(6)
*
          rmax1 = 1.0001d0*r(k)
          imi1 = 3
        else
          rmax1 = r(k)
          imi1 = 3
        endif
*
      endif
*
      gmax1 = 2.0d0*a1*a1/rmax1**3 + 2.0d0*b1/rmax1**2
*
      if(gmax1.gt.0.0d0) then
*
        write(6,*) ' gmax1  >  0.0   iseed,gmax1 = ',iseed,gmax1
        write(6,*) 'm,k,rmax1,r(k),r(m),r(m-1),u(m),u(m-1),e1,a1 =',
     &              m,k,rmax1,r(k),r(m),r(m-1),u(m),u(m-1),e1,a1
        call flush(6)
*
        gmax1 = -1.d10
      endif
*
  27  if(rmax1.lt.rmin1) then
*
        write(6,*) 'rmax1 < rmin1  iseed = ',iseed
        write(6,*) 'r(k),a1,e1,rmax1,rmin1,m = ',
     &             r(k),a1,e1,rmax1,rmin1,m
        call flush(6)
*
      endif
*
      gmax1 = sqrt(-3.0d0*(rmax1 - rmin1)/gmax1)
      gmin1 = sqrt(3.0d0*(rmax1 - rmin1)/gmin1)
*
      if(rmax1.gt.rmaxz) irm1 = 1
*
      if(m.gt.n) then
        kmax1 = n
      else
        kmax1 = m + 1
      endif
*
      if((m.eq.1).or.(imi1.eq.3)) kmax1 = maxk
*
C       print*,ipoin1,k,e1,a1,maxk,rmaxz,gmax1,gmin1,istar,rmin1,
C      &     rmax1,imi1,imk1,kmax1,kmin1,irm1
c      stop
*       second star
*
   30 continue
*
      xgmin(im1) = gmin1
      xgmax(im1) = gmax1
*
      if(ipoin2.eq.1) go to 60
*
      m = mmin(k+1,e2,a2,n)
      if (m.gt.0) goto 43
      i = -1
   40 i = i + 1
      m = k + 1 - i
*
      if(m.eq.0) then
        b1 = 0.0d0
        c1 = u(1)
        go to 45
      endif
*
      q1 = 2.0d0*(e2 - u(m)) - a2**2/r(m)**2
      if(q1.ge.0.0d0) go to 40
*
      if(m.eq.k+1) then
        imk2 = 1
        rmin2 = 0.9999d0*r(k+1)
        write(6,*) '< r(k+1) k+1,im2,r,u,e2,vr,vt=',
     &  k+1,im2,r(k+1),u(k+1),e2,vr(im2),vt(im2)
	call flush(6)
        kmin2 = 1
        go to 47
      endif
*
      if(m.eq.n) then
        b1 = u(m)*r(m)
        c1 = 0.0d0
        go to 45
      endif
*
 43   continue
      dz = 1.0d0/r(m+1) - 1.0d0/r(m)
      b1 = (u(m+1) - u(m))/dz
      c1 = (u(m)/r(m+1) - u(m+1)/r(m))/dz
   45 rmin2 = (-b1 + sqrt(b1*b1 - 2.0d0*a2*a2*(c1 - e2)))/a2/a2
      rmin2 = 1.0d0/rmin2
*
      if(rmin2.gt.r(k+1)) then
*
        xrr = 1.d0 - rmin2/r(k+1)
        xrr = abs(xrr)
*
        if(xrr.gt.1.0d-6) then
          write(6,*) '    rmin2  >  r(k+1)  iseed = ',iseed
          write(6,*) 'm,k1,rm2,r(k+1),r(m),r(m+1),u(m),u(m+1),e2,a2 =',
     &                m,k+1,rmin2,r(k+1),r(m),r(m+1),u(m),u(m+1),e2,a2
          call flush(6)
*
          rmin2 = 0.9999d0*r(k+1)
          imi2 = 1
        else
          rmin2 = r(k+1)
          imi2 = 1
        endif
*
      endif
*
      gmin2 = 2.0d0*a2*a2/rmin2**3 + 2.0d0*b1/rmin2**2
*
      if(gmin2.lt.0.0d0) then
*
        write(6,*) '    gmin2 < 0.0   iseed,gmin2 = ',iseed,gmin2
        write(6,*) 'm,k1,rmin2,r(k+1),r(m),r(m+1),u(m),u(m+1),e2,a2 =',
     &              m,k+1,rmin2,r(k+1),r(m),r(m+1),u(m),u(m+1),e2,a2
        call flush(6)
*
        gmin2 = 1.d10
      endif
*
      if(m.le.1) then
        kmin2 = 1
      else
        kmin2 = m - 1
      endif
*
      if((m.eq.n).or.(imi2.eq.1)) kmin2 = 1
*
 47   continue
*
      m = mmax(k+1,e2,a2,n)
      if (m.gt.0) go to 53

      i = -1
*
   50 i = i + 1
      m = k + 1 + i
*
      if(m.gt.n) then
        b1 = u(m-1)*r(m-1)
        c1 = 0.0d0
        go to 55
      endif
*
      q1 = 2.0d0*(e2 - u(m)) - a2**2/r(m)**2
      if(q1.ge.0.0d0) go to 50
*
      if(m.eq.k+1) then
        if(imk2.eq.1) then
          rnew2 = r(k+1)
          imk2 = 2
          write(6,*) '-- O -- k+1,im2,r,u,e2=',k+1,im2,r(k+1),u(k+1),e2
	  call flush(6)
          rmin2 = rnew2
          rmax2 = rnew2
          go to 60
        endif

        imk2 = 1
        rmax2 = 1.0001d0*r(k+1)
        write(6,*) '> r(k+1) k+1,im2,r,u,e2 = ',k+1,im2,r(k+1),u(k+1),e2
	call flush(6)
        kmax2 = maxk
        go to 60
      endif
*
 53   continue
      dz = 1.0d0/r(m) - 1.0d0/r(m-1)
      b1 = (u(m) - u(m-1))/dz
      c1 = (u(m-1)/r(m) - u(m)/r(m-1))/dz
   55 rmax2 = 0.5d0*(-b1 + sqrt(b1*b1 - 2.0d0*a2*a2*(c1 - e2)))
     &        /(c1 - e2)
*
      if(rmax2.lt.r(k+1)) then
*
        if(imi2.eq.1) then
*
          write(6,*) ' circular orbit   iseed,k+1 =',iseed,k+1
	  call flush(6)
          imi2 = 2
          rnew2 = r(k+1)
          rmin2 = r(k+1)
          rmax2 = r(k+1)
          go to 60
        endif
*
        xrr = 1.0d0 - rmax2/r(k+1)
        xrr = abs(xrr)
*
        if(xrr.gt.1.0d-6) then
          write(6,*) 'rmax2 < r(k+1), rmax2=1.0001*r(k+1),iseed= ',iseed
          write(6,*) 'm,k1,rmax2,r(k+1),r(m),r(m+1),u(m),u(m+1),e2,a2=',
     &              m,k+1,rmax2,r(k+1),r(m),r(m+1),u(m),u(m+1),e2,a2
          call flush(6)
*
          rmax2 = 1.0001d0*r(k+1)
          imi2 = 3
        else
          rmax2 = r(k+1)
          imi2 = 3
        endif
*
      endif
*
      gmax2 = 2.0d0*a2*a2/rmax2**3 + 2.0d0*b1/rmax2**2
*
      if(gmax2.gt.0.0d0) then
*
        write(6,*) '    gmax2  >  0.0   iseed,gmax2 = ',iseed,gmax2
        write(6,*) 'm,k1,rmax2,r(k+1),r(m),r(m-1),u(m),u(m-1),e2,a2 =',
     &              m,k+1,rmax2,r(k+1),r(m),r(m-1),u(m),u(m-1),e2,a2
        call flush(6)
*
        gmax2 = -1.d10
      endif
*
      if(rmax2.lt.rmin2) then
*
        write(6,*) 'rmax2 < rmin2  iseed = ',iseed
        write(6,*) 'r(k+1),a2,e2,rmax2,rmin2,m = ',
     &             r(k+1),a2,e2,rmax2,rmin2,m
        call flush(6)
      endif
*
      gmax2 = sqrt(-3.0d0*(rmax2 - rmin2)/gmax2)
      gmin2 = sqrt(3.0d0*(rmax2 - rmin2)/gmin2)
*
      if(m.gt.n) then
        kmax2 = n
      else
        kmax2 = m + 1
      endif
*
      if((m.eq.2).or.(imi2.eq.3)) kmax2 = maxk
*
      if(rmax2.gt.rmaxz) irm2 = 1
*
   60 continue
*
      xgmin(im2) = gmin2
      xgmax(im2) = gmax2
*
C       print*,ipoin2,k,e2,a2,maxk,rmaxz,gmax2,gmin2,istar,rmin2,
C      &     rmax2,imi2,imk2,kmax2,kmin2,irm2
c      stop

*       determination of the new position
*
*       first star
*
      if(ipoin1.eq.1) go to 80
      xmin(im1) = rmin1
      xmax(im1) = rmax1
      xmaxz(im1) = rmaxz
      if(imi1.eq.2) go to 80
      if(imk1.eq.2) go to 80
*
      if(kmax1-kmin1.lt.2) then
        kmax1 = n
        kmin1 = 1
      endif
*
      ipp = 0
      qmax = 1.9d0
      xcut = 0.01d0
      ycut = 1.d0 - xcut
      rles = rmin1
      rtot = rmax1
      if(irm1.eq.1) then
        rtot = rmaxz
        ycut = 1.0d0
      endif
*
   70 s = ran2(irun)
*
      if(s.lt.xcut) s = xcut
      if(s.gt.ycut) s = ycut
*
      s = -1.0d0 + 2.0d0*s
*
      rnew1 = 0.5d0*(rmin1 + rmax1) +
     &                    0.25d0*(rmax1 - rmin1)*(3.0d0*s - s**3)
*
      ipp = ipp + 1
*
      if(ipp.lt.10000) then
*
        if((rnew1.gt.rtot).or.(rnew1.lt.rmglob)) go to 70
        if((k.eq.1).and.(rnew1.lt.r(1))) go to 70
        if(ikind(im1).lt.0) then
          print*,'newpos-1: i,im1,ikind = ',i,im1,ikind(im1)
        endif
        if(ikind(im1).eq.2.and.rnew1.lt.r(1)) go to 70
c        if(rnew1.lt.r(1)) go to 70
*
        q0 = qmax*ran2(irun)*dmax1(gmin1,gmax1)
c        print*,'q0,s ',q0,s
        ipo = 0
        rnn1 = rnew1
        call potent(ipo,kmin1,kmax1,rnn1,u1)
        v1 = 2.0d0*(e1 - u1) - a1*a1/rnew1/rnew1
*
      else
*
        write(6,*) 'ipp>10000 isee,time,rmz = ',iseed,time,rmaxz
        write(6,*) 'k,im1,rx1,ri1,v1,q0,e1,u1,a1,rnew1 = ',
     &              k,im1,rmax1,rmin1,v1,q0,e1,u1,a1,rnew1
        call flush(6)
        rnew1 = r(k)
        vrp(im1) = vr(im1)
        vtp(im1) = vt(im1)
        go to 80
*
      endif
*
c
c      open(33,file='choose.dat',access='append')
c      write(33,1234) time,ipp,k,im1,xcut,ycut,s,rmin1,rnew1,
c     &               rmax1,rmaxz,e1,a1,u1,v1
c 1234 format(1x,'-1- ',1pe12.4,3i6,1p13e12.4)
c      close(33)
c
      if(v1.le.0.0d0) go to 70
      v1 = sqrt(v1)
      drds = 0.75d0*(rmax1 - rmin1)*(1.0d0 - s*s)
      q1 = drds/v1
*
      if(q1.gt.qmax*dmax1(gmin1,gmax1)) then
c        write(6,*) 'q1 > q  iseed,q1,gmin1,gmax1=',iseed,q1,gmin1,gmax1
        qmax = 2.0d0*qmax
        go to 70
      endif
*
c
c      open(33,file='choose.dat',access='append')
c      write(33,1235) time,ipp,k,im1,qmax,gmin1,gmax1,v1,drds,q0,q1
c 1235 format(1x,'-1a- ',1pe12.4,3i6,1p7e12.4)
c      close(33)
c
      if(q0.gt.q1) go to 70
*
*       check for a nearly radial orbit if the new position is not
*       pick up too close to the apocentre
*
      q1 = a1/rnew1
      vrp(im1) = v1
      vtp(im1) = q1
*
*
   80 continue
*
*       second star
*
      if(ipoin2.eq.1) go to 100
      xmin(im2) = rmin2
      xmax(im2) = rmax2
      xmaxz(im2) = rmaxz
      if(imi2.eq.2) go to 100
      if(imk2.eq.2) go to 100
*
      if(kmax2-kmin2.lt.2) then
        kmax2 = n
        kmin2 = 1
      endif
*
      ipp = 0
      qmax = 1.9d0
      xcut = 0.01d0
      ycut = 1.0d0 - xcut
      rles = rmin2
      rtot = rmax2
      if(irm2.eq.1) then
        rtot = rmaxz
        ytot = 1.0d0
      endif
*
   90 s = ran2(irun)
*
      if(s.lt.xcut) s = xcut
      if(s.gt.ycut) s = ycut
*
      s = -1.0d0 + 2.0d0*s
*
      rnew2 = 0.5d0*(rmin2 + rmax2) +
     &                       0.25d0*(rmax2 - rmin2)*(3.0d0*s - s**3)
*
      ipp = ipp + 1
*
      if(ipp.lt.10000) then
*
        if((rnew2.gt.rtot).or.(rnew2.lt.rmglob)) go to 90
        if((k+1.eq.2).and.(rnew2.lt.r(1))) go to 90
        if(ikind(im2).lt.0) then
          print*,'newpos-2: i,im2,ikind = ',i,im2,ikind(im2)
        endif
        if(ikind(im2).eq.2.and.rnew2.lt.r(1)) go to 90
c        if(rnew2.lt.r(1)) go to 90
*
        q0 = qmax*ran2(irun)*dmax1(gmin2,gmax2)
        ipo = 0
        rnn2 = rnew2
        call potent(ipo,kmin2,kmax2,rnn2,u2)
        v1 = 2.0d0*(e2 - u2) - a2*a2/rnew2/rnew2
*
      else
*
        write(6,*) 'ipp>10000 isee,t,rmz = ',iseed,time,rmaxz
        write(6,*) 'k+1,im2,rx2,ri2,v1,q0,e2,u2,a2,rnew2 = ',
     &              k+1,im2,rmax2,rmin2,v1,q0,e2,u2,a2,rnew2
        call flush(6)
        rnew2 = r(k+1)
        vrp(im2) = vr(im2)
        vtp(im2) = vt(im2)
        go to 100
*
      endif
*
c
c      open(33,file='choose.dat',access='append')
c      write(33,1236) time,ipp,k+1,im2,xcut,ycut,s,rmin2,rnew2,
c     &               rmax2,rmaxz,e2,a2,u2,v1
c 1236 format(1x,'-2- ',1pe12.4,3i6,1p13e12.4)
c      close(33)
c
      if(v1.le.0.0d0) go to 90
      v1 = sqrt(v1)
      drds = 0.75d0*(rmax2 - rmin2)*(1.0d0 - s*s)
      q1 = drds/v1
*
      if(q1.gt.qmax*dmax1(gmin2,gmax2)) then
c        write(6,*) 'q1 > q  iseed,q1,gmin2,gmax2 =',iseed,q1,gmin2,gmax2
        qmax = 2.0d0*qmax
        go to 90
      endif
*
c
c      open(33,file='choose.dat',access='append')
c      write(33,1237) time,ipp,k+1,im2,qmax,gmin2,gmax2,v1,drds,q0,q1
c 1237 format(1x,'-2a- ',1pe12.4,3i6,1p7e12.4)
c      close(33)
c
      if(q0.gt.q1) go to 90
*
*       check for a nearly radial orbit if the new position is not
*       pick up too close to the apocentre
*
      q1 = a2/rnew2
      vrp(im2) = v1
      vtp(im2) = q1
*
*
  100 continue
*
*
      rn1 = rnew1
      rn2 = rnew2
c      print*,rnew1,rnew2
c      stop
      return
*
      end
*
*
*
*
c      subroutine newton(x,r,r2,r1)
*
*       newton iteration for finding minimum and maximum values
*       -------------------------------------------------------
*       of random number s
*       ------------------
*
*
c      real*8 x,r,r2,r1,ri,f,fp,xs,xs1,dx
*
c      integer iter
*
*
c      iter = 0
c      xs = 0.5d0
*
c      ri = (r - r1)/(r2 - r1)
*
c 10   continue
*
c      f = (3.0d0 - 2.0d0*xs)*xs*xs - ri
c      fp = 6.0d0*xs*(1.0d0 - xs)
*
c      xs1 = xs - f/fp
*
c      dx = abs(xs1 - xs)
c      dx = dx/xs
c      xs = xs1
*
c      iter = iter + 1
c      if(iter.gt.1000) stop  'iter > 1000'
*
c      if(dx.gt.1.0d-6) go to 10
*
c      x = xs
*
c      return
c      end
*
*
*
      function mmin(k,e,a,n)
cThis function returns the m-value for rmin in simple case, otherwise zero
      include 'common.h'
      double precision e,a,qk,q1,qstar,q2
      integer k,n,k1,k2,kstar,nit,mmin
      q1 = 2.0d0*(e - u(1)) - a**2/r(1)**2
      qk = 2.0d0*(e - u(k)) - a**2/r(k)**2
      if(q1.lt.0.d0.and.qk.gt.0.d0.and.k.gt.1.and.k.le.n) then
         k1 = 1
         k2 = k
         q2 = qk
         nit = 1
 10      continue
         if (k2.eq.k1 + 1) then
            mmin = k1
            return
         else
            nit = nit + 1
            if (nit.gt.1000) then
               print*,'nit>1000',k,e,a,n,k1,k2,q1,q2
               mmin = 0
               return
            endif
            kstar = (k1+k2)/2
            qstar = 2.0d0*(e - u(kstar)) - a**2/r(kstar)**2
            if (qstar.lt.0.d0) then
               k1 = kstar
               q1 = qstar
               go to 10
            else if (qstar.gt.0.d0) then
               k2 = kstar
               q2 = qstar
               go to 10
            else
               mmin = 0
               print*,'qstar = 0',k,e,a,n,k1,k2,q1,q2
               return
            endif
         endif
      else
         print*,'non-standard',k,e,a,n,q1,qk
         mmin = 0
         return
      endif
      end

      function mmax(k,e,a,n)
cThis function returns the m-value for rmax in simple case, otherwise zero
cNote that it is the index of the star above the maximu
      include 'common.h'
      double precision e,a,qk,qn,qstar,q1,q2
      integer k,n,k1,k2,nit,kstar,mmax
      q1 = 0.d0
      q2 = 0.d0
      qk = 2.0d0*(e - u(k)) - a**2/r(k)**2
      qn = 2.0d0*(e - u(n)) - a**2/r(n)**2
      if(qn.lt.0.d0.and.qk.gt.0.d0.and.k.ge.1.and.k.lt.n) then
         k1 = k
         k2 = n
         q1 = qk
         q2 = qn
         nit = 1
 10      continue
         if (k2.eq.k1 + 1) then
            mmax = k2
            return
         else
            nit = nit + 1
            if (nit.gt.1000) then
               print*,'nit>1000',k,e,a,n,k1,k2,q1,q2
               mmax = 0
               return
            endif
            kstar = (k1+k2)/2
            qstar = 2.0d0*(e - u(kstar)) - a**2/r(kstar)**2
            if (qstar.gt.0.d0) then
               k1 = kstar
               q1 = qstar
               go to 10
            else if (qstar.lt.0.d0) then
               k2 = kstar
               q2 = qstar
               go to 10
            else
               mmax = 0
               print*,'qstar = 0',k,e,a,n,k1,k2,q1,q2
               return
            endif
         endif
      else
         print*,'non-standard',k,e,a,n,q1,qk
         mmax = 0
         return
      endif
      end


text:
C$Procedure  RTRIM ( Right trim )

      INTEGER FUNCTION RTRIM ( STRING )

C$ Abstract
C
C     Return the maximum of 1 and the location of the last non-blank
C     character in the 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     None.
C
C$ Keywords
C
C     CHARACTER
C     STRING
C
C$ Declarations

      CHARACTER*(*)         STRING

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C
C     STRING     I   String to be trimmed.
C
C     The function returns the maximum of 1 and the location of the
C     last non-blank character in STRING.
C
C$ Detailed_Input
C
C     STRING         is a string to be trimmed:  the location of the
C                    last non-blank character is desired.
C
C$ Detailed_Output
C
C     The function returns the maximum of 1 and the location of the
C     last non-blank character in STRING.
C
C     In particular, when STRING is blank, the function returns the
C     value 1.
C
C$ Parameters
C
C     None.
C
C$ Files
C
C     None.
C
C$ Exceptions
C
C     Error free.
C
C$ Particulars
C
C     When writing a character string to a file, we usually are content
C     to omit the trailing blanks.  We'd like to use LASTNB as an upper
C     substring bound, but we have to handle the case where LASTNB
C     returns 0, so we write:
C
C
C        WRITE ( UNIT, '(A)' ),  STRING ( : MAX (1, LASTNB (STRING)) )
C
C
C     This can be simplified using RTRIM:
C
C
C        WRITE ( UNIT, '(A)' ),  STRING ( : RTRIM (STRING) )  )
C
C
C     This routine has a counterpart, LTRIM, which finds the maximum of
C     1 and the position of the first non-blank character of a string.
C
C$ Examples
C
C     1)  Write the non-blank portion of each element of a character
C         cell to file SPUD.DAT:
C
C            DO I = 1,  CARDC (CELL)
C
C               CALL WRLINE ( 'SPUD.DAT',
C           .                  CELL(I) ( LTRIM (CELL) : RTRIM (CELL) ) )
C
C            END DO
C
C         When CELL(I) is blank, the string ' ' will be written.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.J. Bachman   (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.1, 10-MAR-1992 (WLT)
C
C        Comment section for permuted index source lines was added
C        following the header.
C
C-    SPICELIB Version 1.0.0, 02-MAY-1990  (NJB)
C
C-&

C$ Index_Entries
C
C     right trim
C
C-&

C
C     SPICELIB functions
C
      INTEGER               LASTNB

C
C     `Just do it'.
C
      RTRIM = MAX ( 1, LASTNB (STRING) )

      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:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2020 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine resultsk(nk,nactdok,vtu,solk,solt,ipompc,nodempc,
     &  coefmpc,nmpc)
!
!     calculates the turbulence correction (STEP 5) in the nodes
!
      implicit none
!
      integer ipompc(*),nodempc(3,*),nmpc,nk,nactdok(*),i,ist,
     &  node,ndir,index
!
      real*8 coefmpc(*),solk(*),vtu(2,*),fixed_dispk,fixed_dispt,
     &  solt(*)
!
!     extracting the pressure correction from the solution
!
      do i=1,nk
c         write(*,*) i,nactdok(i)
         if(nactdok(i).gt.0) then
            vtu(1,i)=solk(nactdok(i))
            vtu(2,i)=solt(nactdok(i))
c            write(*,*) i,solk(nactdok(i))
c            write(*,*) i,solt(nactdok(i))
         else
            vtu(1,i)=0.d0
            vtu(2,i)=0.d0
         endif
c         write(*,*) 'resultsk ',i,vtu(1,i),vtu(2,i)
      enddo
!
!     inserting the mpc information: it is assumed that the
!     temperature MPC's also apply to the turbulence
!
c      do i=1,nmpc
c         ist=ipompc(i)
c         node=nodempc(1,ist)
c         ndir=nodempc(2,ist)
c         if(ndir.ne.0) cycle
c         index=nodempc(3,ist)
c         fixed_dispk=0.d0
c         fixed_dispt=0.d0
c         if(index.ne.0) then
c            do
c               fixed_dispk=fixed_dispk-coefmpc(index)*
c     &              vtu(1,nodempc(1,index))
c               fixed_dispt=fixed_dispt-coefmpc(index)*
c     &              vtu(2,nodempc(1,index))
c               index=nodempc(3,index)
c               if(index.eq.0) exit
c            enddo
c         endif
c         fixed_dispk=fixed_dispk/coefmpc(ist)
c         vtu(1,node)=fixed_dispk
c         fixed_dispt=fixed_dispt/coefmpc(ist)
c         vtu(2,node)=fixed_dispt
c      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop