text:
From Spring 2001 to Winter 2002 Tiqula Bledsoe was the ASUCD President and ran as a member of the LEAD Slate slate. Before he was president, he was an ASUCD Senate ASUCD senator.

The election that lead to him taking office was even more dramatic than usual for the time. Opponents got his slate disqualified for violating election codes in their campaign, such as spending too much money and posting large fliers. However, after playing the race card (including Bledsoe or a supporter posting fake KKK fliers on campus), the ASUCD Court reversed their decision amidst massive pressure from various ethnic minority groups on campus.

He called himself the Black Caesar of ASUCD.  This, of course, led to some controversy, including some of the hate that the LEAD slate received. Also, this reference caused problems for the California Aggie when they referenced him in their endoftheyear The California Aggie: Spoof edition spoof edition. During his presidential term Bledsoe revoked the authority of the ASUCD Court in 2001.  Later in 2005 the California Aggie featured an http://media.www.californiaaggie.com/media/storage/paper981/news/2005/01/20/Opinion/Editorial.Past.Elections.Problems.Addressed.In.Bills1319260.shtml editorial in support of Senate Bill #6 to revive the ASUCD Court Courts powers. In 2004, Bledsoe http://media.www.californiaaggie.com/media/storage/paper981/news/2004/10/05/Opinion/Uc.More.Segregated.Than.University.Of.Mississippi1318232.shtm wrote a letter to the California Aggie stating that the UC System was more segregated than the http://www.olemiss.edu/ University of Mississippi. Bledsoe was then http://media.www.californiaaggie.com/media/storage/paper981/news/2004/10/07/Opinion/Letters.To.The.Editor1318279.shtml commented on concerning the subject in a letter to the California Aggie editor.

Tiqula Bledsoe however http://www.yolocounty.org/org/bos/agendas/2002/032602/06.pdf received praise from the Yolo County supervisors and State Assemblywoman Lois Wolk in a resolution signed on March 16, 2002.  As ASUCD President, Bledsoe http://www.democraticunderground.com/cgibin/duforum/duboard.cgi?azshow_thread&om1349&forumDCForumID43&archiveyes opposed the Pledge Bill of 2001.

As an ordained minister in the Church of God and Christ, Bledsoe was also routinely called Reverend Tiqula Bledsoe, and jokingly referred to himself by numerous additional ego maniacal nicknames, such as the Mohamed Ali of ASUCD, and even half man, half amazing.

He left the office of ASUCD President with an unusually high debt in his allocated budget, causing ASUCD to have to do some juggling to survive the rest of the fiscal year. It also came to light after he left office that he had attempted to use ASUCD funds for his personal limo to the Presidents Ball (putting on a ball isnt exactly a normal ASUCD activity either, and in this case it ran a huge debt whereas it was supposed to break even), and to pay for parking tickets hed incurred in the state capitol, among other things. Mark Champagne prevented this money from actually being spent  fortunately for all students who are not Tiqula Bledsoe.

After graduating from UC Davis, Tiqula Bledsoe became the Sustainable Development and Fair Growth Fellow at wiki:wikipedia:Greenlining_Institute The Greenlining Institute. Afterwards, he became the Executive Manager of the California State Conference of the http://www.naacp.org/ NAACP. He also went on to earn an Executive wiki:MBA from California State University, Sacramento. In 2005, he became the Executive Vice President of http://www.deebabys.com/ Dee Babys Realty.

His given name is pronounced taiKYUla, if anyone was wondering.


program main
  implicit none
  real(8),parameter :: pi = 4d0*atan(1d0)
  complex(8),parameter :: zI = (0d0,1d0)
  integer,parameter :: Nt = 6000, Nw = 500
  real(8),parameter :: wi = 0.1d0/27.211d0, wf = 20d0/27.211d0, dw = (wf-wi)/Nw
  real(8) :: tt(0:Nt),Act(0:Nt),jt(0:Nt)
  real(8) :: dt,ww,xx,ff
  integer :: it,iw
  complex(8) :: zAw,zEw,zJw,zsigma,zeps
  real(8) :: kick_impulsive

  kick_impulsive = 1d-4

  do it = 0,Nt
    read(*,*)tt(it),Act(it),jt(it)
  end do
  dt = tt(1) - tt(0)
  jt(0) = jt(0)*0.5d0

  open(21,file="zeps.out")
  do iw = 0,Nw
    ww = wi + dw*iw
    zAw = 0d0; zJw = 0d0
    do it = 0,Nt
      xx = (tt(it)-tt(0))/(tt(Nt) - tt(0))
      ff = 1d0 -3d0 * xx**2 + 2d0*xx**3
      zAw = zAw + exp(zI*ww*tt(it)) *Act(it)
      zjw = zjw + exp(zI*ww*tt(it)) *jt(it)*ff
    end do
    zEw = -(-zI*ww)*zAw*dt
    zJw = - zJw*dt
!    zsigma = zJw/zEw
    zsigma = zJw/kick_impulsive
    zeps = 1d0 + 4d0*pi*zi*zsigma/ww

    write(21,"(999e26.16e3)")ww,zeps,zsigma

  end do

  close(21)

end program main


text:
Carlos Canales is a licensed clinical psychologists who enjoys therapy at many levels. In his words, I am attentive to the complexity of treatment, attuned to the multiple layers of experience that occur when two people interact. I enjoy making connections with clients and, together, exploring what brought them to counseling, how their history informs their symptoms or present situation, and how needs are communicated. Basically, how life is experienced and observed. I also value the encounter created in the present moment by two humans invested in making changes towards health.

Services

  Individual Therapy
  Couples Therapy
  Group Therapy

All services can be provided in Spanish.

See also: Counselors and Therapists

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
!! Copyright (c) 1997 Neil N. Carlson
!!
!! This file is part of MFE1 which is released under the MIT license.  See the
!! file LICENSE or visit http://opensource.org/licenses/MIT for details.
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module norm_data

  use mfe_constants, only: wp
  use mfe_types, only: NodeVar
  implicit none
  private

  type(NodeVar), save, public :: ptol
  real(kind=wp), save, public :: rtol

end module norm_data

module norm_procs

  use mfe_constants, only: wp
  use mfe_types, only: NodeVar
  use mfe_data, only: dxmin
  use norm_data
  use common_io, only: element_info
  implicit none
  private

  public :: eval_norm, check_soln

  real(kind=wp), dimension(:), allocatable, save, public :: dx

  contains

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!
   !!  EVAL_NORM
   !!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    function eval_norm (du, key) result (norm)

      type(NodeVar), dimension(:), intent(in) :: du
      integer, intent(in) :: key
      real(kind=wp) :: norm

      real(kind=wp), dimension(size(du)-1) :: del_dx

      ! Weighted max-norm.
      norm = max(maxval(abs(du % x)) / ptol % x, maxval(abs(du % u)) / ptol % u)

      del_dx = du(2:size(du)) % x - du(1:size(du)-1) % x

      ! Relative norm on element lengths.
      norm = max (norm, maxval(abs(del_dx) / dx) / rtol)

    end function eval_norm

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!
   !!  CHECK_SOLN
   !!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    subroutine check_soln (u, key, errc)

      type(NodeVar), dimension(:), intent(in) :: u
      integer, intent(in)  :: key
      integer, intent(out) :: errc

      integer, dimension(1) :: loc

      if (.not. allocated(dx)) then
        allocate (dx(size(u)-1))
      end if

      dx = u(2:size(u)) % x - u(1:size(u)-1) % x
      loc = minloc(dx)

      if (dx(loc(1)) < dxmin) then
        call element_info ("CHECK_SOLN: BAD ELEMENT", loc(1))
        errc = 1
      else
        errc = 0
      end if

    end subroutine check_soln

end module norm_procs


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

!	OpenMP Parallel Region
!	parallel private subroutine call

	program p
	parameter(n=10)
	integer result(n)
	integer expect(n)
	data expect/101,201,301,401,102,202,302,402,103,203/
	call sp1(result,n)
	!print *,result
	call check(result,expect,n)
	end

	subroutine sp1(x,n)
	 integer n
	 integer x(n)
	 integer omp_get_thread_num
	 integer omp_get_num_threads
!$omp   parallel private(iam,np,ipoints)
	 iam = omp_get_thread_num()+1
	 np = omp_get_num_threads()
	 call subdomain(x,iam,n,np)
!$omp   end parallel
	end

	subroutine subdomain(x,iam,n,np)
	integer n
	integer x(n),iam,np
	integer i,j
	j = 0
	do i = iam,n,np
	 j = j + 1
	 x(i) = iam*100 + j
	enddo
	end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
subroutine ale_reaMPI(a)
   use typre
   use Mod_BroadCastBuffer
   use Mod_ALEmov
   use Mod_sldAlemov
   implicit none
   class(AlemovProblem) :: a
   type(BroadCastBuffer) :: BBuffer

   !Initialize Buffer
   call BBuffer%SetMPI(a%MPIcomm,a%MPIsize,a%MPIroot,a%MPIrank)
   call BBuffer%SetLOutputFile(a%lun_memo,a%lun_outpu)
   call BBuffer%Initialize(100,500)

   !Communicate nsi_reaphy
   call BBuffer%Add(a%kfl_isFixedMesh)
   call BBuffer%Add(a%kfl_RemeshingCriteria)
   call BBuffer%Add(a%kfl_doAitken)

   !If we are working with solid ALE and want to print strains
   select type (a)
   type is (sldAlemovProblem)
       call BBuffer%Add(a%kfl_printPrincipalStresses)
       call BBuffer%Add(a%sldale_dispRelax)
       call BBuffer%Add(a%sldale_dispRelax_max)
       call a%solid%SetMPI(a%MPIcomm,a%MPIsize,a%MPIroot,a%MPIrank)
       call a%solid%SpecificReaMPI
   end select

   !BroadCast and DeallocateBuffer
   call BBuffer%BroadCast
   call BBuffer%Dealloc

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
	FUNCTION ZPADC(F)
C
C	THIS FUNCTION RETURNS THE GAIN OF THE DC CHANNEL OF THE WIND Z PREAMP
C	THE FUNCTION ITSELF RETURNS THE MAGNITUDE OF THE GAIN, AND THE PHASE
C	AND THE COMPLEX GAIN ARE RETURNED IN A COMMON BLOCK
C
C	TO BE CHANGED FOR Y,Z:  R11,R32,C27,R9,C9,R10
C	DONE AND TESTED
C
	COMMON /GAINBLK/ PHASE,CGAIN        		  ! PHASE IS IN RADIANS
C	COMMON /COMPBLK/ CVAR(5)
C
	COMPLEX CGAIN,Y1,Z1,Y2,Z2
	DATA R9,C9,R10,CINP /18.E6,84.E-12,2.E8,4.99E-12/
	DATA R33,C16 /2.2E+03,1.202E-9/
	DATA RCOMP,CCOMP /22.1E3,47.E-12/
	DATA R29 /22.1E3/
	DATA TWOPI /6.2831853/
C
C	C16 = 1.203E-9                             ! ADJUSTING EX
C	C9  = 23.6E-12                            ! ADJUSTING EX
C	CINP = 4.99E-12                          ! ADJUSTING EZ
C
C	C16 = CVAR(1)
C	C9 = CVAR(2)
C	CINP = CVAR(3)
C
	W = TWOPI*F
	CGAIN = 1.
C	DIVISION OF SIGNAL BEFORE FOLLOWER
	Y1 = CMPLX(1./R9,W*C9)
	Y2 = CMPLX(1./R10,W*CINP)
	Z1 = 1./Y1
	Z2 = 1./Y2
	CGAIN = CGAIN*Z2/(Z1+Z2)
C	EFFECT OF OUTPUT RC
	WC16 = 1./(W*C16)
C	PRINT*,'W,Z1,WC16',W,Z1,WC16
C	PRINT*,'F,CGAIN,PHASE
	CGAIN = CGAIN*CMPLX(0.,-WC16)/CMPLX(R33,-WC16)
C
C	EFFECT OF INPUT (211) BOARD
C
	Y1 = 1./R29 + CMPLX(0.,W*C16)
	CGAIN = CGAIN/(1. + R33*Y1)
C
	GNIP = AIMAG(CGAIN)
	GNRP = CGAIN
	PHASE = ATAN2(GNIP,GNRP)
	ZPADC = CABS(CGAIN)
	RETURN
	END

ext: for


text:
module parff_parsed_value_m
    implicit none
    private
    public :: parsed_value_t

    type, abstract :: parsed_value_t
    end type
end module

max_stars_repo_name: everythingfunctional/parff

max_issues_repo_name: everythingfunctional/parff

max_forks_repo_name: everythingfunctional/parff


text:
!      ------AS ABAQUS ARRAY SECTION-------------------
       function asabqarray_2s(T,i,j)
        implicit none

        type(Tensor2s), intent(in) :: T
        integer, intent(in) :: i
        integer, intent(in), optional :: j
        real(kind=8), dimension(i) :: asabqarray_2s
        integer :: i1

        asabqarray_2s(1:min(i,4)) = T%a6(1:min(i,4))

        if (i.ge.5) then
         asabqarray_2s(5) = T%a6(6)
        end if

        if (i.ge.6) then
         asabqarray_2s(6) = T%a6(5)
        endif

       end function asabqarray_2s

       function asabqarray_4s(T,i,j)
        implicit none

        type(Tensor4s), intent(in) :: T
        integer, intent(in) :: i,j
        real(kind=8), dimension(i,j) :: asabqarray_4s
        integer, dimension(6) :: imap
        integer :: i1,j1

        imap = (/1,2,3,4,6,5/)

        forall (i1=1:i,j1=1:j) asabqarray_4s(i1,j1) =
     *            T%a6b6(imap(i1),imap(j1))

       end function asabqarray_4s

