text:
!
! Copyright (c) Stanford University, The Regents of the University of
!               California, and others.
!
! All Rights Reserved.
!
! See Copyright-SimVascular.txt for additional details.
!
! 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.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
! IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
! PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
! OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!-----------------------------------------------------------------------
!
!     This routine embodies formulation for solving electrophysiology
!     model equation using operator-splitting method.
!
!-----------------------------------------------------------------------

!     This is for solving 3D electrophysiology diffusion equations
      SUBROUTINE CEP3D (eNoN, nFn, w, N, Nx, al, yl, dl, fN, lR, lK)
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: eNoN, nFn
      REAL(KIND=8), INTENT(IN) :: w, N(eNoN), Nx(3,eNoN), al(tDof,eNoN),
     2   yl(tDof,eNoN), dl(tDof,eNoN), fN(3,nFn)
      REAL(KIND=8), INTENT(INOUT) :: lR(1,eNoN), lK(1,eNoN,eNoN)

      INTEGER a, b, i
      REAL(KIND=8) :: T1, amd, wl, Diso, Dani(nFn), Vd, Vx(3), F(3,3),
     2   C(3,3), Jac, fl(3,nFn), Ls(nFn), D(3,3), DVx(3), DNx(3,eNoN)

      IF (nFn .LT. eq(cEq)%dmn(cDmn)%cep%nFn) err =
     2   "No. of anisotropic conductivies exceed mesh fibers"

      T1   = eq(cEq)%af*eq(cEq)%gam*dt
      amd  = eq(cEq)%am/T1
      wl   = w*T1

      Diso = eq(cEq)%dmn(cDmn)%cep%Diso
      DO i=1, nFn
         IF (i .LE. eq(cEq)%dmn(cDmn)%cep%nFn) THEN
            Dani(i) = eq(cEq)%dmn(cDmn)%cep%Dani(i)
         ELSE
            Dani(i) = Dani(i-1)
         END IF
      END DO

!     Compute the isotropic part of diffusion tensor based on spatial
!     isotropy for electromechanics. This models stretch induced changes
!     in conduction velocities
      Ls(:) = 1D0
      IF (cem%cpld) THEN
!        Get the displacement degrees of freedom
         DO a=1, nEq
            IF (eq(a)%phys .EQ. phys_struct .OR.
     2          eq(a)%phys .EQ. phys_vms_struct) THEN
               i = eq(a)%s
               EXIT
            END IF
         END DO

!        Compute deformation gradient tensor
         F(:,:) = 0D0
         F(1,1) = 1D0
         F(2,2) = 1D0
         F(3,3) = 1D0
         DO a=1, eNoN
            F(1,1) = F(1,1) + Nx(1,a)*dl(i,a)
            F(1,2) = F(1,2) + Nx(2,a)*dl(i,a)
            F(1,3) = F(1,3) + Nx(3,a)*dl(i,a)
            F(2,1) = F(2,1) + Nx(1,a)*dl(i+1,a)
            F(2,2) = F(2,2) + Nx(2,a)*dl(i+1,a)
            F(2,3) = F(2,3) + Nx(3,a)*dl(i+1,a)
            F(3,1) = F(3,1) + Nx(1,a)*dl(i+2,a)
            F(3,2) = F(3,2) + Nx(2,a)*dl(i+2,a)
            F(3,3) = F(3,3) + Nx(3,a)*dl(i+2,a)
         END DO
!        Jacobian
         Jac = MAT_DET(F, 3)

!        Compute Cauchy-Green tensor and its inverse
         C  = MATMUL(TRANSPOSE(F), F)
         C  = MAT_INV(C, 3)

!        Compute fiber stretch
         DO i=1, nFn
            Ls(i) = SQRT(NORM(fN(:,i), MATMUL(C, fN(:,i))))
            fl(:,i) = fN(:,i) / Ls(i)
         END DO
         IF (Ls(1) .LE. 1D0) Ls(1) = 1D0

!        Diffusion tensor - spatial isotropy
         Diso    = Diso * Jac
         Dani(:) = Dani(:) * Jac
         D(:,:)  = Diso * C(:,:)
      ELSE
         D(:,:)  = 0D0
         D(1,1)  = Diso
         D(2,2)  = Diso
         D(3,3)  = Diso
         fl(:,:) = fN(:,:)
      END IF

!     Compute anisotropic components of diffusion tensor
      DO i=1, nFn
         D(1,1) = D(1,1) + Dani(i)*fl(1,i)*fl(1,i)
         D(1,2) = D(1,2) + Dani(i)*fl(1,i)*fl(2,i)
         D(1,3) = D(1,3) + Dani(i)*fl(1,i)*fl(3,i)

         D(2,1) = D(2,1) + Dani(i)*fl(2,i)*fl(1,i)
         D(2,2) = D(2,2) + Dani(i)*fl(2,i)*fl(2,i)
         D(2,3) = D(2,3) + Dani(i)*fl(2,i)*fl(3,i)

         D(3,1) = D(3,1) + Dani(i)*fl(3,i)*fl(1,i)
         D(3,2) = D(3,2) + Dani(i)*fl(3,i)*fl(2,i)
         D(3,3) = D(3,3) + Dani(i)*fl(3,i)*fl(3,i)
      END DO

      i  = eq(cEq)%s
      Vd = 0D0
      Vx = 0D0
      DO a=1, eNoN
         Vd = Vd + N(a)*al(i,a)

         Vx(1) = Vx(1) + Nx(1,a)*yl(i,a)
         Vx(2) = Vx(2) + Nx(2,a)*yl(i,a)
         Vx(3) = Vx(3) + Nx(3,a)*yl(i,a)

         DNx(1,a) = D(1,1)*Nx(1,a) + D(1,2)*Nx(2,a) + D(1,3)*Nx(3,a)
         DNx(2,a) = D(2,1)*Nx(1,a) + D(2,2)*Nx(2,a) + D(2,3)*Nx(3,a)
         DNx(3,a) = D(3,1)*Nx(1,a) + D(3,2)*Nx(2,a) + D(3,3)*Nx(3,a)
      END DO

      DVx(1) = D(1,1)*Vx(1) + D(1,2)*Vx(2) + D(1,3)*Vx(3)
      DVx(2) = D(2,1)*Vx(1) + D(2,2)*Vx(2) + D(2,3)*Vx(3)
      DVx(3) = D(3,1)*Vx(1) + D(3,2)*Vx(2) + D(3,3)*Vx(3)

      DO a=1, eNoN
         lR(1,a) = lR(1,a) + w*(N(a)*Vd + Nx(1,a)*DVx(1) +
     2      Nx(2,a)*DVx(2) + Nx(3,a)*DVx(3))

         DO b=1, eNoN
            lK(1,a,b) = lK(1,a,b) + wl*(N(a)*N(b)*amd +
     2         Nx(1,a)*DNx(1,b) + Nx(2,a)*DNx(2,b) + Nx(3,a)*DNx(3,b))
         END DO
      END DO

      RETURN
      END SUBROUTINE CEP3D
!-----------------------------------------------------------------------
!     This is for solving 2D electrophysiology diffusion equation
      SUBROUTINE CEP2D (eNoN, nFn, w, N, Nx, al, yl, dl, fN, lR, lK)
      USE COMMOD
      USE ALLFUN
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: eNoN, nFn
      REAL(KIND=8), INTENT(IN) :: w, N(eNoN), Nx(2,eNoN), al(tDof,eNoN),
     2   yl(tDof,eNoN), dl(tDof,eNoN), fN(2,nFn)
      REAL(KIND=8), INTENT(INOUT) :: lR(1,eNoN), lK(1,eNoN,eNoN)

      INTEGER a, b, i
      REAL(KIND=8) :: T1, amd, wl, Diso, Dani(nFn), Vd, Vx(2), F(2,2),
     2   C(2,2), Jac, fl(2,nFn), Ls(nFn), D(2,2), DVx(2), DNx(2,eNoN)

      IF (nFn .LT. eq(cEq)%dmn(cDmn)%cep%nFn) err =
     2   "No. of anisotropic conductivies exceed mesh fibers"

      T1   = eq(cEq)%af*eq(cEq)%gam*dt
      amd  = eq(cEq)%am/T1
      wl   = w*T1

      Diso = eq(cEq)%dmn(cDmn)%cep%Diso
      DO i=1, nFn
         IF (i .LE. eq(cEq)%dmn(cDmn)%cep%nFn) THEN
            Dani(i) = eq(cEq)%dmn(cDmn)%cep%Dani(i)
         ELSE
            Dani(i) = Dani(i-1)
         END IF
      END DO

!     Compute the isotropic part of diffusion tensor based on spatial
!     isotropy for electromechanics. This models stretch induced changes
!     in conduction velocities
      Ls(:) = 1D0
      IF (cem%cpld) THEN
         DO a=1, nEq
            IF (eq(a)%phys .EQ. phys_struct .OR.
     2          eq(a)%phys .EQ. phys_vms_struct) THEN
               i = eq(a)%s
               EXIT
            END IF
         END DO

!        Compute deformation gradient tensor
         F(:,:) = 0D0
         F(1,1) = 1D0
         F(2,2) = 1D0
         DO a=1, eNoN
            F(1,1) = F(1,1) + Nx(1,a)*dl(i,a)
            F(1,2) = F(1,2) + Nx(2,a)*dl(i,a)
            F(2,1) = F(2,1) + Nx(1,a)*dl(i+1,a)
            F(2,2) = F(2,2) + Nx(2,a)*dl(i+1,a)
         END DO
!        Jacobian
         Jac = MAT_DET(F, 2)

!        Compute Cauchy-Green tensor and its inverse
         C  = MATMUL(TRANSPOSE(F), F)
         C  = MAT_INV(C, 2)

!        Compute fiber stretch
         DO i=1, nFn
            Ls(i) = SQRT(NORM(fN(:,i), MATMUL(C, fN(:,i))))
            fl(:,i) = fN(:,i) / Ls(i)
         END DO
         IF (Ls(1) .LE. 1D0) Ls(1) = 1D0

!        Diffusion tensor - spatial isotropy
         Diso    = Diso * Jac
         Dani(:) = Dani(:) * Jac
         D(:,:)  = Diso * C(:,:)
      ELSE
         D(:,:)  = 0D0
         D(1,1)  = Diso
         D(2,2)  = Diso
         fl(:,:) = fN(:,:)
      END IF

      DO i=1, nFn
         D(1,1) = D(1,1) + Dani(i)*fl(1,i)*fl(1,i)
         D(1,2) = D(1,2) + Dani(i)*fl(1,i)*fl(2,i)

         D(2,1) = D(2,1) + Dani(i)*fl(2,i)*fl(1,i)
         D(2,2) = D(2,2) + Dani(i)*fl(2,i)*fl(2,i)
      END DO

      i  = eq(cEq)%s
      Vd = 0D0
      Vx = 0D0
      DO a=1, eNoN
         Vd = Vd + N(a)*al(i,a)

         Vx(1) = Vx(1) + Nx(1,a)*yl(i,a)
         Vx(2) = Vx(2) + Nx(2,a)*yl(i,a)

         DNx(1,a) = D(1,1)*Nx(1,a) + D(1,2)*Nx(2,a)
         DNx(2,a) = D(2,1)*Nx(1,a) + D(2,2)*Nx(2,a)
      END DO

      DVx(1) = D(1,1)*Vx(1) + D(1,2)*Vx(2)
      DVx(2) = D(2,1)*Vx(1) + D(2,2)*Vx(2)

      DO a=1, eNoN
         lR(1,a) = lR(1,a) + w*(N(a)*Vd + Nx(1,a)*DVx(1) +
     2      Nx(2,a)*DVx(2))

         DO b=1, eNoN
            lK(1,a,b) = lK(1,a,b) + wl*(N(a)*N(b)*amd +
     2         Nx(1,a)*DNx(1,b) + Nx(2,a)*DNx(2,b))
         END DO
      END DO

      RETURN
      END SUBROUTINE CEP2D