max_stars_repo_name: kengwit/ttb

max_issues_repo_name: kengwit/ttb

max_forks_repo_name: kengwit/ttb


text:
      SUBROUTINE TGLOC04(NTARGS,TARGNAME,NPARMS,PARMS,
     *      KSYSNAT,KCENNAT,TARGUNIT,LUERR,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C THIS ROUTINE IS PART OF THE TOSS TARGET LOCATION PACKAGE, TGLOC.  IT
C COMPUTES UNIT VECTORS IN LOCAL ORBITAL COORDINATES FOR TARGETS
C DEFINED BY AZIMUTH AND ELEVATION IN THOSE COORDINATES.
C
C VAR     DIM    TYPE I/O DESCRIPTION
C ---     ---    ---- --- -----------
C
C NTARGS   1      I*4  I  NUMBER OF TARGETS TO BE PROCESSED.
C                         IF ZERO OR NEGATIVE, NONE ARE PROCESSED.
C
C TARGNAME NTARGS CH*16 I TARGET NAMES;  USED FOR ERROR MESSAGES,
C                         OTHERWISE THE NAME HAS NO SIGNIFICANCE.
C
C NPARMS   1      I*4  I  THE NUMBER OF PARAMETERS SUPPLIED IN THE PARMS
C                         ARRAY.  MUST BE 2 OR GREATER.
C
C PARMS           R*8  I  SAME AS TGLOC'S DESCRIPTION OF PARMS FOR
C    NPARMS,NTARGS        TARGET TYPE 4.
C
C KSYSNAT  1      I*4  O  IDENTIFIES THE NATURAL COORDINATE SYSTEM'S
C                         ORIENTATION.  IS SET TO 12 TO INDICATE THE
C                         ORIENTATION IS LOCAL ORBITAL.
C
C KCENNAT  1      I*4  O  IDENTIFIES THE NATURAL COORDINATE SYSTEM'S
C                         ORIGIN.  IS SET TO 9999 TO INDICATE THE ORIGIN
C                         IS IRRELEVANT.
C
C TARGUNIT 3,NTARGS R*8 O TARGUNIT(-,I) IS THE UNIT VECTOR OF THE I'TH
C                         TARGET.  IN LOCAL ORBITAL COORDINATES.
C
C LUERR    1      I*4   I FORTRAN UNIT NUMBER FOR ERROR MESSAGES.
C                          = 0/NEGATIVE, NO MESSAGES POSSIBLE.
C
C IERR     1      I*4   O ERROR RETURN FLAG.
C                          = 0, NO ERROR.
C                          = OTHERWISE, ERROR.
C
C***********************************************************************
C
C BY C PETRUZZO, GSFC/742, 8/85.
C        MODIFIED....
C
C***********************************************************************
C
      REAL*8 PARMS(NPARMS,1),TARGUNIT(3,1)
      CHARACTER*16 TARGNAME(1)
C
      IERR = 0
      KSYSNAT = 12     ! SYSTEM ORIENTATION IS LOCAL ORBITAL
      KCENNAT = 9999   ! INDEPENDENT OF ORIGIN
      IF(NTARGS.LE.0) GO TO 9999
C
C GET TARGET DIRECTION VECTOR IN 'NATURAL' SYSTEM.
C
      DO ITARG=1,NTARGS
        CALL TGLOC04A(NPARMS,PARMS(1,ITARG),TARGNAME(ITARG),
     *           LUERR,IERR1)
        IF(IERR1.EQ.0) THEN
          AZ = PARMS(1,ITARG)
          EL = PARMS(2,ITARG)
          CALL XYZSPH(1,AZ,EL,1.D0,TARGUNIT(1,ITARG),1.D0)
        ELSE
          IERR = 1
          TARGUNIT(1,ITARG) = 'ERROR'
          TARGUNIT(2,ITARG) = 'ERROR'
          TARGUNIT(3,ITARG) = 'ERROR'
          END IF
        END DO
C
 9999 CONTINUE
      RETURN
      END


text:
C ----------------------------------------------------------------- C
C - This is a parton level only analysis for ttbar production     - C
C - tops bs Ws etc are constructed using MC truth only - no jets  - C
C - The MC truth reconstruction has been tested (see sanity check - C
C - code in the analysis below).                                  - C
C - Since it is parton level you need to comment out from         - C
C - CALL HWDHOB down to CALL HWDHOB inclusive in main-HERWIG.f .  - C
C - Also at some point, the showering went into what looked like  - C
C - an infinite loop after 137K events - gdb said it was in       - C
C - HWHGUP. The same glitch did not occur with *** herwig6520.f *** C
C - I also eliminated the analysis as a possible cause (it occurs - C
C - with HWANAL removed). I did not see anything fishy with the   - C
C - Tevatron, semileptonic event that got caught.                 - C
C ----------------------------------------------------------------- C

c  The next subroutines, open some histograms and prepare them
c      to receive data
c  You can substitute these  with your favourite ones
c  init   :  opens the histograms
c  topout :  closes them
c  pwhgfill  :  fills the histograms with data

      subroutine init_hist_KN
      implicit none
      include  'LesHouches.h'
      include 'pwhg_math.h'
      integer lastnbx,j,l
      character * 20 prefix
      integer nbins
      parameter (nbins=11)
      real * 8 pT_tt_bins(nbins+1)
      data pT_tt_bins/  0d0, 10d0, 25d0, 50d0,100d0,
     1     150d0,200d0,250d0,300d0,400d0,600d0, 900d0/
      real * 8 m_tt_bins(nbins+1)
      data m_tt_bins/ 320d0,360d0,400d0,450d0,500d0,
     1     550d0,600d0,650d0,700d0,800d0,900d0,1000d0/
      external lastnbx

      call inihists

      do j=1,18
         if(j.eq.1) then
            prefix='t'
         elseif(j.eq.2) then
            prefix='tb'
         elseif(j.eq.3) then
            prefix='t-pt5'
         elseif(j.eq.4) then
            prefix='tb-pt5'
         elseif(j.eq.5) then
            prefix='t-pt20'
         elseif(j.eq.6) then
            prefix='tb-pt20'
         elseif(j.eq.7) then
            prefix='btop'
         elseif(j.eq.8) then
            prefix='btop-pt5'
         elseif(j.eq.9) then
            prefix='btop-pt20'
         elseif(j.eq.10) then
            prefix='bbtop'
         elseif(j.eq.11) then
            prefix='bbtop-pt5'
         elseif(j.eq.12) then
            prefix='bbtop-pt20'
         elseif(j.eq.13) then
            prefix='lwp'
         elseif(j.eq.14) then
            prefix='lwp-pt5'
         elseif(j.eq.15) then
            prefix='lwp-pt20'
         elseif(j.eq.16) then
            prefix='lwm'
         elseif(j.eq.17) then
            prefix='lwm-pt5'
         elseif(j.eq.18) then
            prefix='lwm-pt20'
         endif
         l=lastnbx(prefix)
         call bookupeqbins(prefix(1:l)//'_y'  ,0.2d0,-4d0,4d0)
         call bookupeqbins(prefix(1:l)//'_eta',0.2d0,-4d0,4d0)
         call bookupeqbins(prefix(1:l)//'_pt' ,5d0,0d0,400d0)
         call bookupeqbins(prefix(1:l)//'_m'  ,2d0,-0.5d0,201.5d0)
      enddo

      end

      subroutine analysis_KN(dsig0)
      implicit none
      include 'hepevt.h'
      include 'pwhg_math.h'
      include 'LesHouches.h'
      character * 6 whcprg
      common/cwhcprg/whcprg
      integer jpref
      character * 20 prefix(18)
      common/ccccprefix/jpref,prefix
      real * 8  dsig0,dsig
      logical   ini
      data      ini/.true./
      save      ini
      integer   ihep                ! HEPEVT index.
      real * 8 p_top(4),p_tb(4),p_wp(4),p_wm(4),p_lwp(4),p_lwm(4),
     1         p_b(4),p_bb(4),y,eta,pt,mass,ptt
      integer j
      real * 8 prodvec2

      dsig  = dsig0

      if(whcprg.eq.'NLO'.or.whcprg.eq.'LHE') then
         p_top=phep(1:4,3)
         p_tb=phep(1:4,4)
         p_wp=phep(1:4,5)
         p_wm=phep(1:4,6)
         p_lwp=phep(1:4,8)
         p_lwm=phep(1:4,9)
         p_b=phep(1:4,11)
         p_bb=phep(1:4,12)
C --------------------------------------------- C
C - LHE PARTICLE TOP RECONSTRUCTION: MC TRUTH - C
C --------------------------------------------- C
      else
         write(*,*) ' Now only NLO cross sections are implemented'
      endif

      if(nhep.eq.13) then
         ptt=phep(1,13)**2+phep(2,13)**2
         do j=7,12
            ptt=min(ptt,prodvec2(phep(:,j),phep(:,13)))
         enddo
         ptt=sqrt(ptt)
      else
         ptt=0
      endif

      call yetaptmassplot(p_top,dsig,'t')
      if(ptt.gt.5) call yetaptmassplot(p_top,dsig,'t-pt5')
      if(ptt.gt.20) call yetaptmassplot(p_top,dsig,'t-pt20')

      call yetaptmassplot(p_tb,dsig,'tb')
      if(ptt.gt.5) call yetaptmassplot(p_tb,dsig,'tb-pt5')
      if(ptt.gt.20) call yetaptmassplot(p_tb,dsig,'tb-pt20')

      call yetaptmassplot(p_b,dsig,'btop')
      if(ptt.gt.5) call yetaptmassplot(p_b,dsig,'btop-pt5')
      if(ptt.gt.20) call yetaptmassplot(p_b,dsig,'btop-pt20')

      call yetaptmassplot(p_bb,dsig,'bbtop')
      if(ptt.gt.5) call yetaptmassplot(p_bb,dsig,'bbtop-pt5')
      if(ptt.gt.20) call yetaptmassplot(p_bb,dsig,'bbtop-pt20')

      call yetaptmassplot(p_lwp,dsig,'lwp')
      if(ptt.gt.5) call yetaptmassplot(p_lwp,dsig,'lwp-pt5')
      if(ptt.gt.20) call yetaptmassplot(p_lwp,dsig,'lwp-pt20')

      call yetaptmassplot(p_lwm,dsig,'lwm')
      if(ptt.gt.5) call yetaptmassplot(p_lwm,dsig,'lwm-pt5')
      if(ptt.gt.20) call yetaptmassplot(p_lwm,dsig,'lwm-pt20')

      end

C ----------------------------------------------- C
C - SANITY CHECK OF TRUTH SHOWER RECONSTRUCTION - C
C ----------------------------------------------- C
C - p_tmp is the total momentum in minus the total momentum out.
C - If the sum of the absolute values of p_tmp are greater than
C - epsilon (1 GeV) a warning plus debugging output is printed
C - on the screen (see below). At Tevatron with 300K events I
C - see in Herwig messages displayed for epsilon equal to 5 GeV
C - every 30K events - n.b. these violations can be directly
C - and totally attributed exclusively to differences between
C - the Wminus momentum reconstructed and it's true value (see
C - below), what's more 90% of them are associated to the W
C - decaying to strange+charm which then shower - the strange
C - and charm momenta are the things that change quite a bit
C - between W decay and post-showering. Pythia is totally fine
C - for 300K events - no messages. Possibly this is a HW issue then.

      subroutine shower_sanity_checks(p_top,p_tba,p_beam,p_Wp,p_Wm,
     1                                p_b  ,p_bba,p_jet)
      implicit none
      include 'hepevt.h'
      include 'LesHouches.h'
      character * 6 whcprg
      common/cwhcprg/whcprg
      real * 8  p_top(4),p_tba(4),p_beam(4),p_tmp(4)
      real * 8  p_Wp(4) ,p_Wm(4) ,p_b(4)   ,p_bba(4),p_jet(4)
c$$$      real * 8           y_top,eta_top,pt_top,m_top,
c$$$     1                   y_tba,eta_tba,pt_tba,m_tba
c$$$      real * 8           y_jet,eta_jet,pt_jet,m_jet ! Powheg _radiated_ jet
      real * 8  tmp
      integer   n_sanity_violations,n_event
      data      n_sanity_violations,n_event/0,0/
      save      n_sanity_violations,n_event
      integer   ixx

c$$$      call getyetaptmass(p_top,y_top,eta_top,pt_top,m_top)
c$$$      call getyetaptmass(p_tba,y_tba,eta_tba,pt_tba,m_tba)
c$$$      call getyetaptmass(p_jet,y_jet,eta_jet,pt_jet,m_jet)

      n_event=n_event+1

      p_tmp=p_top+p_tba+p_jet+p_beam-phep(1:4,1)-phep(1:4,2)
      tmp=0d0
      do ixx=1,4
         tmp=tmp+abs(p_tmp(ixx))
      enddo
      if(abs(tmp).gt.5d0) then
         n_sanity_violations=n_sanity_violations+1
         if(n_sanity_violations.le.10) then
            write(*,*) ''
            write(*,*) ''
            if(whcprg.eq.'HERWIG') then
               write(*,*) 'HERWIG analysis'
CC                  call HWUEPR
CC                  call HWUBPR
            else if(whcprg.eq.'PYTHIA') then
               write(*,*) 'PYTHIA analysis'
CC                  call pylist(5)
            else if(whcprg.eq.'LHE') then
               write(*,*) 'LHE analysis'
               do ixx=1,nhep
                  write(*,*) 'ihep =',ixx,' id = ',idhep(ixx),
     1                       ' mom. = ',phep(1:4,ixx)
               enddo
            endif
            write(*,*) ''
            write(*,*) 'Consistency err: momentum in != momentum out'
            write(*,*) 'Event number ',n_event
            write(*,*) 'Debugging output follows.'
            write(*,*) 'nup   = ',nup
            write(*,*) 'p_in  = ',phep(1:4,1)+phep(1:4,2)
            write(*,*) 'p_out = ',p_top+p_tba+p_jet+p_beam

            if(whcprg.eq.'PYTHIA') then
               p_tmp=phep(1:4,7)
               write(*,*) 'p_top true = ',p_tmp
            else if(whcprg.eq.'HERWIG') then
               do ixx=1,nhep
                  if(isthep(ixx).eq.155.and.idhep(ixx).eq.6) then
                     p_tmp=phep(1:4,ixx)
                  endif
               enddo
               write(*,*) 'p_top true = ',p_tmp
            endif
            write(*,*) 'p_top      = ',p_top
c$$$            write(*,*) 'm_top      = ',m_top

            do ixx=1,nhep
               if(whcprg.eq.'PYTHIA') then
                  if(isthep(ixx).eq.2.and.idhep(ixx).eq.24) then
                     p_tmp=phep(1:4,ixx)
                     write(*,*) 'p_Wp true  = ',p_tmp
                  endif
               else if(whcprg.eq.'HERWIG') then
                  if(isthep(ixx).eq.155.and.idhep(ixx).eq.24) then
                     p_tmp=phep(1:4,ixx)
                     write(*,*) 'p_Wp true  = ',p_tmp
                  endif
               endif
            enddo
            write(*,*) 'p_Wp       = ',p_Wp

            if(whcprg.eq.'PYTHIA') then
               p_tmp=phep(1:4,11)
               write(*,*) 'p_b true   = ',p_tmp
            else if(whcprg.eq.'HERWIG') then
               do ixx=1,nhep
                  if(isthep(ixx).eq.124.and.idhep(ixx).eq.5) then
                     p_tmp=phep(1:4,ixx)
                  endif
               enddo
               write(*,*) 'p_b true   = ',p_tmp
            endif
            write(*,*) 'p_b        = ',p_b

            if(whcprg.eq.'PYTHIA') then
               p_tmp=phep(1:4,8)
               write(*,*) 'p_tba true = ',p_tmp
            else if(whcprg.eq.'HERWIG') then
               do ixx=1,nhep
                  if(isthep(ixx).eq.155.and.idhep(ixx).eq.-6) then
                     p_tmp=phep(1:4,ixx)
                  endif
               enddo
               write(*,*) 'p_tba true = ',p_tmp
            endif
            write(*,*) 'p_tba      = ',p_tba
c$$$            write(*,*) 'm_tba      = ',m_tba

            do ixx=1,nhep
               if(whcprg.eq.'PYTHIA') then
                  if(isthep(ixx).eq.2.and.idhep(ixx).eq.-24) then
                     p_tmp=phep(1:4,ixx)
                     write(*,*) 'p_Wm true  = ',p_tmp
                  endif
               else if(whcprg.eq.'HERWIG') then
                  if(isthep(ixx).eq.155.and.idhep(ixx).eq.-24) then
                     p_tmp=phep(1:4,ixx)
                     write(*,*) 'p_Wm true  = ',p_tmp
                  endif
               endif
            enddo
            write(*,*) 'p_Wm       = ',p_Wm

            if(whcprg.eq.'PYTHIA') then
               p_tmp=phep(1:4,13)
               write(*,*) 'p_bba true = ',p_tmp
            else if(whcprg.eq.'HERWIG') then
               do ixx=1,nhep
                  if(isthep(ixx).eq.124.and.idhep(ixx).eq.-5) then
                     p_tmp=phep(1:4,ixx)
                  endif
               enddo
               write(*,*) 'p_bba true = ',p_tmp
            endif
            write(*,*) 'p_bba      = ',p_bba

            if(whcprg.eq.'PYTHIA') then
               write(*,*) 'p_jet true = ',phep(1:4,9)
            else if(whcprg.eq.'HERWIG') then
               write(*,*) 'p_jet true = ',phep(1:4,jdahep(1,9))
            endif
            write(*,*) 'p_jet      = ',p_jet

            write(*,*) 'p_beam      = ',p_beam

         else if(n_sanity_violations.eq.11) then
            write(*,*) 'Analysis found 10 momentum conservation'
            write(*,*) 'inconsistencies - no more debug ouput'
            write(*,*) 'will be given on these, just warnings.'
         else if(mod(n_sanity_violations,20).eq.0) then
            write(*,*) 'Analysis has now found',
     1           n_sanity_violations,' momentum inconsistencies.'
         endif
      endif

      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
*
* invert.f subroutine to invert matrix for barotropic streamfunction
* c grid version variable depth last change 20/10/94
* version for non-trivial coastlines started 10/4/00
* wind part removed to wind.f 9/2/1
* variable drag 7/5/2
* variable ds 6/12/4 nre
*
      subroutine invert
#include "ocean.cmn"

      integer i,j,k,l,n,m,im

      real tv,tv1,rat

      n = imax
      m = jmax + 1

      do 710 i=1,n*m
         do 710 j=1,2*n+3
            gap(i,j)=0
 710  continue

c Set equation at Psi points, assuming periodic b.c. in i.
c Cannot solve at both i=0 and i=imax as periodicity => would
c have singular matrix. At dry points equation is trivial.

      do 650 i=1,imax
         do 650 j=0,jmax
            k=i + j*n
            if(max(k1(i,j),k1(i+1,j),k1(i,j+1),k1(i+1,j+1)).le.kmax)then
               tv = (s(j+1)*rh(1,i,j+1) - s(j)*rh(1,i,j))
     &              /(2.*dsv(j)*dphi)
               tv1 = (sv(j)*rh(2,i+1,j) - sv(j)*rh(2,i,j))
     &              /(2.*dphi*dsv(j))

               gap(k,2) = drag(1,i,j)*c(j)*c(j)*rh(1,i,j)/(ds(j)*dsv(j))
     &                  + tv1

               l=n+1
c for periodic boundary in i
               if(i.eq.1)l=2*n+1

               gap(k,l) = drag(2,i,j)*rcv(j)*rcv(j)*rdphi*rdphi
     &             *rh(2,i,j) - tv

               gap(k,n+2) = - (drag(2,i,j)*rh(2,i,j)
     1             + drag(2,i+1,j)*rh(2,i+1,j))/(cv(j)*cv(j)*dphi
     1             *dphi) - (drag(1,i,j)*c(j)*c(j)*rh(1,i,j)/ds(j)
     1         + drag(1,i,j+1)*c(j+1)*c(j+1)*rh(1,i,j+1)/ds(j+1))/dsv(j)

               l=n+3
c for periodic boundary in i
               if(i.eq.imax)l=3

               gap(k,l) = drag(2,i+1,j)*rh(2,i+1,j)
     &                  /(cv(j)*cv(j)*dphi*dphi) + tv
               gap(k,2*n+2) = drag(1,i,j+1)*c(j+1)*c(j+1)*rh(1,i,j+1)
     &                      /(ds(j+1)*dsv(j)) - tv1

            else
               gap(k,n+2) = 1
            endif
 650  continue

* *
* now invert the thing
* *
      do 720 i=1,n*m-1
         im=min(i+n+1,n*m)
         do 725 j=i+1,im
            rat=gap(j,n+2-j+i)/gap(i,n+2)
            ratm(j,j-i)=rat
            if(rat.ne.0)then
               do 730 k=n+2-j+i,2*n+3-j+i
                  gap(j,k)=gap(j,k) - rat*gap(i,k+j-i)
 730           continue
            endif
 725     continue
 720  continue

      end


text:
Modeled on the successful Stories on Stage Sacramento series, Stories on Stage Davis aims to celebrate and share the remarkable writing and acting talent of Davis and the surrounding communities. The series features short fiction by established and emerging writers, introduced by their authors and read by actors.

The series, which is sponsored by the UC Davis Humanities Institute, is held at the Pence Gallery, 212 D Street, in downtown Davis, the second Saturday of every month, from 7:30pm to 9pm. Doors open at 7pm. Wine and light refreshments are available for purchase. Donation: $5. See the http://storiesonstagedavis.com/ Stories on Stage Davis homepage for information on upcoming events and bios of featured writers.

Stories on Stage Davis pays a small honorarium to the writers and actors they work with and is currently looking for both fiction to feature and actors to read. For fiction, short stories or standalone excerpts of novels, 3000 words or less, are preferred, but they will consider great stories up to 5000 words in length. While preference is given to authors living in Davis and surrounding areas, they are happy to look at great writing from anywhere.

Stories on Stage Davis is also accepting written resumes or video links from actors who are interested in reading and interpreting written works. Again, preference is given to local talent.

Submissions from both writers and actors can be uploaded at the http://storiesonstagedavis.com/?page_id6 submissions page at the Stories on Stage Davis website.


! Preconditions use to fail : CODE resource is unavailable
function preconditions()
end function preconditions


text:
Professional photography photographer. Weddings, Anniversaries, Commitment Ceremonies, Bnai Mitzvahs, Family Parties, Corporate Events, Individual & Family Portraits, and other events as requested.
Serving the Regional Sacramento Area. Travel outside Sacramento can be arranged.


text:
! RUN: %python %S/test_errors.py %s  %flang -fopenmp
! REQUIRES: shell
! Check OpenMP clause validity for NONTEMPORAL clause

program omp_simd
  integer i
  integer, allocatable :: a(:)

  allocate(a(10))

  !$omp simd nontemporal(a)
  do i = 1, 10
    a(i) = i
  end do
  !$omp end simd

  !$omp parallel do simd nontemporal(a)
  do i = 1, 10
    a(i) = i
  end do
  !$omp end parallel do simd

  !$omp parallel do simd nontemporal(a)
  do i = 1, 10
    a(i) = i
  end do
  !$omp end parallel do simd

  !ERROR: NONTEMPORAL clause is not allowed on the DO SIMD directive
  !$omp do simd nontemporal(a)
  do i = 1, 10
    a(i) = i
  end do
  !$omp end do simd

  !$omp taskloop simd nontemporal(a)
  do i = 1, 10
    a(i) = i
  end do
  !$omp end taskloop simd

  !$omp teams
  !$omp distribute parallel do simd nontemporal(a)
  do i = 1, 10
    a(i) = i
  end do
  !$omp end distribute parallel do simd
  !$omp end teams

  !$omp teams
  !$omp distribute simd nontemporal(a)
  do i = 1, 10
    a(i) = i
  end do
  !$omp end distribute simd
  !$omp end teams

  !$omp target parallel do simd nontemporal(a)
  do i = 1, 10
    a(i) = i
  end do
  !$omp end target parallel do simd

  !$omp target simd nontemporal(a)
  do i = 1, 10
    a(i) = i
  end do
  !$omp end target simd

  !$omp teams distribute simd nontemporal(a)
  do i = 1, 10
    a(i) = i
  end do
  !$omp end teams distribute simd

  !$omp teams distribute parallel do simd nontemporal(a)
  do i = 1, 10
    a(i) = i
  end do
  !$omp end teams distribute parallel do simd

  !$omp target teams distribute parallel do simd nontemporal(a)
  do i = 1, 10
    a(i) = i
  end do
  !$omp end target teams distribute parallel do simd

  !$omp target teams distribute simd nontemporal(a)
  do i = 1, 10
    a(i) = i
  end do
  !$omp end target teams distribute simd

end program omp_simd

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm


text:
c*************************************************************************
c                          ANAL_ENERGY_WRITE_HJS.F
c*************************************************************************
c Writes the energy of the total system (massive bodies) wrt time in
c   the hierarchical Jacobi case
c
c      Input:
c            t             ==>  current time
c            nbod          ==>  number of massive bodies (int scalar)
c            mass          ==>  mass of bodies (real array)
c            umat          ==>  Conversion matrix Jacobi => Barycentric
c                                (2D real array)
c            xj,yj,zj      ==>  current position in Jacobi coord
c                               (real arrays)
c            vxj,vyj,vzj   ==>  current velocity in Jacobi coord
c                               (real arrays)
c            iu            ==>  unit to write to (int scalar)
c            fopenstat     ==>  The status flag for the open
c                                statements of the output files.
c                                          (character*80)
c            eoff          ==> An energy offset that is added to the energy
c                                      (real*8 scalar)
c
c Remarks: Adapted from anal_energy_write & anal_energy
c Authors:  Herve Beust
c Date:    Sep 4, 2006

      subroutine anal_energy_write_hjs(t,nbod,umat,mass,xj,yj,zj,
     &     vxj,vyj,vzj,iu,fopenstat,diro,eoff)

      include '../swift.inc'

c...  Inputs:
      integer nbod,iu
      real*8 mass(nbod),t,umat(NPLMAX,NPLMAX),eoff
      real*8 xj(nbod),yj(nbod),zj(nbod)
      real*8 vxj(nbod),vyj(nbod),vzj(nbod)
      character*80 fopenstat,diro

c...  Internals
      integer i1st
      real*8 energy,eltot(3),ke,pot
      real*8 elx,ely,elz
      real*8 xx,yy,zz,rr2
      real*8 xb(NPLMAX),yb(NPLMAX),zb(NPLMAX)
      real*8 vxb(NPLMAX),vyb(NPLMAX),vzb(NPLMAX)
      integer i,j
      data i1st/0/
      save i1st

c----
c...  Executable code

c Compute and print initial ke,pot,energy and ang. mom.

      call coord_g2b(nbod,umat,mass,xj,yj,zj,vxj,vyj,vzj,
     &      xb,yb,zb,vxb,vyb,vzb)

      eltot(1)=(yb(nbod)*vzb(nbod)-zb(nbod)*vyb(nbod))*mass(nbod)
      eltot(2)=(zb(nbod)*vxb(nbod)-xb(nbod)*vzb(nbod))*mass(nbod)
      eltot(3)=(xb(nbod)*vyb(nbod)-yb(nbod)*vxb(nbod))*mass(nbod)

      ke = 0.5*mass(nbod)*(vxb(nbod)**2 + vyb(nbod)**2 + vzb(nbod)**2)
      pot = 0.d0

      do i = 1,nbod-1

         elx=(yb(i)*vzb(i)-zb(i)*vyb(i))*mass(i)
         ely=(zb(i)*vxb(i)-xb(i)*vzb(i))*mass(i)
         elz=(xb(i)*vyb(i)-yb(i)*vxb(i))*mass(i)
         eltot(1) = eltot(1) + elx
         eltot(2) = eltot(2) + ely
         eltot(3) = eltot(3) + elz

         ke = ke + 0.5*mass(i)*(vxb(i)**2 + vyb(i)**2 + vzb(i)**2)
         do j = i+1,nbod

	    xx = xb(i) - xb(j)
	    yy = yb(i) - yb(j)
	    zz = zb(i) - zb(j)
	    rr2 = xx**2 + yy**2 + zz**2
	    pot = pot - mass(i)*mass(j)/(sqrt(rr2))

         enddo
      enddo

      energy = ke + pot + eoff

      call io_energy_write(i1st,t,energy,eltot,iu,diro,fopenstat)

      if(i1st.eq.0) then
         i1st=1
      endif

      return
      end     !anal_energy_write_hjs
c--------------------------------------------------------------------------

max_stars_repo_name: Simske/exostriker


text:
!!!      subroutine thermal(ix)
!!!          use CN_Objects_manager
!!!          dimension ix(4,*)
!!!        do i=1,ElemCountAct
!!!            iElemconnect(i,1)=ix(1,i)
!!!            iElemconnect(i,2)=ix(2,i)
!!!            iElemconnect(i,3)=ix(3,i)
!!!            iElemconnect(i,4)=ix(4,i)
!!!        enddo
!!!
!!!!!!!!!!!!!c     this subroutine is to solve the thermal FEM problem
!!!!!!      integer ele,nod,ii,connect2(99999,4)
!!!!!!      dimension u(2*nod), yold(nod), zold(nod), y(nod), z(nod),
!!!!!!     1   yzold(2*nod),yz(2*nod),ix(4,*)
!!!!!!
!!!!!!      integer id(2*nod),maxint!,ix(4,ele)
!!!!!!      common /WMLthermali/connect(99999,4)
!!!!!!      common /WMLthermal/thermalflag,thermalconstraint(99999)
!!!!!!      common/bk08/kprint,nstep,ite,ilimit,newstf
!!!!!!      integer thermalconstraint,nstep
!!!!!!      integer connect
!!!!!!
!!!!!!!!!!!!!!!c      ele=fourele/4
!!!!!!!!!!!!!!!c      write(*,*)'in thermal.f'
!!!!!!!!!!!!!!!c      write(*,*)'ele=',ele
!!!!!!!!!!!!!!!c      write(*,*)'nod=',nod
!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!c      do i=1,ele
!!!!!!!!!!!!!!!c      	 connect(i,1)=ix(1,i)
!!!!!!!!!!!!!!!c        connect(i,2)=ix(2,i)
!!!!!!!!!!!!!!!c	 connect(i,3)=ix(3,i)
!!!!!!!!!!!!!!!c	 connect(i,4)=ix(4,i)
!!!!!!!!!!!!!!!c      enddo
!!!!!!
!!!!!!      do i=1,nod
!!!!!!      yzold(2*i-1)=yold(i)
!!!!!!      yzold(2*i)=zold(i)
!!!!!!      enddo
!!!!!!
!!!!!!      ii=1
!!!!!!      do i=1,2*nod
!!!!!!      if (id(i).eq.0) then
!!!!!!      yz(i)=yzold(i)
!!!!!!      ii=ii+1
!!!!!!      else
!!!!!!      yz(i)=yzold(i)+u(id(ii))   !ii could be replaced with do loop variable i
!!!!!!      ii=ii+1
!!!!!!      endif
!!!!!!      enddo
!!!!!!
!!!!!!      do i=1,nod
!!!!!!      y(i)=yz(2*i-1)
!!!!!!      z(i)=yz(2*i)
!!!!!!      enddo
!!!!!!
!!!!!!
!!!!!!
!!!!!!!!!!!!!      do i=1,2*nod
!!!!!!!!!!!!!c      write(*,*) u(i),y(i),z(i),id(i)
!!!!!!!!!!!!!      enddo
!!!!!!!!!!!!!c      write(*,*) 'before thermalFEMcall'
!!!!!!      maxint=4
!!!!!!      call thermalFEM(y,z,ele,nod,maxint)
!!!!!!
!!!
!!!      return
!!!      end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


text:
      subroutine mul66(rt,rtt,r)
      real*8 rt(6,6),r(6,6),rtt(6,6)
      do i=1,6
        do j=1,6
          do ia=1,6
            rtt(i,ia)=rt(i,j)*r(j,ia)+rtt(i,ia)
          end do
        end do
      end do
      end

      program test
      real*8 xj(6,6),w(6,6),w1(6,6)
      parameter(idump=0)
      integer i,j

      do i=1,6
        do j=1,6
          xj(i,j) = 0.0d0
          w1(i,j) = 0.0d0
          w(i,j) = i * 10.0d0 + j;
        end do
      end do

      xj(1,2) =  1.0d0
      xj(2,1) = -1.0d0
      xj(3,4) =  1.0d0
      xj(4,3) = -1.0d0
      xj(5,6) =  1.0d0
      xj(6,5) = -1.0d0

      call mul66(xj,w1,w)

      if (idump.ne.0) then
        write(6,*) 'w1 after call to mul66'
        do i = 1,6
          do j = 1,6
            write(6,'(D15.7)') w1(i,j)
          end do
        end do
      end if

      if (w1(1,1).ne.21.0d0) then
        call abort()
      end if

      end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
SUBROUTINE matsweep(g)

!-----------------------------------------------------------------
!
!  Directs the sweep across all computational cells for all
!   angles in all quadrants for constructing the matrices
!
!  For each cell, calls the weight, conab, cong, gammas routines
!
!  The calls jima and afcm to set construct the matrices
!
!-----------------------------------------------------------------

USE invar
USE solvar
IMPLICIT NONE
INTEGER, INTENT(IN) :: g
INTEGER :: oct, xs, ys, xe, ye, incx, incy, i, j, n, m
REAL*8 :: x, y, mu, eta, omega, sig, pi, tmp
REAL*8, DIMENSION(0:anord) :: tsxs

tmp = 0.0
pi = 2.0*ACOS(tmp)

! Allocate X-xmat and Y-ymat and Z-zmat and their 'old' counterparts
ALLOCATE(xmat(nmom,nx,ny), ymat(nmom,nx,ny,nx))
ALLOCATE(xold(nmom,nx,ny), yold(nmom,nx,ny))

! Allocate Boundary Condition Variables
ALLOCATE(ybcx(nx,nx),xbcx(ny,nx))
ALLOCATE(ybcy(nx),ybcxo(nx),ybcyo(nx))
ALLOCATE(xbcy(ny),xbcxo(ny),xbcyo(ny))

! Allocate gamma sub-matrices
ALLOCATE(gaa(nmom), gxa(nmom), gya(nmom))

! Start the sweep
! Start the loops in the quadrants
DO oct = 1, 4
   IF (oct == 1) THEN
      xs = 1
      xe = nx
      incx = 1
      ys = 1
      ye = ny
      incy = 1
   ELSE IF (oct == 2) THEN
      xs = nx
      xe = 1
      incx = -1
      ys = 1
      ye = ny
      incy = 1
   ELSE IF (oct == 3) THEN
      xs = nx
      xe = 1
      incx = -1
      ys = ny
      ye = 1
      incy = -1
   ELSE IF (oct == 4) THEN
      xs = 1
      xe = nx
      incx = 1
      ys = ny
      ye = 1
      incy = -1
   END IF

   ! Start the loop in all direction
   DO n = 1, apo
      mu  = ang(n,1)
      eta = ang(n,2)
      omega = omg(n)
      IF (incy == -1) omega = pi - omega

      ! Reset the Y matrix at beginning of new angle
      ymat = 0.0

      ! Reset the y-directed BC variables at start of a new angle
      ybcx = 0.0
      xbcx = 0.0

      ! Start the plane sweeps
      DO j = ys, ye, incy
         y = dy(j)

         ! Reset the X matrix at beginning of new row
         xmat = 0.0

         ! Reset the x-directed BC variables at start of new row
         ybcy = 0.0
         xbcy = 0.0

         ! Start the sweep in each row
         DO i = xs, xe, incx
            x = dx(i)

            ! Set up the spatial parameters
            m = mat(i,j)
            sig = sigt(m,g)
            tsxs = sigs(:,m,g,g)
            ex = mu/x
            ey = eta/y

            ! Get the elements of gamma
            CALL gammas(incx,mu,sig,omega)

            ! Call jima to set up matrices for phi and S vectors
            CALL jima(i,j,xs,ys,xe,ye,incx,incy,n,oct,mu,omega,tsxs)

            ! Call afcm to set up matrices for psi_in vector
            CALL afcm(i,j,xs,ys,xe,ye,incx,incy,n,oct,mu,omega)

            ! All matrices should be constructed, can exit sweep and solve
         ! End of the mesh sweep
         END DO                ! Rows x
      END DO                   ! Columns y
   END DO                      ! Angles
END DO                         ! End the loop over all octants

DEALLOCATE(xmat,xold,ymat,yold)
DEALLOCATE(ybcx,ybcy,ybcxo,ybcyo)
DEALLOCATE(xbcx,xbcy,xbcxo,xbcyo)
DEALLOCATE(gaa,gxa,gya)

RETURN
END SUBROUTINE matsweep

max_stars_repo_name: lanl/pidots

max_issues_repo_name: lanl/pidots

max_forks_repo_name: lanl/pidots


text:
      SUBROUTINE RCOVB
C
C     RCOVB PERFORMS THE BACK-SUBSTITUTIONS TO OBTAIN THE G-SET
C     DISPLACEMENTS OF A SUBSTRUCTURE WHOSE LEVEL IS LOWER THAN OR
C     EQUAL TO THAT OF THE FINAL SOLUTION STRUCTURE (FSS).
C     FOR EACH SUBSTRUCTURE WHOSE DISPLACEMENTS ARE RECOVERED,
C     AN SOLN ITEM IS CREATED BY EDITING THE SOLN ITEM OF THE FSS.
C
      EXTERNAL         ANDF
      LOGICAL          MODAL
      INTEGER          MCBTRL(7)  ,DRY        ,STEP       ,FSS       ,
     1                 RFNO       ,UINMS      ,SCHK       ,UA        ,
     2                 SSNM1      ,SYSBUF     ,RSP        ,RDP       ,
     3                 RECT       ,UPPER      ,LOWER      ,SYM       ,
     4                 HMCB       ,UBMCB      ,UAOMCB     ,UAMCB     ,
     5                 TFLAG      ,SIGNAB     ,SIGNC      ,SCRM      ,
     6                 UGV        ,UI(5)      ,SCR2       ,SCR3      ,
     7                 SCR5       ,NAME(2)    ,BLANK      ,UVEC      ,
     8                 POVE       ,HORG       ,SCR1       ,GMASK     ,
     9                 PAO        ,UB         ,IZ(1)      ,SOFSIZ    ,
     O                 SOF1       ,SOF2       ,SOF3       ,BUF1      ,
     1                 BUF2       ,RC         ,SSNM(2)    ,EQSS      ,
     2                 BUF(1)     ,ANDF       ,UIMPRO     ,ENERGY    ,
     3                 RMASK      ,FILE       ,RD         ,RDREW     ,
     4                 WRT        ,WRTREW     ,REW        ,EOFNRW    ,
     5                 BUF3       ,BUF4
C     INTEGER          SCR6       ,SCR7       ,SRD        ,SWRT
      DOUBLE PRECISION DZ(1)
      CHARACTER        UFM*23     ,UWM*25     ,UIM*29     ,SFM*25    ,
     1                 SWM*27
      COMMON /XMSSG /  UFM        ,UWM        ,UIM        ,SFM       ,
     1                 SWM
      COMMON /BLANK /  DRY        ,LOOP       ,STEP       ,FSS(2)    ,
     1                 RFNO       ,NEIGV      ,LUI        ,UINMS(2,5),
     2                 NOSORT     ,UTHRES     ,PTHRES     ,QTHRES
      COMMON /RCOVCR/  ICORE      ,LCORE      ,BUF1       ,BUF2      ,
     1                 BUF3       ,BUF4       ,SOF1       ,SOF2      ,
     2                 SOF3
      COMMON /RCOVCM/  MRECVR     ,UA         ,PA         ,QA        ,
     1                 IOPT       ,SSNM1(2)   ,ENERGY     ,UIMPRO    ,
     2                 RANGE(2)   ,IREQ       ,LREQ       ,LBASIC
      COMMON /SYSTEM/  SYSBUF     ,NOUT
      COMMON /NAMES /  RD         ,RDREW      ,WRT        ,WRTREW    ,
     1                 REW        ,NOREW      ,EOFNRW     ,RSP       ,
     2                 RDP        ,CSP        ,CDP        ,SQUARE    ,
     3                 RECT       ,DIAG       ,UPPER      ,LOWER     ,
     4                 SYM
      COMMON /MPYADX/  HMCB(7)    ,UBMCB(7)   ,UAOMCB(7)  ,UAMCB(7)  ,
     1                 MPYZ       ,TFLAG      ,SIGNAB     ,SIGNC     ,
     2                 MPREC      ,SCRM
      COMMON /ZZZZZZ/  Z(1)
      EQUIVALENCE      (BUF(1)    ,Z(1))
      EQUIVALENCE      (Z(1)      ,IZ(1)      ,DZ(1))
      DATA    NAME  /  4HRCOV,4HB          /
      DATA    UGV   ,  SCR1,SCR2,SCR3,SCR5 /
     1        106   ,  301, 302, 303, 305  /
      DATA    UI    /  204, 205, 206, 207, 208 /
      DATA    UVEC  ,  POVE,HORG,EQSS / 4HUVEC,4HPOVE,4HHORG,4HEQSS /
      DATA    IB    ,  SCHK / 1,  3   /
      DATA    SCR6  ,  SCR7,SRD,SWRT  / 306,307, 1,2 /
      DATA    RMASK /  469762048  /
      DATA    GMASK /  268435456  /
      DATA    MMASK /  134217728  /
      DATA    BLANK /  4H         /
C
C     INITIALIZE
C
      LCOREZ= KORSZ(Z) - LREQ
      SOF1  = LCOREZ - SYSBUF + 1
      SOF2  = SOF1 - SYSBUF - 1
      SOF3  = SOF2 - SYSBUF
      BUF1  = SOF3 - SYSBUF
      BUF2  = BUF1 - SYSBUF
      BUF3  = BUF2 - SYSBUF
      BUF4  = BUF3 - SYSBUF
      LCORE = BUF4 - 1
      IF (LCORE .LE. 0) GO TO 9008
      CALL SOFOPN (Z(SOF1),Z(SOF2),Z(SOF3))
      UA    = 0
      PAO   = 0
      TFLAG = 0
      SIGNAB= 1
      SIGNC = 1
      MPREC = 0
      SCRM  = SCR5
C
C     FIND OUT HOW MANY UI FILES THERE ARE AND WHICH ONES
C
      DO 10 I = 1,5
      IZ(1) = UI(I)
      CALL RDTRL (IZ)
      IF (IZ(1) .LT. 0) UINMS(1,I) = 0
   10 CONTINUE
C
C     IF UINMS(1,I) = 0         THEN  FILE UI(I) IS PURGED
C     IF UINMS(1,I) = BLANK     THEN  FILE UI(I) IS AVAILABLE AND NOT
C                                     IN USE
C     IF UINMS(1,I) = OTHER     THEN  FILE UI(I) CONTAINS UGV FOR
C                                     SUBSTRUCTURE -OTHER-
C
      SSNM(1) = SSNM1(1)
      SSNM(2) = SSNM1(2)
C
C     IF SSNM IS THE FINAL SOLUTION STRUCTURE (FSS), NO RECOVERY IS
C     NECESSARY.
C
      IF (SSNM(1).NE.FSS(1) .OR. SSNM(2).NE.FSS(2)) GO TO 190
      UA = UGV
      GO TO 508
C
C     SEARCH THE SOF FOR A DISPLACEMENT MATRIX OF SSNM OR A HIGHER
C     LEVEL SUBSTRUCTURE FROM WHICH THE REQUESTED DISPLACEMENTS CAN BE
C     RECOVERED
C
  190 JLVL = 1
  200 CALL SOFTRL (SSNM,UVEC,MCBTRL)
      RC  = MCBTRL(1)
      IF (RC .EQ. 1) GO TO 270
      IF (RC.EQ.2 .AND. DRY.LT.0) GO TO 270
      IF (RC .EQ. 3) GO TO 210
      IF (RC .EQ. 5) CALL SMSG (3,UVEC,SSNM)
      IF (RC .EQ. 4) GO TO 500
      WRITE (NOUT,63070) UWM,SSNM1,SSNM
      GO TO 9200
C
C     NO UVEC AT THIS LEVEL.  SAVE SSNM IN A STACK AT TOP OF OPEN CORE
C     AND SEARCH FOR UVEC OF THE NEXT HIGHER LEVEL
C
  210 LASTSS = 2*JLVL - 1
      IZ(LASTSS  ) = SSNM(1)
      IZ(LASTSS+1) = SSNM(2)
      JLVL = JLVL + 1
      CALL FNDNXL (Z(LASTSS),SSNM)
      IF (SSNM(1) .NE. BLANK) GO TO 230
      WRITE (NOUT,63060) UWM,IZ(LASTSS),IZ(LASTSS+1)
      GO TO 9200
  230 IF (SSNM(1).NE.IZ(LASTSS) .OR. SSNM(2).NE.IZ(LASTSS+1)) GO TO 240
      WRITE (NOUT,63080) UWM,SSNM1,SSNM
      GO TO 9200
C
C     IF SSNM IS NOT THE FSS, LOOK FOR UVEC ON THE SOF.  IF DRY RUN,
C     EXIT.  IF IT IS THE FSS, SET UA=UGV.  IF UGV IS NOT PURGED GO TO
C     BEGIN BACK-SUBSTITUTION.  OTHERWISE, GIVE IT THE SAME TREATMENT
C     AS IF IT WERE NOT THE FSS.
C
  240 IF (SSNM(1).NE.FSS(1) .OR. SSNM(2).NE.FSS(2)) GO TO 200
      IF (DRY .LT. 0) GO TO 500
      UA = UGV
      MCBTRL(1) = UA
      CALL RDTRL (MCBTRL)
      IF (MCBTRL(1) .GT. 0) GO TO 340
      GO TO 200
C
C     FOUND A UVEC ON SOF FOR THIS LEVEL.  SEE IF IT HAS ALREADY BEEN
C     PUT ON A UI FILE.  (IF DRY RUN, EXIT)
C
  270 IF (DRY.LT.0 .OR. JLVL.EQ.1) GO TO 500
      DO 280 I = 1,5
      UA = UI(I)
      IF (SSNM(1).EQ.UINMS(1,I) .AND. SSNM(2).EQ.UINMS(2,I)) GO TO 340
  280 CONTINUE
C
C     DATA BLANK /4H    /
C
C     IT DOES NOT RESIDE ON ANY UI FILE.  FIND A UI FILE TO USE.
C
      J = 0
      DO 290 I = 1,5
      IF (UINMS(1,I) .EQ. 0) GO TO 290
      J = J + 1
      IF (UINMS(1,I) .EQ. BLANK) GO TO 310
  290 CONTINUE
      GO TO 297
C
C     ALL UI FILES SEEM TO BE IN USE.  DO ANY REALLY EXIST
C
  297 IF (J .EQ. 0) GO TO 320
C
C     AT LEAST ONE EXISTS.  RE-USE THE ONE WITH OLDEST DATA
C
      I = LUI + 1
      IF (I .GT. 5) I = 1
      J = I
  300 IF (UINMS(1,I) .NE. 0) GO TO 310
C
C     NO FILE THERE.  TRY NEXT ONE.
C
      I = I + 1
      IF (I .GT. 5) I = 1
      IF (I .EQ. J) GO TO 320
      GO TO 300
C
C     FOUND A UI FILE TO USE
C
  310 LUI = I
      UA  = UI(I)
      UINMS(1,I) = SSNM(1)
      UINMS(2,I) = SSNM(2)
      GO TO 330
C
C     ALL UI FILES ARE PURGED.  USE SCR1 INSTEAD
C
  320 UA = SCR1
C
C     COPY UVEC FROM SOF TO UA
C
  330 CALL MTRXI (UA,SSNM,UVEC,0,RC)
C
C     TOP OF BACK-SUBSTITUTION LOOP
C
  340 UB = UA
      UAOMCB(1) = 0
      ICORE = LASTSS  + 2
      IDPCOR= ICORE/2 + 1
C
C     CHECK IF THE EQSS ITEM IS THERE FOR THIS SUBSTRUCTURE
C
      CALL SFETCH (Z(LASTSS),EQSS,SCHK,RC)
      IF (RC .NE. 1) GO TO 6317
C
C     COMPUTE TIME TO RECOVER THIS LEVEL AND CHECK TIME-TO-GO
C
C     (A DETAILED TIME CHECK SHOULD BE CODED LATER.  FOR THE PRESENT,
C     JUST CHECK TO SEE IF TIME HAS RUN OUT NOW.)
C
      CALL TMTOGO (I)
      IF (I .LE. 0) GO TO 6309
C
C     CHECK REMAINING SPACE ON SOF.  FIRST CALCULATE HOW MUCH SPACE
C     THE RECOVERED DISPLACEMENT MATRIX WILL TAKE (ASSUMING IT IS FULL).
C
      MCBTRL(1) = UB
      CALL RDTRL (MCBTRL)
      I = MCBTRL(2)
C
C     NO. OF COLUMNS IN DISPLACEMENT MATRIX IN I
C
      CALL SOFTRL (Z(LASTSS),HORG,MCBTRL)
      RC = MCBTRL(1)
      ITEM = HORG
      IF (RC .GT. 1) GO TO 6317
      NROW = MCBTRL(3)
      J = I*NROW
C
C     NOW CHECK SPACE
C
      IF (SOFSIZ(I) .LT. J) GO TO 6310
C
C     CREATE THE SOLUTION ITEM FOR THE RECOVERED SUBSTRUCTURE.
C
      CALL RCOVLS (Z(LASTSS))
      IF (IOPT .LT. 0) GO TO 9000
C
C     FIND A UI FILE FOR DISPLACEMENTS
C
      J = 0
      DO 420 I = 1,5
      IF (UINMS(1,I) .EQ. 0) GO TO 420
      J = J + 1
      IF (UINMS(1,I) .EQ. BLANK) GO TO 440
  420 CONTINUE
C
C     NO UNUSED UI FILES ARE AVAILABLE.  IF TWO OR MORE UI FILES ARE
C     NOT PURGED, USE THE ONE WITH OLDEST DATA.  OTHERWISE, USE SCR2.
C     MAKE SURE WE DON T ASSIGN THE SAME FILE AS THE HIGHER
C     LEVEL DISPLACEMENTS ARE ON (UB)
C
      IF (J .LT. 2) GO TO 450
      I = LUI + 1
      IF (I .GT. 5) I = 1
      J = I
  430 IF (UINMS(1,I).NE.0 .AND. UI(I).NE.UB) GO TO 440
      I = I + 1
      IF (I .GT. 5) I = 1
      IF (I .EQ. J) GO TO 450
      GO TO 430
C
C     FOUND A UI FILE
C
  440 LUI = I
      UA  = UI(I)
      UINMS(1,I) = IZ(LASTSS  )
      UINMS(2,I) = IZ(LASTSS+1)
      GO TO 455
  450 UA = SCR2
C
C     IF THE RECOVERED SUBSTRUCTURE WAS NOT REDUCED GENERATE THE
C     DISPLACEMENTS DIRECTLY.
C     IF THE SUBSTRUCTURE WAS REDUCED AND THE UIMPROVED FLAG IS SET
C     AND THIS IS A NON-STATICS RUN GENERATE THE IMPROVED DISPLACEMENTS.
C     IF THE SUBSTRUCTURE WAS IN A GUYAN REDUCTION AND THIS IS A
C     STATICS RUN GENERATE THE LOADS ON THE OMMITED POINTS.
C
C     INCLUDE THE CHECK ON THE POVE ITEM ALSO TO BE COMPATABLE WITH
C     PREVIOUS SOFS WITH NO TYPE BITS
C
  455 CALL SOFTRL (Z(LASTSS),POVE,MCBTRL)
      IPOVE = MCBTRL(1)
      CALL FDSUB (SSNM,IDIT)
      RC = 4
      IF (IDIT .LT. 0) GO TO 6317
      CALL FMDI (IDIT,IMDI)
      MODAL = .FALSE.
      IF (ANDF(BUF(IMDI+IB),MMASK) .NE. 0) MODAL = .TRUE.
      IF (ANDF(BUF(IMDI+IB),RMASK).NE.0 .AND. UIMPRO.NE.0 .AND.
     1    RFNO.GT.2) GO TO 470
      IF (ANDF(BUF(IMDI+IB),GMASK).NE.0 .AND. RFNO.LE.2) GO TO 480
      IF (ANDF(BUF(IMDI+IB),RMASK).EQ.0 .AND. IPOVE.EQ.1 .AND.
     1    RFNO.LE.2) GO TO 480
      GO TO 490
C
C     IF THE USER REQUESTED AN IMPROVED VECTOR AND THIS IS A NONSTATICS
C     RUN THEN GENERATE IT.
C
  470 CALL RCOVUI (UB,Z(LASTSS),MODAL)
      IF (IOPT .LT. 0) GO TO 9000
      GO TO 495
C
C     GENERATE THE LOADS ON THE OMITED POINTS FOR REDUCED SUBSTRUCTURES
C     IF THIS IS A STATICS RUN
C
  480 CALL RCOVUO (0,UAOMCB(1),Z(LASTSS))
      IF (IOPT .LT. 0) GO TO 9000
C
C     MULIPLY AND ADD TO GET DISPLACEMENTS OF LOWER-LIVEL SUBSTRUCTURE.
C
C     COPY H OR G TRANSFORMATION MATRIX TO SCR3
C
      CALL SOFOPN (Z(SOF1),Z(SOF2),Z(SOF3))
  490 ITEM = HORG
      CALL MTRXI (SCR3,Z(LASTSS),HORG,0,RC)
      IF (RC .NE. 1) GO TO 6317
C
C     SETUP FOR MPYAD
C
      CALL SOFCLS
      HMCB(1) = SCR3
      UBMCB(1)= UB
      CALL RDTRL (HMCB)
      CALL RDTRL (UBMCB)
      IF (UAOMCB(1) .NE. 0) CALL RDTRL (UAOMCB)
      CALL MAKMCB (UAMCB,UA,HMCB(3),RECT,UBMCB(5))
      MPYZ = LCOREZ - ICORE - 7
      CALL MPYAD (DZ(IDPCOR),DZ(IDPCOR),DZ(IDPCOR))
      CALL WRTTRL (UAMCB)
C
C     COPY RECOVERED DISPLACEMENTS TO SOF
C
  495 CALL SOFOPN (Z(SOF1),Z(SOF2),Z(SOF3))
      CALL MTRXO (UA,Z(LASTSS),UVEC,0,RC)
C
C     END OF BACK-SUBSTITUTION LOOP
C     CLOSE AND REOPEN THE SOF TO GET ANY CONTROL BLOCKS WRITTEN TO
C     FILE
C
      CALL SOFCLS
      CALL SOFOPN (Z(SOF1),Z(SOF2),Z(SOF3))
      SSNM(1) = IZ(LASTSS)
      SSNM(2) = IZ(LASTSS+1)
      LASTSS  = LASTSS - 2
      JLVL    = JLVL - 1
      WRITE (NOUT,63120) UIM,JLVL,SSNM
      IF (JLVL .GT. 1) GO TO 340
C
C     NORMAL COMPLETION OF MODULE EXECUTION
C
  508 CONTINUE
  500 DO 510 I = 1,5
      IF (UINMS(1,I) .EQ. 0) UINMS(1,I) = BLANK
  510 CONTINUE
      CALL SOFCLS
      RETURN
C
C     ERROR PROCESSING
C
 6309 WRITE (NOUT,63090) SFM,IZ(LASTSS),IZ(LASTSS+1),SSNM,SSNM1
      N = -37
      GO TO 9100
 6310 WRITE (NOUT,63100) SWM,IZ(LASTSS),IZ(LASTSS+1),SSNM,SSNM1
      GO TO 9200
 6317 IF (RC .EQ. 2) RC = 3
      CALL SMSG (RC-2,ITEM,Z(LASTSS))
 9000 WRITE (NOUT,63170) SWM,SSNM1
      GO TO 9200
 9008 N = 8
 9100 CALL SOFCLS
      CALL MESAGE (N,FILE,NAME)
 9200 IOPT = -1
      GO TO 500
C
C     FORMAT STATEMENTS
C
63060 FORMAT (A25,' 6306, ATTEMPT TO RECOVER DISPLACEMENTS FOR NON-',
     1       'EXISTANT SUBSTRUCTURE ',2A4)
63070 FORMAT (A25,' 6307, WHILE ATTEMPTING TO RECOVER DISPLACEMENTS ',
     1       'FOR SUBSTRUCTURE ',2A4,1H,, /32X,'THE DISPLACEMENTS FOR ',
     2       'SUBSTRUCTURE ',2A4,' WERE FOUND TO EXIST IN DRY RUN ',
     3       'FORM ONLY.')
63080 FORMAT (A25,' 6308, NO SOLUTION AVAILABLE FROM WHICH DISPLACE',
     1       'MENTS FOR SUBSTRUCTURE ',2A4, /32X,'CAN BE RECOVERED.  ',
     2       'HIGHEST LEVEL SUBSTRUCTURE FOUND WAS ',2A4)
63090 FORMAT (A25,' 6309, INSUFFICIENT TIME REMAINING TO RECOVER DIS',
     1       'PLACEMENTS OF SUBSTRUCTURE ',2A4, /32X,'FROM THOSE OF ',
     2       'SUBSTRUCTURE ',2A4,'.  (PROCESSING USER RECOVER REQUEST',
     3       /32X,'FOR SUBSTRUCTURE ',2A4,1H))
63100 FORMAT (A27,' 6310, INSUFFICIENT SPACE ON SOF TO RECOVER DIS',
     1       'PLACEMENTS OF SUBSTRUCTURE ',2A4, /32X,' FROM THOSE OF ',
     2       'SUBSTRUCTURE ',2A4,' WHILE PROCESSING USER RECOVER ',
     3       'REQUEST', /32X,'FOR SUBSTRUCTURE ',2A4)
63120 FORMAT (A29,' 6312, LEVEL',I4,' DISPLACEMENTS FOR SUBSTRUCTURE ',
     1       2A4, /36X,'HAVE BEEN RECOVERED AND SAVED ON THE SOF.')
63170 FORMAT (A25,' 6317, RECOVER OF DISPLACEMENTS FOR SUBSTRUCTURE ',
     1       2A4,' ABORTED.')
      END


text:
C
C	$Id: gwiiws.f,v 1.6 2008-07-27 00:21:07 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE GWIIWS
C
C  Initialization of constants.
C
C
C  Initialize workstation state list and other internal constants.
C
      include 'gwiwsl.h'
      include 'gwiins.h'
      include 'gwienu.h'
C
C  Internal constants.
C
C  MIN and MAX VDC address limits in the X and Y directions
C  (the VDC extent).
C
      MINXVD = 0
      MAXXVD = 32767
      MINYVD = 0
      MAXYVD = 32767
C
C  NDC to VDC mapping coefficients, X and Y directions.
C
      MXOFF = 0
      MXSCAL = 32767
      MYOFF = 0
      MYSCAL = 32767
C
C  WSL items.
C
C  State is inactive (0=GINACT)
C
      MSTATE = GINACT
C
C  Display surface is empty (1=GEMPTY)
C
      MDEMPT = GEMPTY
C
C  Requested workstation window in NDC is unit square.
C
      RWINDO(1) = 0.0
      RWINDO(2) = 1.0
      RWINDO(3) = 0.0
      RWINDO(4) = 1.0
C
C  Current workstation window in NDC is unit square.
C
      CWINDO(1) = 0.0
      CWINDO(2) = 1.0
      CWINDO(3) = 0.0
      CWINDO(4) = 1.0
C
C  Requested workstation viewport is real address space.
C
      RWKVP(1) = 0.
      RWKVP(2) = 32767.
      RWKVP(3) = 0.
      RWKVP(4) = 32767.
C
C  Current workstation viewport is real address space.
C
      CWKVP(1) = 0.
      CWKVP(2) = 32767.
      CWKVP(3) = 0.
      CWKVP(4) = 32767.
C
C  WSL clipping control parameters.
C
C  Clipping indicator is on (1=GCLIP)
C
      MRCLIP = GCLIP
C
C  Clipping rectangle.
C
      MRCREC(1) = MINXVD
      MRCREC(2) = MINYVD
      MRCREC(3) = MAXXVD
      MRCREC(4) = MAXYVD
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
! { dg-options "-O3 -ffast-math" }

      DIMENSION FPQR(25,25,25)
      INTEGER P,Q,R
            DO 130 R=1,N4MAX,2
               IF(P.GT.1) THEN
                  FPQR(P,Q,R)= RM2*FPQR(P,Q,R-2)*REC(P+Q+R-2)
               END IF
  130       RM2= RM2+TWO
      END


text:
cc ------------ dpmjet3.4 - authors: S.Roesler, R.Engel, J.Ranft -------
cc -------- phojet1.12-40 - authors: S.Roesler, R.Engel, J.Ranft -------
cc                                                      - oct'13 -------
cc ----------- pythia-6.4 - authors: Torbjorn Sjostrand, Lund'10 -------
cc ---------------------------------------------------------------------
cc                                  converted for use with FLUKA -------
cc                                                      - oct'13 -------

C...PYGANO
C...Evaluates the parton distributions of the anomalous photon,
C...inhomogeneously evolved from a scale P2 (where it vanishes) to Q2.
C...KF=0 gives the sum over (up to) 5 flavours,
C...KF<0 limits to flavours up to abs(KF),
C...KF>0 is for flavour KF only.
C...ALAM is the 4-flavour Lambda, which is automatically converted
C...to 3- and 5-flavour equivalents as needed.
C...Adapted from SaSgam library, authors G.A. Schuler and T. Sjostrand.

      SUBROUTINE PYGANO(KF,X,Q2,P2,ALAM,XPGA,VXPGA)

C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)

C...Local arrays and data.
      DIMENSION XPGA(-6:6), VXPGA(-6:6), ALAMSQ(3:5)
      DATA PMC/1.3D0/, PMB/4.6D0/, AEM/0.007297D0/, AEM2PI/0.0011614D0/

C...Reset output.
      DO 100 KFL=-6,6
        XPGA(KFL)=0D0
        VXPGA(KFL)=0D0
  100 CONTINUE
      IF(Q2.LE.P2) RETURN
      KFA=ABS(KF)

C...Calculate Lambda; protect against unphysical Q2 and P2 input.
      ALAMSQ(3)=(ALAM*(PMC/ALAM)**(2D0/27D0))**2
      ALAMSQ(4)=ALAM**2
      ALAMSQ(5)=(ALAM*(ALAM/PMB)**(2D0/23D0))**2
      P2EFF=MAX(P2,1.2D0*ALAMSQ(3))
      IF(KF.EQ.4) P2EFF=MAX(P2EFF,PMC**2)
      IF(KF.EQ.5) P2EFF=MAX(P2EFF,PMB**2)
      Q2EFF=MAX(Q2,P2EFF)
      XL=-LOG(X)

C...Find number of flavours at lower and upper scale.
      NFP=4
      IF(P2EFF.LT.PMC**2) NFP=3
      IF(P2EFF.GT.PMB**2) NFP=5
      NFQ=4
      IF(Q2EFF.LT.PMC**2) NFQ=3
      IF(Q2EFF.GT.PMB**2) NFQ=5

C...Define range of flavour loop.
      IF(KF.EQ.0) THEN
        KFLMN=1
        KFLMX=5
      ELSEIF(KF.LT.0) THEN
        KFLMN=1
        KFLMX=KFA
      ELSE
        KFLMN=KFA
        KFLMX=KFA
      ENDIF

C...Loop over flavours the photon can branch into.
      DO 110 KFL=KFLMN,KFLMX

C...Light flavours: calculate t range and (approximate) s range.
        IF(KFL.LE.3.AND.(KFL.EQ.1.OR.KFL.EQ.KF)) THEN
          TDIFF=LOG(Q2EFF/P2EFF)
          S=(6D0/(33D0-2D0*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
     &    LOG(P2EFF/ALAMSQ(NFQ)))
          IF(NFQ.GT.NFP) THEN
            Q2DIV=PMB**2
            IF(NFQ.EQ.4) Q2DIV=PMC**2
            SNFQ=(6D0/(33D0-2D0*NFQ))*LOG(LOG(Q2DIV/ALAMSQ(NFQ))/
     &      LOG(P2EFF/ALAMSQ(NFQ)))
            SNFP=(6D0/(33D0-2D0*(NFQ-1)))*LOG(LOG(Q2DIV/ALAMSQ(NFQ-1))/
     &      LOG(P2EFF/ALAMSQ(NFQ-1)))
            S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNFP-SNFQ)
          ENDIF
          IF(NFQ.EQ.5.AND.NFP.EQ.3) THEN
            Q2DIV=PMC**2
            SNF4=(6D0/(33D0-2D0*4))*LOG(LOG(Q2DIV/ALAMSQ(4))/
     &      LOG(P2EFF/ALAMSQ(4)))
            SNF3=(6D0/(33D0-2D0*3))*LOG(LOG(Q2DIV/ALAMSQ(3))/
     &      LOG(P2EFF/ALAMSQ(3)))
            S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNF3-SNF4)
          ENDIF

C...u and s quark do not need a separate treatment when d has been done.
        ELSEIF(KFL.EQ.2.OR.KFL.EQ.3) THEN

C...Charm: as above, but only include range above c threshold.
        ELSEIF(KFL.EQ.4) THEN
          IF(Q2.LE.PMC**2) GOTO 110
          P2EFF=MAX(P2EFF,PMC**2)
          Q2EFF=MAX(Q2EFF,P2EFF)
          TDIFF=LOG(Q2EFF/P2EFF)
          S=(6D0/(33D0-2D0*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
     &    LOG(P2EFF/ALAMSQ(NFQ)))
          IF(NFQ.EQ.5.AND.NFP.EQ.4) THEN
            Q2DIV=PMB**2
            SNFQ=(6D0/(33D0-2D0*NFQ))*LOG(LOG(Q2DIV/ALAMSQ(NFQ))/
     &      LOG(P2EFF/ALAMSQ(NFQ)))
            SNFP=(6D0/(33D0-2D0*(NFQ-1)))*LOG(LOG(Q2DIV/ALAMSQ(NFQ-1))/
     &      LOG(P2EFF/ALAMSQ(NFQ-1)))
            S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNFP-SNFQ)
          ENDIF

C...Bottom: as above, but only include range above b threshold.
        ELSEIF(KFL.EQ.5) THEN
          IF(Q2.LE.PMB**2) GOTO 110
          P2EFF=MAX(P2EFF,PMB**2)
          Q2EFF=MAX(Q2,P2EFF)
          TDIFF=LOG(Q2EFF/P2EFF)
          S=(6D0/(33D0-2D0*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
     &    LOG(P2EFF/ALAMSQ(NFQ)))
        ENDIF

C...Evaluate flavour-dependent prefactor (charge^2 etc.).
        CHSQ=1D0/9D0
        IF(KFL.EQ.2.OR.KFL.EQ.4) CHSQ=4D0/9D0
        FAC=AEM2PI*2D0*CHSQ*TDIFF

C...Evaluate parton distributions (normalized to unit momentum sum).
        IF(KFL.EQ.1.OR.KFL.EQ.4.OR.KFL.EQ.5.OR.KFL.EQ.KF) THEN
          XVAL= ((1.5D0+2.49D0*S+26.9D0*S**2)/(1D0+32.3D0*S**2)*X**2 +
     &    (1.5D0-0.49D0*S+7.83D0*S**2)/(1D0+7.68D0*S**2)*(1D0-X)**2 +
     &    1.5D0*S/(1D0-3.2D0*S+7D0*S**2)*X*(1D0-X)) *
     &    X**(1D0/(1D0+0.58D0*S)) * (1D0-X**2)**(2.5D0*S/(1D0+10D0*S))
          XGLU= 2D0*S/(1D0+4D0*S+7D0*S**2) *
     &    X**(-1.67D0*S/(1D0+2D0*S)) * (1D0-X**2)**(1.2D0*S) *
     &    ((4D0*X**2+7D0*X+4D0)*(1D0-X)/3D0 - 2D0*X*(1D0+X)*XL)
          XSEA= 0.333D0*S**2/(1D0+4.90D0*S+4.69D0*S**2+21.4D0*S**3) *
     &    X**(-1.18D0*S/(1D0+1.22D0*S)) * (1D0-X)**(1.2D0*S) *
     &    ((8D0-73D0*X+62D0*X**2)*(1D0-X)/9D0 +
     &    (3D0-8D0*X**2/3D0)*X*XL + (2D0*X-1D0)*X*XL**2)

C...Threshold factors for c and b sea.
          SLL=LOG(LOG(Q2EFF/ALAM**2)/LOG(P2EFF/ALAM**2))
          XCHM=0D0
          IF(Q2.GT.PMC**2.AND.Q2.GT.1.001D0*P2EFF) THEN
            SCH=MAX(0D0,LOG(LOG(PMC**2/ALAM**2)/LOG(P2EFF/ALAM**2)))
            XCHM=XSEA*(1D0-(SCH/SLL)**3)
          ENDIF
          XBOT=0D0
          IF(Q2.GT.PMB**2.AND.Q2.GT.1.001D0*P2EFF) THEN
            SBT=MAX(0D0,LOG(LOG(PMB**2/ALAM**2)/LOG(P2EFF/ALAM**2)))
            XBOT=XSEA*(1D0-(SBT/SLL)**3)
          ENDIF
        ENDIF

C...Add contribution of each valence flavour.
        XPGA(0)=XPGA(0)+FAC*XGLU
        XPGA(1)=XPGA(1)+FAC*XSEA
        XPGA(2)=XPGA(2)+FAC*XSEA
        XPGA(3)=XPGA(3)+FAC*XSEA
        XPGA(4)=XPGA(4)+FAC*XCHM
        XPGA(5)=XPGA(5)+FAC*XBOT
        XPGA(KFL)=XPGA(KFL)+FAC*XVAL
        VXPGA(KFL)=VXPGA(KFL)+FAC*XVAL
  110 CONTINUE
      DO 120 KFL=1,5
        XPGA(-KFL)=XPGA(KFL)
        VXPGA(-KFL)=VXPGA(KFL)
  120 CONTINUE

      RETURN
      END

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
! ###################################################################
! Copyright (c) 2013-2021, Marc De Graef Research Group/Carnegie Mellon University
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
!     - Redistributions of source code must retain the above copyright notice, this list
!        of conditions and the following disclaimer.
!     - Redistributions in binary form must reproduce the above copyright notice, this
!        list of conditions and the following disclaimer in the documentation and/or
!        other materials provided with the distribution.
!     - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
!        of its contributors may be used to endorse or promote products derived from
!        this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER 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.
! ###################################################################
module HDFsupportTest

contains

!--------------------------------------------------------------------------
!
! SUBROUTINE:HDFsupportTest_ExecuteTest
!
!> @author y
!
!> @brief
!
!> @note
!
!> @date
!--------------------------------------------------------------------------
subroutine HDFsupportTest_ExecuteTest(result) &
bind(c, name='HDFsupportTest_ExecuteTest') ! this routine is callable from a C/C++ program
!DEC$ ATTRIBUTES DLLEXPORT :: HDFsupportTest_ExecuteTest

use local
use HDF5
use typedefs
use HDFsupport
use,INTRINSIC :: ISO_C_BINDING

IMPLICIT NONE

integer(C_INT32_T),INTENT(OUT)               :: result

result = 0

end subroutine HDFsupportTest_ExecuteTest

!--------------------------------------------------------------------------
!
! SUBROUTINE:HDFsupportTest_TestRWChar
!
!> @author y
!
!> @brief
!
!> @note
!
!> @date
!--------------------------------------------------------------------------
subroutine HDFsupportTest_TestRWChar(result) &
bind(c, name='HDFsupportTest_TestRWChar') ! this routine is callable from a C/C++ program
!DEC$ ATTRIBUTES DLLEXPORT :: HDFsupportTest_TestRWChar

use local
use HDF5
use typedefs
use HDFsupport
use,INTRINSIC :: ISO_C_BINDING

IMPLICIT NONE

integer(C_INT32_T),INTENT(OUT)                :: result
character(fnlen)                              :: HDFfilename, tpn
character(len=1)                              :: EMsoftnativedelimiter
type(HDFobjectStackType)                      :: HDF_head
integer(kind=irg)                             :: hdferr
integer                                       :: dim1, dim2, dim3, dim4
integer(HSIZE_T)                              :: dims1(1), dims2(2), dims3(3), dims4(4)
integer*4                                     :: i1, i2, i3, i4
character(fnlen)                              :: dataset
character, DIMENSION(:), ALLOCATABLE          :: array1
character, DIMENSION(:, :), ALLOCATABLE       :: array2
character, DIMENSION(:, :, :), ALLOCATABLE    :: array3
character, DIMENSION(:, :, :, :), ALLOCATABLE :: array4

! determine the pathname delimiter character
EMsoftnativedelimiter = EMsoft_getEMsoftnativedelimiter()

! nullify the push/pop stack pointer
nullify(HDF_head%next)

! initialize the fortran HDF interface
CALL h5open_EMsoft(hdferr)

tpn = EMsoft_toNativePath(EMsoft_getEMsoftTestingPath())
HDFfilename = trim(tpn)//EMsoftnativedelimiter//'HDFsupport_Char.h5'//CHAR(0)
hdferr =  HDF_createFile(trim(HDFfilename), HDF_head)
if(hdferr.lt.0) then
  result = hdferr
  return
end if
dim1 = 10
dim2 = 10
dim3 = 10
dim4 = 20

dims1 = (/ dim1 /)
dims2 = (/ dim1, dim2 /)
dims3 = (/ dim1, dim2, dim3 /)
dims4 = (/ dim1, dim2, dim3, dim4 /)

write (*,*) "Before Allocate"
ALLOCATE (array1(dim1))
ALLOCATE (array2(dim1,dim2))
ALLOCATE (array3(dim1,dim2,dim3))
ALLOCATE (array4(dim1,dim2,dim3,dim4))

write (*,*) "Before Loops..."

do i1=0,dim1
  array1(i1+1) = char(i1)
  do i2=0,dim2
    array2(i1+1,i2+1) = char(i1 * i2)
    do i3=0,dim3
      array3(i1+1,i2+1,i3+1) = char(i1 * i2 * i3)
      do i4=0,dim4
        array4(i1+1,i2+1,i3+1,i4+1) = char(i1 * i2 * i3 * i4)
      end do
    end do
  end do
end do

write (*,*) "Before Writing"

dataset = 'char_1D'//CHAR(0)
hdferr = HDF_writeDatasetCharArray1D(dataset, array1, dims1, HDF_head)
if(hdferr.lt.0) then
  result = -1000
  DEALLOCATE (array1)
  DEALLOCATE (array2)
  DEALLOCATE (array3)
  DEALLOCATE (array4)
  return
end if

write (*,*) "Before DEALLOCATE"

DEALLOCATE (array1)
 DEALLOCATE (array2)
! DEALLOCATE (array3)
! DEALLOCATE (array4)

write (*,*) "Before HDF_pop"

call HDF_pop(HDF_head,.TRUE.)

write (*,*) "Before h5close_EMsoft"

! and close the fortran hdf interface
call h5close_EMsoft(hdferr)
if(hdferr.lt.0) then
  result = hdferr
end if

end subroutine HDFsupportTest_TestRWChar

end module HDFsupportTest

!HDF_writeEMheader
!HDF_push
!HDF_pop
!HDF_stackdump
!HDF_handleError
!-HDF_createFile
!-HDF_openFile
!-HDF_createGroup
!-HDF_openGroup
!-HDF_openDataset
!-HDF_writeDatasetTextFile
!-HDF_readfromTextfile
!HDF_readDatasetStringArray
!-HDF_extractDatasetTextfile
!HDF_writeDatasetStringArray
!HDF_writeDatasetCharArray1D
!HDF_writeDatasetCharArray2D
!HDF_writeDatasetCharArray3D
!HDF_writeDatasetCharArray4D
!-HDF_writeDatasetInteger
!-HDF_writeDatasetInteger1byteArray1D
!-HDF_writeDatasetIntegerArray1D
!-HDF_writeDatasetIntegerArray2D
!-HDF_writeDatasetIntegerArray3D
!-HDF_writeDatasetIntegerArray4D
!-HDF_writeDatasetFloat
!-HDF_writeDatasetDouble
!-HDF_writeDatasetFloatArray1D
!-HDF_writeDatasetFloatArray2D
!-HDF_writeDatasetFloatArray3D
!-HDF_writeDatasetFloatArray4D
!-HDF_writeDatasetDoubleArray1D
!-HDF_writeDatasetDoubleArray2D
!-HDF_writeDatasetDoubleArray3D
!-HDF_writeDatasetDoubleArray4D

!HDF_readDatasetCharArray1D
!HDF_readDatasetCharArray2D
!HDF_readDatasetCharArray3D
!HDF_readDatasetCharArray4D
!-HDF_readDatasetInteger
!-HDF_readDatasetIntegerArray1D
!-HDF_readDatasetIntegerArray2D
!-HDF_readDatasetIntegerArray3D
!-HDF_readDatasetIntegerArray4D
!-HDF_readDatasetFloat
!-HDF_readDatasetFloatArray1D
!-HDF_readDatasetFloatArray2D
!-HDF_readDatasetFloatArray3D
!-HDF_readDatasetFloatArray4D
!-HDF_readDatasetDouble
!-HDF_readDatasetDoubleArray1D
!-HDF_readDatasetDoubleArray2D
!-HDF_readDatasetDoubleArray3D
!-HDF_readDatasetDoubleArray4D

!!!!! HyperSlab Routines

max_stars_repo_name: josephtessmer/EMsoft

max_stars_repo_licenses: Intel Unlicense MIT-CMU BSD-3-Clause

max_issues_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense

max_forks_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense


text:
! This file is part of mctc-lib.
!
! 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.

module mctc_io_read_vasp
   use mctc_env_accuracy, only : wp
   use mctc_env_error, only : error_type, fatal_error
   use mctc_io_convert, only : aatoau
   use mctc_io_resize, only : resize
   use mctc_io_structure, only : structure_type, new
   use mctc_io_structure_info, only : structure_info
   use mctc_io_symbols, only : to_number, symbol_length
   use mctc_io_utils, only : getline
   implicit none
   private

   public :: read_vasp

   logical, parameter :: debug = .false.

contains

subroutine read_vasp(self, unit, error)

   !> Instance of the molecular structure data
   type(structure_type), intent(out) :: self

   !> File handle
   integer, intent(in) :: unit

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   logical :: selective, cartesian
   integer :: i, j, k, nn, ntype, natoms, izp, stat
   integer, allocatable :: ncount(:)
   real(wp) :: ddum, latvec(3), scalar, coord(3), lattice(3, 3)
   real(wp), allocatable :: xyz(:, :)
   character(len=:), allocatable :: line, comment
   character(len=2*symbol_length), allocatable :: args(:), args2(:)
   character(len=symbol_length), allocatable :: sym(:)
   type(structure_info) :: info

   selective = .false. ! Selective dynamics
   cartesian = .true.  ! Cartesian or direct
   lattice = 0
   stat = 0

   ntype = 0
   ! first line contains the symbols of different atom types
   call getline(unit, line, stat)
   if (stat /= 0) then
      call fatal_error(error, "Unexpected end of input encountered")
      return
   end if
   if (debug) print'(">", a)', line

   call parse_line(line, args, ntype)
   call move_alloc(line, comment)

   ! this line contains the global scaling factor,
   call getline(unit, line, stat)
   if (stat /= 0) then
      call fatal_error(error, "Unexpected end of input encountered")
      return
   end if
   if (debug) print'(">", a)', line
   read(line, *, iostat=stat) ddum
   if (stat /= 0) then
      call fatal_error(error, "Cannot read scaling factor from input")
      return
   end if
   ! the Ang->au conversion is included in the scaling factor
   if (debug) print'("->", g0)', ddum
   scalar = ddum*aatoau

   ! reading the lattice constants
   do i = 1, 3
      call getline(unit, line, stat)
      if (stat /= 0) then
         call fatal_error(error, "Unexpected end of lattice vectors encountered")
         return
      end if
      if (debug) print'("->", a)', line
      read(line, *, iostat=stat) latvec
      if (stat /= 0) then
         call fatal_error(error, "Cannot read lattice vectors from input")
         return
      end if
      lattice(:, i) = latvec * scalar
   end do
   ! Either here are the numbers of each element,
   ! or (>vasp.5.1) here are the element symbols
   call getline(unit, line, stat)
   if (stat /= 0) then
      call fatal_error(error, "Unexpected end of input encountered")
      return
   end if
   if (debug) print'(">", a)', line

   ! try to verify that first element is actually a number
   i = max(verify(line, ' '), 1)
   j = scan(line(i:), ' ') - 2 + i
   if (j < i) j = len_trim(line)

   ! CONTCAR files have additional Element line here since vasp.5.1
   if (verify(line(i:j), '1234567890') /= 0) then
      call parse_line(line, args, ntype)
      call getline(unit, line, stat)
      if (debug) print'("->", a)', line
      if (stat /= 0) then
         call fatal_error(error, "Unexpected end of input encountered")
         return
      end if
   else
      deallocate(comment)
   end if
   call parse_line(line, args2, nn)
   if (nn /= ntype) then
      call fatal_error(error, 'Number of atom types mismatches the number of counts')
      return
   end if

   allocate(ncount(nn), source = 0)
   do i = 1, nn
      read(args2(i), *, iostat=stat) ncount(i)
      izp = to_number(args(i))
      if (izp < 1 .or. ncount(i) < 1) then
         call fatal_error(error, "Unknown element '"//trim(args(i))//"' encountered")
         return
      end if
   end do

   natoms = sum(ncount)
   allocate(sym(natoms))
   allocate(xyz(3, natoms))

   k = 0
   do i = 1, nn
      do j = 1, ncount(i)
         k = k+1
         sym(k) = trim(args(i))
      end do
   end do

   call getline(unit, line, stat)
   if (stat /= 0) then
      call fatal_error(error, "Could not read POSCAR")
      return
   end if
   if (debug) print'(">", a)', line
   line = adjustl(line)
   if (line(:1).eq.'s' .or. line(:1).eq.'S') then
      selective = .true.
      call getline(unit, line, stat)
      if (debug) print'("->", a)', line
      if (stat /= 0) then
         call fatal_error(error, "Unexpected end of input encountered")
         return
      end if
      line = adjustl(line)
   end if

   cartesian = (line(:1).eq.'c' .or. line(:1).eq.'C' .or. &
      &         line(:1).eq.'k' .or. line(:1).eq.'K')
   do i = 1, natoms
      call getline(unit, line, stat)
      if (stat /= 0) then
         call fatal_error(error, "Unexpected end of geometry encountered")
         return
      end if
      if (debug) print'("-->", a)', line
      read(line, *, iostat=stat) coord
      if (stat /= 0) then
         call fatal_error(error, "Cannot read geometry from input")
         return
      end if

      if (cartesian) then
         xyz(:, i) = coord*scalar
      else
         xyz(:, i) = matmul(lattice, coord)
      end if

   end do

   ! save information about this POSCAR for later
   info = structure_info(scale=ddum, selective=selective, cartesian=cartesian)
   call new(self, sym, xyz, lattice=lattice, info=info)
   if (allocated(comment)) self%comment = comment

end subroutine read_vasp

subroutine parse_line(line, args, nargs)
   character(len=*), intent(in) :: line
   character(len=2*symbol_length), allocatable, intent(out) :: args(:)
   integer, intent(out) :: nargs
   integer, parameter :: p_initial_size = 50
   integer :: istart, iend
   allocate(args(p_initial_size), source=repeat(' ', 2*symbol_length))
   istart = 1
   iend = 1
   nargs = 0
   do while(iend < len_trim(line))
      istart = verify(line(iend:), ' ') - 1 + iend
      iend = scan(line(istart:), ' ') - 1 + istart
      if (iend < istart) iend = len_trim(line)
      if (nargs >= size(args)) then
         call resize(args)
      end if
      nargs = nargs + 1
      args(nargs) = trim(line(istart:iend))
   end do
end subroutine parse_line

end module mctc_io_read_vasp

max_stars_repo_name: loriab/mctc-lib

max_issues_repo_name: loriab/mctc-lib

max_forks_repo_name: loriab/mctc-lib


! Copyright (c) 2016 Alex Kramer <kramer.alex.kramer@gmail.com>
! See the LICENSE.txt file at the top-level of this distribution.

module vd

  use ieee_arithmetic

  use numerics, only: numerics_cmplx_phase, numerics_d1, numerics_d2
  use log, only: log_log_critical, log_stderr

  use precision, only: ip, fp

  implicit none

  private

  public :: vd_get_local_quantities
  public :: vd_get_indices
  public :: vd_validate_quantum_update

  public :: vd_obj
  type vd_obj
     ! Base class for VD point

     ! Total probability flux through virtual detector
     real(fp) :: net_flux

     ! Semi-classical binning flag
     logical :: semi_classical
     ! Number of standard deviations to include in quantum VD binning
     integer(ip) :: vd_np_stdev

     ! Whether detector is disjoint or not
     ! Non-disjoint detectors allocate their own momentum distribution count
     !   array
     ! Disjoint detectors are passed a pointer to an overall momentum
     !   distribution count array
     logical :: vd_disjoint

     ! Temporal grid step
     real(fp) :: dt

     ! Units
     real(fp) :: hbar
     real(fp) :: m

   contains
     procedure :: init_vd_base => vd_obj_init
  end type vd_obj

contains

  subroutine vd_obj_init(this, dt, sc, vd_disjoint, hbar, m, vd_np_stdev)
    ! Initialize base class. This helps to eliminate code reuse in subtypes
    !
    ! This method is exposed as this%init_vd_base
    ! In an inherited type, it is exposed as this%vd_obj%init_vd_base
    !
    ! this :: vd_obj instance
    ! dt :: temporal grid step
    ! sc :: semi-classical flag
    ! vd_disjoint :: disjoint VD flag
    ! hbar :: units
    ! m :: particle mass
    ! vd_np_stdev :: REQUIRED if sc = False, number of standard deviations to
    !   include in quantum VD binning
    class(vd_obj), intent(inout) :: this
    real(fp), intent(in) :: dt
    logical, intent(in) :: sc
    logical, intent(in) :: vd_disjoint
    real(fp), intent(in) :: hbar
    real(fp), intent(in) :: m
    integer(ip), intent(in), optional :: vd_np_stdev

    this%dt = dt
    this%semi_classical = sc
    this%hbar = hbar
    this%m = m
    this%vd_disjoint = vd_disjoint

    if (present(vd_np_stdev)) then
       this%vd_np_stdev = vd_np_stdev
    end if

    this%net_flux = 0.0_fp
  end subroutine vd_obj_init

  subroutine vd_get_indices(nx, nxl_ext, nxr_ext, nxl_vd, nxr_vd, xl_min, &
       xl_max, xr_min, xr_max)
    ! Get VD indices relative to total spatial grid
    !
    ! This method also checks whether the external and VD region sizes are
    ! coherent, and stops with an error if not.
    !
    ! nx :: spatial grid size
    ! nxl_ext :: number of left external points
    ! nxr_ext :: number of right external points
    ! nxl_vd :: number of left VD points
    ! nxr_vd :: number of right VD points
    ! xl_min :: left VD minimum index (returned value)
    ! xl_max :: left VD maximum index (returned value)
    ! xr_min :: right VD minimum index (returned value)
    ! xr_max :: right VD maximum index (returned value)
    integer, intent(in) :: nx
    integer, intent(in) :: nxl_vd, nxl_ext
    integer, intent(in) :: nxr_vd, nxr_ext
    integer, intent(out) :: xl_min, xl_max
    integer, intent(out) :: xr_min, xr_max

    ! Make sure external and VD region sizes make sense
    call vd_validate_spatial_input(nx, nxl_ext, nxr_ext, nxl_vd, nxr_vd)

    ! We have to be careful here, since Fortran arrays are 1-indexed
    xl_min = nxl_ext + 1
    xl_max = nxl_ext + nxl_vd

    ! Account for no virtual detectors
    ! We use a trick so that iterations like vd_xl_min - 1, vd-xl_mmax + 1
    ! don't actually iterate.
    ! Note that this means we *cannot* use array slices, we must use a do loop
    if (xl_min > xl_max) then
       xl_min = 0
       xl_max = -nx
    end if

    xr_max = nx - nxr_ext
    xr_min = xr_max - (nxr_vd - 1)

    if (xr_min > xr_max) then
       xr_min = 0
       xr_max = -nx
    end if

  end subroutine vd_get_indices

  subroutine vd_get_local_quantities(psi_arr, dx, m, hbar, p_mu, p_var, j)
    ! Get local p_mu, p_var, j quantities
    !
    ! psi_arr :: 3-element slice of wavefunction surround the virtual detector,
    !   in the desired spatial component direction. We require three elements
    !   because we use three-point finite differencing for derivatives.
    ! dx :: spatial component grid step
    ! m :: mass
    ! hbar :: units
    ! p_mu :: local component momentum value (to be returned)
    ! p_var :: local component momentum variance (to be returned)
    ! j :: local component flux (to be returned)
    complex(fp), intent(in) :: psi_arr(3)
    real(fp), intent(in) :: dx
    real(fp), intent(in) :: m
    real(fp), intent(in) :: hbar
    real(fp), intent(inout) :: p_mu
    real(fp), intent(inout) :: p_var
    real(fp), intent(inout) :: j

    real(fp) :: eps_fp

    integer(ip) :: i_x

    real(fp) :: phi_arr(3), mag_arr(3), p_arr(3), p_var_arr(3)
    complex(fp) :: z

    eps_fp = epsilon(1.0_fp)

    do i_x = 1, 3
       z = psi_arr(i_x)
       phi_arr(i_x) = numerics_cmplx_phase(z) * hbar
       mag_arr(i_x) = abs(z)**2
    end do

    ! Get local (mean) momentum
    call numerics_d1(phi_arr, p_arr, dx)

    p_mu = p_arr(2)

    ! Get local current
    j = mag_arr(2) / m * p_mu

    ! Calculate local momentum variance
    ! Note that this is numerically unstable!
    call numerics_d2(log(mag_arr), p_var_arr, dx)
    p_var = - hbar**2 / 4.0_fp * p_var_arr(2)

    ! There are cases in which the values of mag_arr are less than machine
    ! precision so that the derivative operation above is completely
    ! numerically unstable. In these cases, we attempt to solve this by setting
    ! the variance equal to the machine epsilon. This can be justified
    ! (informally) using L'Hopital's rule:
    !
    ! d^2/dx^2 log(rho) = 1 / rho^2 ( rho * rho'' - rho''^2)
    !
    ! If rho <= eps locally then rho' and rho'' <= rho, so
    ! d^2/dx^2 log(rho) <= eps, which means we can't have a better estimate
    ! than setting p_var to be the machine epsilon.
    if (ieee_is_nan(p_var) .or. (p_var .lt. eps_fp)) then
       p_var = eps_fp
    end if

  end subroutine vd_get_local_quantities

  subroutine vd_validate_spatial_input(nx, nxl_ext, nxr_ext, nxl_vd, nxr_vd)
    ! Validate spatial grid setup and exit abnormally if errors are found
    !
    ! nx :: number of spatial grid points
    ! nxl_ext :: number of left external grid points
    ! nxr_ext :: number of right external grid points
    ! nxl_vd :: number of left virtual detector points
    ! nxr_vd :: number of right virtuald detector points
    integer(ip) :: nx
    integer(ip) :: nxl_ext
    integer(ip) :: nxr_ext
    integer(ip) :: nxl_vd
    integer(ip) :: nxr_vd

    logical :: sane
    character(:), allocatable :: error_msg

    ! sane = (nxl_ext .gt. 0) .and. (nxr_ext .gt. 0)
    ! if (.not. sane) then
    !    error_msg = "No external region defined in numerical grid; "// &
    !         "stopping abnormally."
    !    call log_log_critical(error_msg, log_stderr)
    !    stop 0
    ! end if

    sane = (nxl_vd .gt. 0) .or. (nxr_vd .gt. 0)
    if (.not. sane) then
       error_msg = "At least one virtual detectors must be present in "// &
            "numerical grid; stopping abnormally."
       call log_log_critical(error_msg, log_stderr)
       stop 0
    end if

    sane = .true.

    if ( (nxl_vd .gt. 0 .and. nxl_ext .eq. 0) .or. &
         (nxr_vd .gt. 0 .and. nxr_ext .eq. 0) ) then
       sane = .false.
       error_msg = "At least 1 external point required for a virtual "// &
            "detector to be present; stopping abnormally."
       call log_log_critical(error_msg, log_stderr)
       stop 0
    end if

    sane = (nx .gt. nxl_ext + nxr_ext + nxl_vd + nxr_vd)
    if (.not. sane) then
       error_msg = "Numerical grid size greater than number of allocated "// &
            "points; no internal region possible; exiting abnormally."
       call log_log_critical(error_msg, log_stderr)
       stop 0
    end if

  end subroutine vd_validate_spatial_input

  subroutine vd_validate_quantum_update(i_p_min, i_p_max, np, valid)
    integer(ip), intent(inout) :: i_p_min, i_p_max
    integer(ip), intent(in) :: np
    logical, intent(inout) :: valid

    valid = .true.

    ! If p_min > this%p_max or p_max < this%p_min we don't do anything
    if ( (i_p_min .eq. 0_ip) .or. (i_p_max .eq. -1_ip) ) then

       valid = .false.

    else

       ! If p_mmin < this%p_min, we start at this%p_min
       if (i_p_min .eq. -1_ip) then
          i_p_min = 1_ip
       end if

       ! If p_max > this%p_max, we end at this%p_max
       if (i_p_max .eq. 0_ip) then
          i_p_max = np
       end if

    end if
  end subroutine vd_validate_quantum_update

end  module vd