!-----------------------------------------------------------------------
!     This is for solving 1D electrophysiology diffusion equation
!     for Purkinje fibers
      PURE SUBROUTINE CEP1D (eNoN, insd, w, N, Nx, al, yl, lR, lK)
      USE COMMOD
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: eNoN, insd
      REAL(KIND=8), INTENT(IN) :: w, N(eNoN), Nx(insd,eNoN),
     2   al(tDof,eNoN), yl(tDof,eNoN)
      REAL(KIND=8), INTENT(INOUT) :: lR(1,eNoN), lK(1,eNoN,eNoN)

      INTEGER a, b, i
      REAL(KIND=8) :: T1, amd, wl, Td, Tx, Diso, DNx(eNoN)

      T1   = eq(cEq)%af*eq(cEq)%gam*dt
      amd  = eq(cEq)%am/T1
      Diso = eq(cEq)%dmn(cDmn)%cep%Diso
      i    = eq(cEq)%s
      wl   = w*T1

      Td = 0D0
      Tx = 0D0
      DO a=1, eNoN
         Td = Td + N(a)*al(i,a)
         Tx = Tx + Nx(1,a)*yl(i,a)
         DNx(a) = Diso*Nx(1,a)
      END DO

      DO a=1, eNoN
         lR(1,a) = lR(1,a) + w*(N(a)*Td + Nx(1,a)*Diso*Tx)

         DO b=1, eNoN
            lK(1,a,b) = lK(1,a,b) + wl*(N(a)*N(b)*amd + Nx(1,a)*DNx(b))
         END DO
      END DO

      RETURN
      END SUBROUTINE CEP1D
!#######################################################################
      PURE SUBROUTINE BCEP (eNoN, w, N, h, lR)
      USE COMMOD
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: eNoN
      REAL(KIND=8), INTENT(IN) :: w, N(eNoN), h
      REAL(KIND=8), INTENT(INOUT) :: lR(dof,eNoN)

      INTEGER :: a
      REAL(KIND=8) f

      f = w*h

!     Here the loop is started for constructing left and right hand side
      DO a=1, eNoN
         lR(1,a) = lR(1,a) + N(a)*f
      END DO

      RETURN
      END SUBROUTINE BCEP
!#######################################################################

ext: f


program main

        use ISO_FORTRAN_ENV

implicit none

        ! Rearrangement using a + b + c = p; a^2 + b^2 = c^2 gives us;
        ! b = (p^2 - 2pa)/(2p - 2a)
        ! So if we iterate up to p <= 1000, and across all a; for a < b;
        ! the solution for b should be natural if there is a solution

        integer (int64), parameter :: MAXIMUM = 1000;
        integer (int64) :: a, p, num_sol, max_sol, largest;

        num_sol = 0
        largest = 0
        max_sol = 0

        OUTER: &
        do p = 2, MAXIMUM, 1

                num_sol = 0

                INNER: &
                do a = 2, p/2, 1

                        if(mod(p**2 - 2*p*a, (2*(p - a))) == 0) then
                                num_sol = num_sol + 1;
                        endif

                end do inner

                if(num_sol > max_sol) then
                        max_sol = num_sol;
                        largest = p;
                endif

        end do outer

        print*,(largest);

end program main


text:
Davis Well 27, according to the Davis Water Wells water well map.

I believe chlorination happens here (there is a shack that says Danger, 12.5% Liquid Chlorine). There is also a backup 300 KW diesel genset on site. Neat!

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! test_genalg.f90 --
!     Test for the genetic_algorithms module
!
!     Note:
!     The first step (find the minimum for f) does not work as
!     expected, the algorithm gets stuck on a suboptimum (at least
!     with one random sequence). Probably it is necessary to tune
!     the parameters to get a better performance
!     (Intended minimum: (-1.1, 0, 0, ... )
!
!     The second step does find the absolute minimum, but this is
!     a much simpler function.
!
module minimize_function
    use select_precision
contains

real(wp) function f( coords )
    real(wp), dimension(:), intent(in) :: coords

    real(wp)                           :: r

    r = sqrt( sum(coords**2) )
    f = (r**2 - 1.0) ** 2 - coords(1)

end function f

real(wp) function g( coords )
    real(wp), dimension(:), intent(in) :: coords

    g = sum(coords**2) - 1.0

end function g

end module minimize_function

program test_genalg
    use minimize_function
    use genetic_algorithms
    use select_precision

    implicit none

    type(GA_PARAMETERS)            :: params
    integer, parameter             :: ncoords = 10
    real(wp), dimension(2,ncoords) :: range
    real(wp), dimension(ncoords)   :: vector
    real(wp)                       :: value
    integer                        :: i

    !
    ! Set the parameters for the procedure
    !
    call set_parameters( params, update = .false., verbose = .true., &
             number_iterations = 1000 )

    !
    ! We can use the default procedure, so let's do that
    ! (Otherwise we would use the reversed communication example)
    !
    range(1,:) = -5.0
    range(2,:) =  5.0

    call find_minimum( params, range, vector, f, value )

    write(*,*) 'Minimum:         ',value
    write(*,*) 'Solution vector: ',vector
    write(*,*) 'f(vector):       ',f(vector)

    call find_minimum( params, range, vector, g, value )

    write(*,*) 'Minimum:         ',value
    write(*,*) 'Solution vector: ',vector
    write(*,*) 'g(vector):       ',g(vector)

end program

max_issues_repo_name: TerribleDev/Fortran-docker-mvc

max_forks_repo_name: TerribleDev/Fortran-docker-mvc


text:
yeeeeeeeeeeehaaaaaaaaaaaaaaaaaaaa

wooooooopieeeeeeeeeeee

euuuuuuuuunhhhhhhhhhhhh

solitarizationism

Signed up 200510

What it is? I could`nt understand
                                          Bol


module ps_nloc3_gth_module

  implicit none
  private
  public :: init_ps3_gth

contains

  subroutine init_ps3_gth( norb, no, lo, ro, hnl, gk, viodgk )
    implicit none
    integer,intent(in) :: norb(:)
    integer,intent(in) :: no(:,:)
    integer,intent(in) :: lo(:,:)
    real(8),intent(in) :: ro(:,:)
    real(8),intent(in) :: hnl(:,0:,:)
    real(8),intent(in) :: gk(:)
    real(8),intent(inout) :: viodgk(:,:,:)
    integer :: ielm,Nelm,iorb,l,i,mpgk,n
    real(8) :: q,pi,c

    pi = acos(-1.0d0)

    Nelm = size(norb)
    mpgk = size(viodgk,1)

    do ielm = 1, Nelm
      do iorb = 1, norb(ielm)

        n = no(iorb,ielm)
        l = lo(iorb,ielm)
        c = sqrt(hnl(n,l,ielm))

        select case( l )
        case( 0 )

          select case( n )
          case( 1 )
            do i = 1, mpgk
              q = gk(i)
              viodgk(i,iorb,ielm) = c*4.0d0*ro(iorb,ielm) &
              *sqrt(2.0d0*ro(iorb,ielm))*pi**1.25d0 &
              *exp(-0.5d0*(q*ro(iorb,ielm))**2)
            end do
          case( 2 )
            do i = 1, mpgk
              q = gk(i)
              viodgk(i,iorb,ielm) = c*8.0d0*ro(iorb,ielm) &
              *sqrt(2.0d0*ro(iorb,ielm)/15.0d0)*pi**1.25d0 &
              *exp(-0.5d0*(q*ro(iorb,ielm))**2) &
              *(3.0d0-(q*ro(iorb,ielm))**2)
            end do
          end select

        case( 1 )

          do i = 1, mpgk
            q = gk(i)
            viodgk(i,iorb,ielm) = c*8.0d0*ro(iorb,ielm)**2 &
            *sqrt(ro(iorb,ielm)/3.0d0)*pi**1.25d0 &
            *exp(-0.5d0*(q*ro(iorb,ielm))**2) &
            *q
          end do

        end select

      end do !iorb
    end do !ielm

  end subroutine init_ps3_gth

end module ps_nloc3_gth_module


text:
c
c   Dint – version 2.0  is licensed under the Apache License, Version 2.0 (the "License");
c   you may not use Dint – version 2.0 except in compliance with the License.
c   You may obtain a copy of the License at
c       http://www.apache.org/licenses/LICENSE-2.0
c   The license is also given in the LICENSE file.
c   Unless required by applicable law or agreed to in writing, software
c   distributed under the License is distributed on an "AS IS" BASIS,
c   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
c   See the License for the specific language governing permissions and limitations under the License.
c
c -------------------------------------------------------------------------------------------
c  Dint : Direct Nonadiabatic Trajectories A code for non-Born–Oppenheimer molecular dynamics
c
c  version 2.0
c
c  A. W. Jasper
c  Argonne National Laboratory
c
c  Rui Ming Zhang
c  Tsinghua University
c
c  and
c
c  D. G. Truhlar
c  University of Minnesota
c
c  copyright  2020
c  Donald G. Truhlar and Regents of the University of Minnesota
c----------------------------------------------------------------------------------------------

      subroutine ange(xx,pp,mm,natom,eig,bigj,bigjtot,erot,erottot)

c Compute rotational energies around principal axes

      implicit none
      include 'param.f'

      integer natom
      double precision xx(3,mnat),pp(3,mnat),mm(mnat),
     &   bigj(3),bigjtot,xxp(3,mnat),ppp(3,mnat)

c local
      integer i,j,info,ismall
      double precision erot(3),erottot,
     & temp1,temp2,
     & mom(3,3),ap(6),eig(3),rot(3,3),work(9),temp3,
     & mom2,rr

c single atom
      if (natom.eq.1) then
         bigj(1) = 0.d0
         bigj(2) = 0.d0
         bigj(3) = 0.d0
         bigjtot = 0.d0
         erot(1) = 0.d0
         erot(2) = 0.d0
         erot(3) = 0.d0
         erottot = 0.d0
         return
      endif

c I modified the polyatomic part to work for linears instead of this
c diatom
c      if (natom.eq.2) then
c      rr=0.d0
c      do j=1,3
c      rr=rr+(xx(j,1)-xx(j,2))**2
c      enddo
c      mom2=mm(1)*mm(2)/(mm(1)+mm(2))*rr   ! remember rr=R**2
c      write(6,*)"mom2=",mom2
c      return
c      endif

c natom >= 2
c compute moment of intertia matrix mom
      do i=1,3
      do j=1,3
          mom(i,j) = 0.d0
      enddo
      enddo

      do i=1,natom
         mom(1,1)=mom(1,1)+mm(i)*(xx(2,i)**2+xx(3,i)**2)
         mom(2,2)=mom(2,2)+mm(i)*(xx(1,i)**2+xx(3,i)**2)
         mom(3,3)=mom(3,3)+mm(i)*(xx(1,i)**2+xx(2,i)**2)
         mom(1,2)=mom(1,2)-mm(i)*(xx(1,i)*xx(2,i))
         mom(1,3)=mom(1,3)-mm(i)*(xx(1,i)*xx(3,i))
         mom(2,3)=mom(2,3)-mm(i)*(xx(2,i)*xx(3,i))
      enddo
      mom(2,1)=mom(1,2)
      mom(3,1)=mom(1,3)
      mom(3,2)=mom(2,3)

c      write(6,*)"mom"
c      write(6,*)(mom(1,j),j=1,3)
c      write(6,*)(mom(2,j),j=1,3)
c      write(6,*)(mom(3,j),j=1,3)

c      call angmom(xx,pp,mm,natom,bigj,bigjtot)
c      write(6,*)"j",bigj,bigjtot

c     diagonalize the mom matrix
      do i=1,3
      do j=i,3
        ap(i+(j-1)*j/2)=mom(i,j)
      enddo
      enddo
      call dspev( 'v','u',3,ap,eig,rot,3,work,info )

c      write(6,*)"rot"
c      write(6,*)(rot(1,j),j=1,3)
c      write(6,*)(rot(2,j),j=1,3)
c      write(6,*)(rot(3,j),j=1,3)
c      write(6,*)"eig"
c      write(6,*)(eig(j),j=1,3)

c rotate to diagonalize mom
      do i=1,natom
         temp1 = xx(1,i)
         temp2 = xx(2,i)
         temp3 = xx(3,i)
         xxp(1,i)=temp1*rot(1,1)+temp2*rot(2,1)+temp3*rot(3,1)
         xxp(2,i)=temp1*rot(1,2)+temp2*rot(2,2)+temp3*rot(3,2)
         xxp(3,i)=temp1*rot(1,3)+temp2*rot(2,3)+temp3*rot(3,3)
         temp1 = pp(1,i)
         temp2 = pp(2,i)
         temp3 = pp(3,i)
         ppp(1,i)=temp1*rot(1,1)+temp2*rot(2,1)+temp3*rot(3,1)
         ppp(2,i)=temp1*rot(1,2)+temp2*rot(2,2)+temp3*rot(3,2)
         ppp(3,i)=temp1*rot(1,3)+temp2*rot(2,3)+temp3*rot(3,3)
      enddo

      call angmom(xxp,ppp,mm,natom,bigj,bigjtot)

c rotational energies
c      write(6,*)"Ia = ",(eig(j),j=1,3)," au"
      erottot = 0.d0
      ismall=0
      do i=1,3
        if (eig(i).lt.1.d-10) then
          ismall=ismall+1
          erot(i) = 0.d0
        else
          erot(i) = bigj(i)**2/(2.d0*eig(i))
        endif
        erottot=erottot+erot(i)
      enddo
c      write(6,*)"Erot = ",(erot(j)*autoev,j=1,3)," eV"
c      write(6,*)"Ji = ",(bigj(j),j=1,3)
      if (ismall.eq.1) write(6,*)"Found & skipped 1 small Ia,",
     & " must be a linear molecule"
      if (ismall.gt.1) then
        write(6,*)"Found ",ismall," small Ia! Problem."
        write(6,*)"The values are:",(eig(j),j=1,3)," au"
        stop
      endif

c      write(6,*)"ange"
c      write(6,*)bigj(1),erot(1)*autoev,eig(1)
c      write(6,*)bigj(2),erot(2)*autoev,eig(2)
c      write(6,*)bigj(3),erot(3)*autoev,eig(3)

      return

      end

max_stars_repo_name: kevinmooreiii/DiNT

max_issues_repo_name: kevinmooreiii/DiNT

max_forks_repo_name: kevinmooreiii/DiNT


! { dg-do compile }
! { dg-options "-O -floop-nest-optimize" }
      SUBROUTINE ZTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, &
                   B, LDB )
      CHARACTER*1        SIDE, UPLO, TRANSA, DIAG
      INTEGER            M, N, LDA, LDB
      complex(kind((1.0d0,1.0d0)))         ALPHA
      complex(kind((1.0d0,1.0d0)))         A( LDA, * ), B( LDB, * )
      EXTERNAL           XERBLA
      INTRINSIC          CONJG, MAX
      LOGICAL            LSIDE, NOCONJ, NOUNIT, UPPER
      INTEGER            I, INFO, J, K, NROWA
      complex(kind((1.0d0,1.0d0)))         TEMP
      complex(kind((1.0d0,1.0d0)))         ONE
      PARAMETER        ( ONE  = ( 1.0D+0, 0.0D+0 ) )
      complex(kind((1.0d0,1.0d0)))         ZERO
      PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
      LSIDE  =  scan( SIDE  , 'Ll' )>0
      IF( LSIDE )THEN
         NROWA = M
      ELSE
         NROWA = N
      END IF
      NOCONJ =  scan( TRANSA, 'Tt' )>0
      NOUNIT =  scan( DIAG  , 'Nn' )>0
      UPPER  =  scan( UPLO  , 'Uu' )>0
      INFO   = 0
      IF( N.EQ.0 ) &
   RETURN
      IF( ALPHA.EQ.ZERO )THEN
         DO 20, J = 1, N
            DO 10, I = 1, M
               B( I, J ) = ZERO
   10       CONTINUE
   20    CONTINUE
         RETURN
      END IF
               DO 160, J = 1, N
                  DO 150, I = 1, M
                     TEMP = B( I, J )
                     IF( NOCONJ )THEN
                        IF( NOUNIT ) &
                     TEMP = TEMP*A( I, I )
                        DO 130, K = I + 1, M
                           TEMP = TEMP + A( K, I )*B( K, J )
  130                   CONTINUE
                     ELSE
                        IF( NOUNIT ) &
                     TEMP = TEMP*CONJG( A( I, I ) )
                        DO 140, K = I + 1, M
                           TEMP = TEMP + CONJG( A( K, I ) )*B( K, J )
  140                   CONTINUE
                     END IF
                     B( I, J ) = ALPHA*TEMP
  150             CONTINUE
  160          CONTINUE
      RETURN
      END


module funcs
use, intrinsic:: iso_fortran_env, only: sp=>real32, dp=>real64
implicit none (type, external)

interface addtwo
procedure addtwo_s, addtwo_d, addtwo_i
end interface addtwo

contains

elemental real(sp) function addtwo_s(x) result(y)
real(sp), intent(in) :: x

y = x + 2
end function addtwo_s

elemental real(dp) function addtwo_d(x) result(y)
real(dp), intent(in) :: x

y = x + 2
end function addtwo_d

elemental integer function addtwo_i(x) result(y)
integer, intent(in) :: x

y = x + 2
end function addtwo_i

end module funcs

program test2

use funcs
implicit none (type, external)

real(sp) :: twos = 2._sp
real(dp) :: twod = 2._dp
integer :: twoi = 2

print *, addtwo(twos), addtwo(twod), addtwo(twoi)

end program


text:
      subroutine desugar03(x, y, *, *)

C     Check desugaring of labelled DO loops with a forward reference

      read *, i

      if(i.lt.0) go to 400

c     labelled loop
 400  do 500 i = 1, n
         print *, i
 500  continue

      end


text:
module interpolate_module

  contains

      function interpolate(r, npts_model, model_r, model_var)

!     given the array of model coordinates (model_r), and variable (model_var),
!     find the value of model_var at point r using linear interpolation.
!     Eventually, we can do something fancier here.

      double precision, intent(in   ) :: r
      integer         , intent(in   ) :: npts_model
      double precision, intent(in   ) :: model_r(npts_model), model_var(npts_model)
      double precision                :: interpolate

      ! Local variables
      integer                         :: i, id
      double precision                :: slope,minvar,maxvar

!     find the location in the coordinate array where we want to interpolate
      do i = 1, npts_model
         if (model_r(i) .ge. r) exit
      enddo

      if (i .eq. npts_model+1) i = i-1
      id = i

      if (id .eq. 1) then

         slope = (model_var(id+1) - model_var(id))/(model_r(id+1) - model_r(id))
         interpolate = slope*(r - model_r(id)) + model_var(id)

         ! safety check to make sure interpolate lies within the bounding points
         minvar = min(model_var(id+1),model_var(id))
         maxvar = max(model_var(id+1),model_var(id))
         interpolate = max(interpolate,minvar)
         interpolate = min(interpolate,maxvar)

      else if (id .eq. npts_model) then

         slope = (model_var(id) - model_var(id-1))/(model_r(id) - model_r(id-1))
         interpolate = slope*(r - model_r(id)) + model_var(id)

         ! safety check to make sure interpolate lies within the bounding points
         minvar = min(model_var(id),model_var(id-1))
         maxvar = max(model_var(id),model_var(id-1))
         interpolate = max(interpolate,minvar)
         interpolate = min(interpolate,maxvar)

      else

         if (r .ge. model_r(id)) then

            slope = (model_var(id+1) - model_var(id))/(model_r(id+1) - model_r(id))
            interpolate = slope*(r - model_r(id)) + model_var(id)

            ! safety check to make sure interpolate lies within the bounding points
            minvar = min(model_var(id+1),model_var(id))
            maxvar = max(model_var(id+1),model_var(id))
            interpolate = max(interpolate,minvar)
            interpolate = min(interpolate,maxvar)

         else

            slope = (model_var(id) - model_var(id-1))/(model_r(id) - model_r(id-1))
            interpolate = slope*(r - model_r(id)) + model_var(id)

            ! safety check to make sure interpolate lies within the bounding points
            minvar = min(model_var(id),model_var(id-1))
            maxvar = max(model_var(id),model_var(id-1))
            interpolate = max(interpolate,minvar)
            interpolate = min(interpolate,maxvar)

          end if

      endif

      end function interpolate

end module interpolate_module

max_stars_repo_name: memmett/BoxLib

max_issues_repo_name: memmett/BoxLib

max_forks_repo_name: memmett/BoxLib


text:
      SUBROUTINE PZBMATGEN( ICTXT, AFORM, AFORM2, BWL, BWU, N,
     $                     MB, NB, A,
     $                     LDA, IAROW, IACOL, ISEED,
     $                     MYROW, MYCOL, NPROW, NPCOL )
*
*
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 15, 1997
*
*     .. Scalar Arguments ..
*     .. Scalar Arguments ..
      CHARACTER*1        AFORM, AFORM2
      INTEGER            IACOL, IAROW, ICTXT,
     $                   ISEED, LDA, MB, MYCOL, MYROW, N,
     $                   NB, NPCOL, NPROW, BWL, BWU
*     ..
*     .. Array Arguments ..
      COMPLEX*16         A( LDA, * )
*     ..
*
*  Purpose
*  =======
*
*  PZBMATGEN : Parallel Complex Double precision Band MATrix GENerator.
*  (Re)Generate a distributed Band matrix A (or sub-matrix of A).
*
*  Arguments
*  =========
*
*  ICTXT   (global input) INTEGER
*          The BLACS context handle, indicating the global context of
*          the operation. The context itself is global.
*
*  AFORM   (global input) CHARACTER*1
*          if AFORM = 'L' : A is returned as a hermitian lower
*            triangular matrix, and is diagonally dominant.
*          if AFORM = 'U' : A is returned as a hermitian upper
*            triangular matrix, and is diagonally dominant.
*          if AFORM = 'G' : A is returned as a general matrix.
*          if AFORM = 'T' : A is returned as a general matrix in
*            tridiagonal-compatible form.
*
*  AFORM2  (global input) CHARACTER*1
*          if the matrix is general:
*            if AFORM2 = 'D' : A is returned diagonally dominant.
*            if AFORM2 != 'D' : A is not returned diagonally dominant.
*          if the matrix is symmetric or hermitian:
*            if AFORM2 = 'T' : A is returned in tridiagonally-compatible
*              form (a transpose form).
*            if AFORM2 != 'T' : A is returned in banded-compatible form.
*
*  M       (global input) INTEGER
*          The number of nonzero rows in the generated distributed
*           band matrix.
*
*  N       (global input) INTEGER
*          The number of columns in the generated distributed
*          matrix.
*
*  MB      (global input) INTEGER
*          The row blocking factor of the distributed matrix A.
*
*  NB      (global input) INTEGER
*          The column blocking factor of the distributed matrix A.
*
*  A       (local output) COMPLEX*16, pointer into the local memory
*          to an array of dimension ( LDA, * ) containing the local
*          pieces of the distributed matrix.
*
*  LDA     (local input) INTEGER
*          The leading dimension of the array containing the local
*          pieces of the distributed matrix A.
*
*  IAROW   (global input) INTEGER
*          The row processor coordinate which holds the first block
*          of the distributed matrix A.
*            A( DIAG_INDEX, I ) = A( DIAG_INDEX, I ) + BWL+BWU
*
*  IACOL   (global input) INTEGER
*          The column processor coordinate which holds the first
*          block of the distributed matrix A.
*
*  ISEED   (global input) INTEGER
*          The seed number to generate the distributed matrix A.
*
*  MYROW   (local input) INTEGER
*          The row process coordinate of the calling process.
*
*  MYCOL   (local input) INTEGER
*          The column process coordinate of the calling process.
*
*  NPROW   (global input) INTEGER
*          The number of process rows in the grid.
*
*  NPCOL   (global input) INTEGER
*          The number of process columns in the grid.
*
*  Notes
*  =====
*
*  This code is a simple wrapper around PZMATGEN, for band matrices.
*
*  =====================================================================
*
*  Code Developer: Andrew J. Cleary, University of Tennessee.
*    Current address: Lawrence Livermore National Labs.
*  This version released: August, 2001.
*
*  =====================================================================
*
*     ..
*     .. Parameters ..
      DOUBLE PRECISION   ONE, ZERO
      PARAMETER          ( ONE = 1.0D+0 )
      PARAMETER          ( ZERO = 0.0D+0 )
      COMPLEX*16         CONE, CZERO
      PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
      PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER           DIAG_INDEX, I, J, M_MATGEN, NQ, N_MATGEN,
     $                  START_INDEX
*     ..
*     .. External Subroutines ..
      EXTERNAL           PZMATGEN
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ICEIL, NUMROC
      EXTERNAL           ICEIL, NUMROC, LSAME
*     ..
*     .. Executable Statements ..
*
*
      IF( LSAME( AFORM, 'L' ).OR.LSAME( AFORM, 'U' ) ) THEN
         M_MATGEN = BWL + 1
         N_MATGEN = N
         START_INDEX = 1
         IF( LSAME( AFORM, 'L' ) ) THEN
            DIAG_INDEX = 1
         ELSE
            DIAG_INDEX = BWL + 1
         ENDIF
      ELSE
         M_MATGEN = BWL + BWU + 1
         N_MATGEN = N
         DIAG_INDEX = BWU + 1
         START_INDEX = 1
      ENDIF
*
      NQ = NUMROC( N, NB, MYCOL, IACOL, NPCOL )
*
*
*     Generate a random matrix initially
*
      IF( LSAME( AFORM, 'T' ) .OR.
     $  ( LSAME( AFORM2, 'T' ) ) ) THEN
*
          CALL PZMATGEN( ICTXT, 'T', 'N',
     $                        N_MATGEN, M_MATGEN,
     $                        NB, M_MATGEN, A( START_INDEX, 1 ),
     $                        LDA, IAROW, IACOL,
     $                        ISEED, 0, NQ, 0, M_MATGEN,
     $                        MYCOL, MYROW, NPCOL, NPROW )
*
      ELSE
*
          CALL PZMATGEN( ICTXT, 'N', 'N',
     $                        M_MATGEN, N_MATGEN,
     $                        M_MATGEN, NB, A( START_INDEX, 1 ),
     $                        LDA, IAROW, IACOL,
     $                        ISEED, 0, M_MATGEN, 0, NQ,
     $                        MYROW, MYCOL, NPROW, NPCOL )
*
*        Zero out padding at tops of columns
*
         DO 1000 J=1,NB
*
            DO 2000 I=1, LDA-M_MATGEN
*
*              Indexing goes negative; BMATGEN assumes that space
*              has been preallocated above the first column as it
*              has to be if the matrix is to be input to
*              Scalapack's band solvers.
*
               A( I-LDA+M_MATGEN, J ) = CZERO
*
 2000       CONTINUE
*
 1000    CONTINUE
*
      ENDIF
*
      IF( LSAME( AFORM2, 'D' ).OR.
     $  ( LSAME( AFORM, 'L' ).OR.LSAME( AFORM, 'U' ) ) ) THEN
*
*       Loop over diagonal elements stored on this processor.
*
*
       DO 330 I=1, NQ
         IF( LSAME( AFORM, 'T' ) .OR.
     $     ( LSAME( AFORM2, 'T' ) ) ) THEN
             IF( NPROW .EQ. 1 ) THEN
                A( I, DIAG_INDEX ) = DCMPLX( DBLE( A( I, DIAG_INDEX ) )
     $                               + DBLE( 2*( BWL+BWU+1 ) ) )
             ENDIF
          ELSE
             IF( NPROW .EQ. 1 ) THEN
                A( DIAG_INDEX, I ) = DCMPLX( DBLE( A( DIAG_INDEX, I ) )
     $                               + DBLE( 2*( BWL+BWU+1 ) ) )
             ENDIF
          END IF
  330  CONTINUE
*
*
      ELSE
*
*       Must add elements to keep condition of matrix in check
*
        DO 380 I=1, NQ
*
          IF( NPROW .EQ. 1 ) THEN
*
            IF( MOD(I+MYCOL*NB,2) .EQ. 1 ) THEN
                A( DIAG_INDEX+1, I ) =
     $                         DCMPLX( DBLE( A( DIAG_INDEX+1, I ) )
     $                         + DBLE( 2*( BWL+BWU+1 ) ) )
*
            ELSE
*
                A( DIAG_INDEX-1, I ) =
     $                          DCMPLX( DBLE( A( DIAG_INDEX-1, I ) )
     $                          + DBLE( 2*( BWL+BWU+1 ) ) )
            ENDIF
*
          ENDIF
*
  380   CONTINUE
*
      END IF
*
      RETURN
*
*     End of PZBMATGEN
*
      END

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack


text:
*DECK DHEFA
      SUBROUTINE DHEFA (A, LDA, N, IPVT, INFO, JOB)
      INTEGER LDA, N, IPVT(*), INFO, JOB
      DOUBLE PRECISION A(LDA,*)
C-----------------------------------------------------------------------
C     This routine is a modification of the LINPACK routine DGEFA and
C     performs an LU decomposition of an upper Hessenberg matrix A.
C     There are two options available:
C
C          (1)  performing a fresh factorization
C          (2)  updating the LU factors by adding a row and a
C               column to the matrix A.
C-----------------------------------------------------------------------
C     DHEFA factors an upper Hessenberg matrix by elimination.
C
C     On entry
C
C        A       DOUBLE PRECISION(LDA, N)
C                the matrix to be factored.
C
C        LDA     INTEGER
C                the leading dimension of the array  A .
C
C        N       INTEGER
C                the order of the matrix  A .
C
C        JOB     INTEGER
C                JOB = 1    means that a fresh factorization of the
C                           matrix A is desired.
C                JOB .ge. 2 means that the current factorization of A
C                           will be updated by the addition of a row
C                           and a column.
C
C     On return
C
C        A       an upper triangular matrix and the multipliers
C                which were used to obtain it.
C                The factorization can be written  A = L*U  where
C                L  is a product of permutation and unit lower
C                triangular matrices and  U  is upper triangular.
C
C        IPVT    INTEGER(N)
C                an integer vector of pivot indices.
C
C        INFO    INTEGER
C                = 0  normal value.
C                = k  if  U(k,k) .eq. 0.0 .  This is not an error
C                     condition for this subroutine, but it does
C                     indicate that DHESL will divide by zero if called.
C
C     Modification of LINPACK, by Peter Brown, LLNL.
C     Written 7/20/83.  This version dated 6/20/01.
C
C     BLAS called: DAXPY, IDAMAX
C-----------------------------------------------------------------------
      INTEGER IDAMAX, J, K, KM1, KP1, L, NM1
      DOUBLE PRECISION T
C
      IF (JOB .GT. 1) GO TO 80
C
C A new facorization is desired.  This is essentially the LINPACK
C code with the exception that we know there is only one nonzero
C element below the main diagonal.
C
C     Gaussian elimination with partial pivoting
C
      INFO = 0
      NM1 = N - 1
      IF (NM1 .LT. 1) GO TO 70
      DO 60 K = 1, NM1
         KP1 = K + 1
C
C        Find L = pivot index
C
         L = IDAMAX (2, A(K,K), 1) + K - 1
         IPVT(K) = L
C
C        Zero pivot implies this column already triangularized
C
         IF (A(L,K) .EQ. 0.0D0) GO TO 40
C
C           Interchange if necessary
C
            IF (L .EQ. K) GO TO 10
               T = A(L,K)
               A(L,K) = A(K,K)
               A(K,K) = T
   10       CONTINUE
C
C           Compute multipliers
C
            T = -1.0D0/A(K,K)
            A(K+1,K) = A(K+1,K)*T
C
C           Row elimination with column indexing
C
            DO 30 J = KP1, N
               T = A(L,J)
               IF (L .EQ. K) GO TO 20
                  A(L,J) = A(K,J)
                  A(K,J) = T
   20          CONTINUE
               CALL DAXPY (N-K, T, A(K+1,K), 1, A(K+1,J), 1)
   30       CONTINUE
         GO TO 50
   40    CONTINUE
            INFO = K
   50    CONTINUE
   60 CONTINUE
   70 CONTINUE
      IPVT(N) = N
      IF (A(N,N) .EQ. 0.0D0) INFO = N
      RETURN
C
C The old factorization of A will be updated.  A row and a column
C has been added to the matrix A.
C N-1 is now the old order of the matrix.
C
  80  CONTINUE
      NM1 = N - 1
C
C Perform row interchanges on the elements of the new column, and
C perform elimination operations on the elements using the multipliers.
C
      IF (NM1 .LE. 1) GO TO 105
      DO 100 K = 2,NM1
        KM1 = K - 1
        L = IPVT(KM1)
        T = A(L,N)
        IF (L .EQ. KM1) GO TO 90
          A(L,N) = A(KM1,N)
          A(KM1,N) = T
  90    CONTINUE
        A(K,N) = A(K,N) + A(K,KM1)*T
 100    CONTINUE
 105  CONTINUE
C
C Complete update of factorization by decomposing last 2x2 block.
C
      INFO = 0
C
C        Find L = pivot index
C
         L = IDAMAX (2, A(NM1,NM1), 1) + NM1 - 1
         IPVT(NM1) = L
C
C        Zero pivot implies this column already triangularized
C
         IF (A(L,NM1) .EQ. 0.0D0) GO TO 140
C
C           Interchange if necessary
C
            IF (L .EQ. NM1) GO TO 110
               T = A(L,NM1)
               A(L,NM1) = A(NM1,NM1)
               A(NM1,NM1) = T
  110       CONTINUE
C
C           Compute multipliers
C
            T = -1.0D0/A(NM1,NM1)
            A(N,NM1) = A(N,NM1)*T
C
C           Row elimination with column indexing
C
               T = A(L,N)
               IF (L .EQ. NM1) GO TO 120
                  A(L,N) = A(NM1,N)
                  A(NM1,N) = T
  120          CONTINUE
               A(N,N) = A(N,N) + T*A(N,NM1)
         GO TO 150
  140    CONTINUE
            INFO = NM1
  150    CONTINUE
      IPVT(N) = N
      IF (A(N,N) .EQ. 0.0D0) INFO = N
      RETURN
C----------------------- End of Subroutine DHEFA -----------------------
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


SUBROUTINE read_inflow_sct_step(inflow_file)

!-------------------------------------------------------------
!
! Reads fixed inflow BC
!
!-------------------------------------------------------------
USE invar
IMPLICIT NONE

CHARACTER(30), INTENT(IN) :: inflow_file
INTEGER :: octant,n
INTEGER :: ix,iy,iz,jx,jy,jz,dir
INTEGER :: sgn_mu,sgn_eta,sgn_xi

ALLOCATE( frbc(nx,nz,apo,4,1,1) )
ALLOCATE( babc(nx,nz,apo,4,1,1) )
ALLOCATE( lebc(ny,nz,apo,4,1,1) )
ALLOCATE( ribc(ny,nz,apo,4,1,1) )
ALLOCATE( bobc(nx,ny,apo,4,1,1) )
ALLOCATE( tobc(nx,ny,apo,4,1,1) )

OPEN(UNIT=12, FILE=inflow_file,STATUS = "OLD", ACTION = "READ",FORM='UNFORMATTED')

DO octant=1,8
   DO n=1,apo
      READ(12) sgn_mu,sgn_eta,sgn_xi
      ! >> front/back
      IF(sgn_xi>0 .and. sgn_mu>0) THEN
         dir=1
      ELSE IF (sgn_xi>0 .and. sgn_mu<0) THEN
         dir=2
      ELSE IF (sgn_xi<0 .and. sgn_mu>0) THEN
         dir=3
      ELSE
         dir=4
      END IF
      !
      IF (sgn_eta>0) THEN
         DO iz=1,nz
            DO ix=1,nx
               READ(12) frbc(ix,iz,n,dir,1,1)
            END DO
         END DO
      ELSE
         DO iz=1,nz
            DO ix=1,nx
               READ(12) babc(ix,iz,n,dir,1,1)
            END DO
         END DO
      END IF
      ! >> left/right
      IF(sgn_eta>0 .and. sgn_xi>0) THEN
         dir=1
      ELSE IF (sgn_eta>0 .and. sgn_xi<0) THEN
         dir=2
      ELSE IF (sgn_eta<0 .and. sgn_xi>0) THEN
         dir=3
      ELSE
         dir=4
      END IF
      !
      IF (sgn_mu>0) THEN
        DO iy=1,ny
          DO iz=1,nz
            READ(12) lebc(iy,iz,n,dir,1,1)
          END DO
        END DO
      ELSE
        DO iy=1,ny
          DO iz=1,nz
             READ(12) ribc(iy,iz,n,dir,1,1)
          END DO
        END DO
      END IF
      ! >> bottom/top
      IF(sgn_mu>0 .and. sgn_eta>0) THEN
         dir=1
      ELSE IF (sgn_mu>0 .and. sgn_eta<0) THEN
         dir=2
      ELSE IF (sgn_mu<0 .and. sgn_eta>0) THEN
         dir=3
      ELSE
         dir=4
      END IF
      !
      IF (sgn_xi>0) THEN
        DO ix=1,nx
          DO iy=1,ny
            READ(12) bobc(ix,iy,n,dir,1,1)
          END DO
        END DO
      ELSE
        DO ix=1,nx
          DO iy=1,ny
            READ(12) tobc(ix,iy,n,dir,1,1)
          END DO
        END DO
      END IF
      !
   END DO
END DO

!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Needs to be removed
! lebc=1.0d0
! ribc=1.0d0
! frbc=0.5d0
! babc=0.5d0
! bobc=0.25d0
! tobc=0.25d0
! Needs to be removed
!!!!!!!!!!!!!!!!!!!!!!!!!!!

CLOSE(UNIT=12)

END SUBROUTINE read_inflow_sct_step