text:
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!  BIGSTICK configuration interaction shell-model code
!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  BBASIS_HAIKUS.f
!
!  third set of routines to generate the basis:
!  HAIKUS -- sorting
!
! SUBROUTINES:
! sort_haikus -- sorts haikus by parity, Jz
! test_haiku_qns
!	calls: convert_haiku, jz_of_occ, par_of_occ
!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc
! NOTE: it is not clear that sorting is really needed, as I later
! reorganize the haikus into blocks. This only slows down the code
! for very large cases, so it is acceptable for now
!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

!---------------------------------------------------------------------------
!  list of haikus
!---------------------------------------------------------------------------
module haikus
  implicit none

  type haikulist
     integer          :: nhaiku
     integer, pointer :: par(:),jz(:)
     integer, pointer :: hsd(:,:)
  end type haikulist

  type haikuWlist
     integer          :: nw        ! # of ws allowed
     integer, pointer :: wlist(:)  ! list of allowed ws
     integer, pointer :: wmap(:)   ! map of w-value onto wlist
     type (haikulist), pointer :: of_w(:)
  end type haikuWlist

  type haikuNlist
     type ( haikuWlist ),pointer :: of_nh(:)
  end type haikuNlist

  type ( haikuNlist   ) :: haiku(-2:2)

contains
!=====
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!  BIGSTICK configuration-interaction shell-model code
!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  BBASISLIB2.f
!
!  first set of routines to generate the basis:
!  HAIKUS
!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  NOTES:  The goal of these routines is to generate the "haikus."
!  Haikus are half-Slater-determinants: left-haikus are constructed from
!  s.p. states with jz < 0, and right-haikus with jz >= 0. Combining
!  them, one constructs the proton and neutron Slater determinants.
!  There are usually far fewer haikus than proton and neutron SDs.
!
!  Haikus are constructed recursively using the "parent" information for
!  templates, which are already organized by nh (particle number) and W.
!  (After this point the templates can be discarded.)
!  Once constructed, they are sorted by parity and jz and then organized
!  into "blocks" which are labeled by the quantum numbers:  nh, wh, parh, and jzh.
!C  The location of an individual haiku within a block is its "address."
!C  Later, "hops" will encode the application of creation/destruction operators
!C  on haikus, listed by initial and final block/address. After this point, the
!C  representation of haikus as binary words can be discarded.  Because we
!C  only use the haikus briefly, we don't have to worry about fetching them from
!C  from memory efficiently and can point to a list out-of-order.
!C
!C  Step 1: For each nh, generate all haikus for each allowed W, using templates
!C  Step 2: For each nh, wh, sort by parity and then by jz
!C  Step 3: Compute haiku blocks (lists of nh-wh-parh-jzh and # of haikus)
!C  Step 4: create block_lists
!C
!C  In order to know how to effectively store haikus, one must think ahead to how
!c  the proton / neutron SDs will be constructed from haikus.
!C
!C  The proton / neutron SDs are organized by sectors, labeled by Jz, parity, W.
!C  Each sector is defined by a list of paired haiku blocks.
!C
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!C
!C  SUBROUTINES :
!C   ** haiku master **:  master routine to call others
!C     calls subroutines get_templates_ws
!C                    :  count_create_haikus:
!C        		haiku_qns :
!C  (** test_haikus **   for testing only)
!   ** get_templates_ws **  extracts info on Ws from templates
!   ** count_create_haikus ** recursively creates haikus
!C     calls subroutines  convert_haiku
!                        add_haiku_bit  :
!     calls functions    bit_occ
!   **  add_haiku_bit **  adds a bit to a haiku
!     call subroutine  convert_loc_word
!   **  convert_loc_word ** converts the loc of a state to a word and specific bit
!   **  haiku_qns   **   computes q#s (jz, parity) of a set of haikus
!     calls subroutine convert_haiku
!     call functions  jzh_of_occ
!                     par_of_occ
!   ** convert_haiku **  converts a haiku in binary form to an occupation array
!     call function  bit_occ
!
!  (** convert_haiku_array **  converts a haiku in binary form to an occupation array
!                        does not assume a fix # of particles; seldom used)

!  FUNCTIONS:
!   ** bit_occ **        sees if a bit (state) in a haiku is occupied
!   ** jzh_of_occ **   computes jzh of an occ
!   ** par_of_occ0 **   computes parity of an occ
!   ** parmult  ** multiples two parities of the form 1,2
!
!  (OBSOLETE OR SELDOM USED FUNCTIONS )
!   ** nh_of_occ **     computes nh of an occ
!   ** jzh_of_occ0 **   computes jzh of an occ, checking ALL bits
!   ** wh_of_occ0 **   computes jzh of an occ, checking ALL bits
!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

!  subroutine haiku_master
!
!  coordinates creation of haikus
!
!  INPUT: ith = species/handness
!
!  SUBROUTINES CALLED:
!    get_templates_ws : extracs information on Ws allowed for fixed nh
!    count_create_haikus: creates the haikus recursively
!    haiku_qns : computes haiku quantum #s
!
!
subroutine haiku_master(ith)

use bitstuff
use haiku_info
use templates
use verbosity
implicit none

integer ith
integer it
integer nh
integer nw,whmin,whmax
integer, pointer :: wlist(:),wmap(:)
integer iw,wh
integer nhaiku,ntemp
integer, allocatable :: occ(:)
integer :: aerr

it = abs(ith)

allocate(occ(nhsps(abs(ith))), stat=aerr)
if(aerr /= 0) call memerror("haiku_master 1")
allocate(haiku(ith)%of_nh(0:maxNh(ith)), stat=aerr)
if(aerr /= 0) call memerror("haiku_master 2")

!------------- DERIVE THE HAIKU for zero particles

haiku(ith)%of_nh(0)%nw = 1
allocate(haiku(ith)%of_nh(0)%wlist(1), stat=aerr)
if(aerr /= 0) call memerror("haiku_master 3")
allocate(haiku(ith)%of_nh(0)%wmap(0:0), stat=aerr)
if(aerr /= 0) call memerror("haiku_master 4")
allocate(haiku(ith)%of_nh(0)%of_w(1), stat=aerr)
if(aerr /= 0) call memerror("haiku_master 5")
haiku(ith)%of_nh(0)%wlist(1) = 0
haiku(ith)%of_nh(0)%wmap(0) = 1
haiku(ith)%of_nh(0)%of_w(1)%nhaiku = 1

allocate(haiku(ith)%of_nh(0)%of_w(1)%hsd(nword(it),1 ), stat=aerr)
if(aerr /= 0) call memerror("haiku_master 6")
haiku(ith)%of_nh(0)%of_w(1)%hsd = 0
call haiku_qns(ith,0,1,occ)
!------------- loop over nh
do nh = 1,maxNh(ith)
!------------- FIGURE OUT W INFORMATION --------------

   call get_template_ws(ith,nh,nw,wlist,whmin,whmax,wmap)

!------------- ALLOCATE ------------------------
   haiku(ith)%of_nh(nh)%nw = nw
   allocate(haiku(ith)%of_nh(nh)%wlist(nw), stat=aerr)
   if(aerr /= 0) call memerror("haiku_master 10")
   allocate(haiku(ith)%of_nh(nh)%wmap(whmin:whmax), stat=aerr)
   if(aerr /= 0) call memerror("haiku_master 11")
   do iw = 1,nw
      haiku(ith)%of_nh(nh)%wlist(iw) = wlist(iw)
   enddo !iw
   do iw = whmin,whmax
      haiku(ith)%of_nh(nh)%wmap(iw) = wmap(iw)
   enddo  !iw

   allocate(haiku(ith)%of_nh(nh)%of_w(nw), stat=aerr)
   if(aerr /= 0) call memerror("haiku_master 12")

   do iw = 1,nw
      wh = wlist(iw)
!-------------- COUNT UP # OF HAIKUS
      call count_create_haikus(.false.,ith,nh,wh,iw,nhaiku)
!------------- ALLOCATE --------------------
      haiku(ith)%of_nh(nh)%of_w(iw)%nhaiku = nhaiku

      allocate(haiku(ith)%of_nh(nh)%of_w(iw)%hsd(nword(it),nhaiku), stat=aerr)
      if(aerr /= 0) call memerror("haiku_master 13")

!--------------- STORE THE HAIKUS ----------------------
      call count_create_haikus(.true.,ith,nh,wh,iw,nhaiku)
!--------------- COMPUTE QUANTUM NUMBERS --------------------
      call haiku_qns(ith,nh,iw,occ)
!----------------- SORT -------------------------------
   enddo  !iw
enddo ! nh
do nh = 1,maxNh(ith)
    nw = haiku(ith)%of_nh(nh)%nw
enddo  ! nh
deallocate(occ)
!---------------------- WRITE OUT IF VERBOSE --------------

if(verbose_haikus)then
   nhaiku = 0
   print*,' '
   print*,' species ',ith
   do nh = 0,maxNh(ith)
          ntemp = 0
         nw = haiku(ith)%of_nh(nh)%nw
         do iw = 1,nw
           ntemp = ntemp + haiku(ith)%of_nh(nh)%of_w(iw)%nhaiku
         enddo  !iw
         print*,' For nh ',nh,' there are ',ntemp,' haikus '
         nhaiku = nhaiku + ntemp
   enddo
   print*,' There are a total of ',nhaiku,' haikus '
   print*,' '
endif
return
end subroutine haiku_master
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  for testing purposes only
!
!  CALLS:
!   convert_haiku_array
!   nh_of_occ
!   wh_of_occ
!
      subroutine test_haikus(ith)

      use haiku_info
      implicit none
!      include 'binterfaces.inc'

      integer ith
      integer nh
      integer wh,iwh,nw
      integer nhaiku,ih

      integer,allocatable :: occ(:)
      integer,pointer  :: hsd(:)
      integer :: aerr

!      integer :: nh_of_occ,wh_of_occ

      allocate(occ(nhsps(abs(ith))), stat=aerr)
      if(aerr /= 0) call memerror("test_haikus 1")

      do nh = 0,maxNh(ith)
        nw = haiku(ith)%of_nh(nh)%nw
        do iwh = 1,nw
          wh = haiku(ith)%of_nh(nh)%wlist(iwh)
          nhaiku = haiku(ith)%of_nh(nh)%of_w(iwh)%nhaiku
          do ih = 1,nhaiku
           hsd => haiku(ith)%of_nh(nh)%of_w(iwh)%hsd(:,ih)
           call convert_haiku_array(ith,hsd,occ)
103      format(60i1)
           if(nh /= nh_of_occ(ith,occ))then
              print*,' nh problem ',nh,nh_of_occ(ith,occ),iwh,ih

           endif
           if(wh /= wh_of_occ(ith,nh,occ))then
              print*,' wh problem ',wh,ih
           endif

          enddo

        enddo ! iwh
      enddo ! nh

      deallocate(occ)
      return
      end subroutine test_haikus

! subroutine get_template_ws MOVED in 7.6.9 to bbasis_template.f90 containing the module templates
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  subroutine count_create_haikus
!
!  counts up and creates haikus recursively, using "parents" in templates
!
!  INPUT:
!   create : logical flag to create (.true.) or just count up (.false.)
!   ith = species/handness
!   nh = # of particles
!   wh = w-value of haikus
!   iwh = label of which wh for haikus
!
! OUTPUT (if create = .false.)
!  nhaiku = # of haikus for this nh,wh
!
! SUBROUTINES CALLED:
!  convert_haiku  : converts a haiku (binary word) to occupation array (occ)
!  add_haiku_bit  : directly adds a bit to a haiku
!
! FUNCTIONS CALLED:
!  bit_occ = .true. if a single particle state is occupied
!

      subroutine count_create_haikus(create,ith,nh,wh,iwh,nhaiku)

      use spsections
      use templates
      use haiku_info
      implicit none

!.....................................................
      logical create
      integer ith
      integer nh
      integer wh,iwh
      integer nhaiku
!......................
      integer it
      integer i
      integer parent,wparent,nparent,iwparent
      integer  ih
      integer padd
      integer pstart,pfin
      integer p
      integer iword, ibit
      integer hstart,hfin
      integer,pointer :: hsd(:),hsd2(:)
!      integer addbit
!      integer nh_of_occ
      integer occ(nhsps(abs(ith)))
      it = abs(ith)
      nhaiku = 0

      do i = 1,ntemplate(ith,nh)
        if(template(ith)%of_nh(nh)%list(i)%w < wh)cycle
        if(template(ith)%of_nh(nh)%list(i)%w > wh)exit

       if(.not.create)template(ith)%of_nh(nh)%list(i)%hstart = nhaiku+1

        padd    = template(ith)%of_nh(nh)%list(i)%add
        parent = template(ith)%of_nh(nh)%list(i)%parent

        if(nh == 1)then
           wparent = 0
           iwparent = 1
           nparent = 1
           hstart = 1
           hfin = 1
        else
           wparent = template(ith)%of_nh(nh-1)%list(parent)%w
           iwparent = haiku(ith)%of_nh(nh-1)%wmap(wparent)
           hfin = template(ith)%of_nh(nh-1)%list(parent)%hfin
           hstart = template(ith)%of_nh(nh-1)%list(parent)%hstart
        endif

        pstart = section(ith)%list(padd)%start
        pfin = section(ith)%list(padd)%fin
        do ih = hstart,hfin

          hsd => haiku(ith)%of_nh(nh-1)%of_w(iwparent)%hsd(:,ih)

          do p = pfin,pstart,-1  ! loop over and find open space
             if(.not.bit_occ(hsd,p))then
                nhaiku = nhaiku+1
                if(create)then
                   haiku(ith)%of_nh(nh)%of_w(iwh)%hsd(:,nhaiku)=      haiku(ith)%of_nh(nh-1)%of_w(iwparent)%hsd(:,ih)
       call add_haiku_bit(haiku(ith)%of_nh(nh)%of_w(iwh)%hsd(:,nhaiku),nword(it),p)

!------------------------ TESTING ------------------------
          hsd2 => haiku(ith)%of_nh(nh)%of_w(iwh)%hsd(:,nhaiku)
           call convert_haiku(ith,nh,hsd2,occ)

                endif
              else
                exit
             endif
          enddo  ! p
        enddo  ! ih
        if(.not.create)then  ! fill in # of haikus per template
          template(ith)%of_nh(nh)%list(i)%hfin = nhaiku
        endif
101     format('so far ',6i3)
      enddo ! i

      return
      end subroutine count_create_haikus

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  function bit_occ
!
!  determines if the isps state is occupied
!  NB: this function is called a LOT; some tweaking to optimize has been done
!  but possibly could be improved (however not currently a significant bottleneck)
!
!  INPUT:
!   hsd(:) = pointer of binary words that make up the haiku
!   isps = which s.p. state
!
!  OUTPUT:
!   bit_occ = .true. if there is a bit at location isps
!
      function bit_occ(hsd,isps)
      use bitstuff
      implicit none
      logical bit_occ
      integer,pointer :: hsd(:)
      integer isps

      integer iword,i
      integer imask

      iword = (isps-1)/(max_bit_word) + 1
      i = isps - (iword -1)*(max_bit_word)
      imask = ishft(1,i-1)
	  bit_occ = .not.(iand(imask,hsd(iword)) == 0)   ! should be faster than an if statement
!      if(iand(imask,hsd(iword)) == 0)then
!        bit_occ = .false.
!      else
!        bit_occ = .true.
!      endif
      return
      end function bit_occ

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  subroutine add_haiku_bit
!    adds a bit to a haiku to represent an occupied state
!
! INPUT:
!   hsd(:) = array of binary words = haiku
!   n = nwords in haiku
!   isps = location of state
!
!  subroutines called: convert_loc_word
!
      subroutine add_haiku_bit(hsd,n,isps)
      implicit none
      integer :: hsd(n)
      integer isps
      integer n
      integer i
      integer iword
      integer imask
      call convert_loc_word(isps,iword,i)
      imask = ishft(1,i-1)
      hsd(iword) = hsd(iword) + imask
      return
      end subroutine add_haiku_bit

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  function addbit
!
!  adds a bit in the isps'th location to a binary word (haiku)
!
!  INPUT: isps  = which s.p. state
!
!  OUTPUT: addbit = binary word to be added directly
!
!  SUBROUTINES CALLED:
!   convert_loc_word: computes where the word, bit of the state is located
!
      integer function addbit(isps)
      implicit none
      integer isps

      integer iword
      integer i
      call convert_loc_word(isps,iword,i)
      addbit = ishft(1,i-1)
      return
      end function addbit

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  subroutine convert_loc_word
!
!  converts a location in abstract single-particle space to
!  an absolute location in terms of words and bits
!
!  each haiku is stored in binary over several words
!
!  the location of a single-particle state is which bit on which word
!  We have to map an abstract index for the state (e.g., 29, 44, 83)
!
!  for example, suppose we have 31 bits per "word"
!  If isps = 29, then iword = 1,i = 29
!  If isps = 44, then iword = 2,i = 13
!  If sps = 83,  then iword = 3,i = 11
!
!  INPUT: isps  = abstract location
!
!  OUTPUT:
!         iword = which word
!         i     = which bit of iword
!

      subroutine convert_loc_word(isps,iword,i)
      use bitstuff
      implicit none
      integer isps,iword,i

      iword = (isps-1)/(max_bit_word) + 1
      i = isps - (iword -1)*(max_bit_word)
      return

      end subroutine convert_loc_word
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!      subroutine haiku_qns
!
!   computes the quantum numbers jz, wh, parh for a list of haikus
!   assumes the haikus have already been segregated by nh, w
!
!  INPUT:
!    ith  = species/handness
!    nh = # of particles, iwh = which W value
!    occ(:)  an array used repeatedly
!
!  SUBROUTINES CALLED:
!    convert_haiku : converts a haiku (hsd) to occupation array (occ)
!
!  FUNCTIONS CALLED:
!    jzh_of_occ, par_of_occ
!
      subroutine haiku_qns(ith,nh,iwh,occ)

      use haiku_info
      implicit none

      integer ith
      integer nh
      integer iwh
      integer :: occ(:)
!............................................
      integer nhaiku
      integer ih
      integer,pointer :: hsd(:)
      integer :: aerr
!      integer jzh_of_occ,par_of_occ

      nhaiku = haiku(ith)%of_nh(nh)%of_w(iwh)%nhaiku
      allocate(haiku(ith)%of_nh(nh)%of_w(iwh)%par(nhaiku), stat=aerr)
      if(aerr /= 0) call memerror("haiku_qns 1")
      allocate(haiku(ith)%of_nh(nh)%of_w(iwh)%jz(nhaiku), stat=aerr)
      if(aerr /= 0) call memerror("haiku_qns 2")
!$omp parallel do private(ih,hsd,occ)
      do ih = 1,nhaiku
         hsd => haiku(ith)%of_nh(nh)%of_w(iwh)%hsd(:,ih)
         call convert_haiku(ith,nh,hsd,occ)
         haiku(ith)%of_nh(nh)%of_w(iwh)%jz(ih) = jzh_of_occ(ith,nh,occ)
         haiku(ith)%of_nh(nh)%of_w(iwh)%par(ih) =  par_of_occ(ith,nh,occ)
      enddo
!$omp end parallel do
      return
      end subroutine haiku_qns

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  subroutine convert_haiku_array
!
!  converts a haiku  (binary word) to an array
!
!  INPUT:
!   ith = species/handness
!   hsd(:) = pointer to binary words that make up the haiku
!
!  OUTPUT:
!    occ(:)  array = 0. or 1 designating occupied s.p. states
!
!  FUNCTIONS CALLED:
!  bit_occ: computes whether or not a specific bit is occupied
!
      subroutine convert_haiku_array(ith,hsd,occ)

      use haiku_info
      implicit none

      integer :: ith
      integer,pointer :: hsd(:)
      integer :: occ(nhsps(abs(ith)))
      integer isps

      occ = 0
      do isps = 1,nhsps(abs(ith))
        if(bit_occ(hsd,isps))occ(isps) =1
      enddo
      return
      end subroutine convert_haiku_array
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  subroutine convert_haiku
!  converts a haiku  (binary word) to an array ASSUMING fixed nh
!  faster because it exits
!  note this has a different definition of occ from convert_haiku_array
!
!  INPUT:
!   ith = species/handness
!   nh = # of paricles
!   hsd(:) = pointer to binary words that make up the haiku
!
!  OUTPUT:
!    occ(:)  array = 0. or 1 designating occupied s.p. states
!
!
      subroutine convert_haiku_OLD(ith,nh,hsd,occ)
      use bitstuff
      use haiku_info
      implicit none

      integer :: ith
      integer :: nh
      integer,pointer :: hsd(:)
      integer :: occ(:)

      integer isps
      integer n
      integer i, iword, imask
      occ = 0
      n = 0
      do isps = 1,nhsps(abs(ith))
         iword = (isps-1)/(max_bit_word) + 1
         i = isps - (iword -1)*(max_bit_word)
         imask = ishft(1,i-1)
         if(iand(imask,hsd(iword)) /= 0)then
            n = n + 1
            occ(n) = isps
            if(n==nh)return
         endif
      enddo
      if(n < nh)then
        print*,' problem converting ',n,nh
        stop
      endif
      return
      end subroutine convert_haiku_OLD

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!
!  subroutine convert_haiku NEW
!  converts a haiku  (binary word) to an array ASSUMING fixed nh
!  NEW VERSION March 2013 to speed up usage for certain cases
!  reduces division which is slow
!  note this has a different definition of occ from convert_haiku_array
!
!  INPUT:
!   ith = species/handness
!   nh = # of paricles
!   hsd(:) = pointer to binary words that make up the haiku
!
!  OUTPUT:
!    occ(:)  array = 0. or 1 designating occupied s.p. states
!
!
      subroutine convert_haiku(ith,nh,hsd,occ)
      use bitstuff
      use haiku_info
      implicit none

      integer :: ith
      integer :: nh
      integer,pointer :: hsd(:)
      integer :: occ(:)

      integer isps
      integer n
      integer i, iword, imask,ilast,jhsd
      occ = 0
      n = 0

      isps = 0
      do iword = 1,Nword( abs(ith) )
         jhsd=hsd(iword)
!............. SET UP SIZE OF INNER DO LOOP .....
         if(iword == Nword( abs(ith) ))then
            ilast = nhsps(abs(ith))-(Nword(abs(ith)) -1)*(max_bit_word)
         else
            ilast = max_bit_word
         end if
         imask = 1
         do i = 1,ilast
            isps = isps + 1
            if(iand(imask,jhsd) /= 0)then
               n = n + 1
               occ(n) = isps
            endif
            imask = ishft(imask,1)  ! shift by 1
         end do
         if(n==nh)return

      end do ! iword

      if(n < nh)then
        print*,' problem converting ',n,nh
        stop
      endif
      return
      end subroutine convert_haiku

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  function nh_of_occ0  (mostly obsolete)
!
!  computes the nh of a haiku, using an occupation array
!  only used for testing, to verify that an occupation array has nh particles
!
!  input: ith = species/handness
!      occ(:) = occupation of s.p. states
!
!
      function nh_of_occ(ith,occ)
      use spsections
      use haiku_info

      implicit none
      integer :: ith
      integer :: occ(nhsps(abs(ith)))
      integer nh_of_occ
      integer isps

      nh_of_occ = 0
      do isps = 1,nhsps(abs(ith))
         nh_of_occ = nh_of_occ + occ(isps)
      enddo
      return
      end function nh_of_occ
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  function wh_of_occ0  (mostly obsolete)
!
!  computes the w of a haiku, using an occupation array
!  the difference from wh_of_occ is does not use nh to stop counting
!
!  input: ith = species/handness
!      occ(:) = occupation of s.p. states
!
!
      function wh_of_occ0(ith,occ)
      use spstate
      use haiku_info

      implicit none
      integer :: ith
      integer :: occ(nhsps(abs(ith)))
      integer wh_of_occ0
      integer isps

      wh_of_occ0 = 0
      do isps = 1,nhsps(abs(ith))
         wh_of_occ0 = wh_of_occ0 + occ(isps)*hspsqn(ith,isps)%w
      enddo
      return
      end function wh_of_occ0

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  function jzh_of_occ0  (mostly obsolete)
!
!  computes the jz of a haiku, using an occupation array
!  the difference from jzh_of_occ is does not use nh to stop counting
!
!  input: ith = species/handness
!      occ(:) = occupation of s.p. states
!
!
      function jzh_of_occ0(ith,occ)
      use spstate
      use haiku_info

      implicit none
      integer :: ith
      integer :: occ(nhsps(abs(ith)))
      integer jzh_of_occ0
      integer isps

      jzh_of_occ0 = 0
      do isps = 1,nhsps(abs(ith))
         jzh_of_occ0 = jzh_of_occ0 + occ(isps)*hspsqn(ith,isps)%m
      enddo
      return
      end function jzh_of_occ0

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  function wh_of_occ
!
!  computes the wh of a haiku, using an occupation array
!
!  input: ith = species/handness
!         nh = # of particles
!      occ(:) = occupation of s.p. states
!
!
      function wh_of_occ(ith,nh,occ)
      use spstate
      use haiku_info

      implicit none
      integer :: ith
      integer :: nh
      integer :: occ(nhsps(abs(ith)))
      integer wh_of_occ
      integer isps

      wh_of_occ = 0
      do isps = 1,nh
         wh_of_occ = wh_of_occ + hspsqn(ith,occ(isps))%w
      enddo
      return
      end function wh_of_occ

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  function jzh_of_occ
!
!  computes the jz of a haiku, using an occupation array
!
!  input: ith = species/handness
!         nh = # of particles
!      occ(:) = occupation of s.p. states
!
!
      function jzh_of_occ(ith,nh,occ)
      use spstate
      use haiku_info

      implicit none
      integer :: ith
      integer :: nh
      integer :: occ(nhsps(abs(ith)))
      integer jzh_of_occ
      integer isps

      jzh_of_occ = 0
      do isps = 1,nh
         jzh_of_occ = jzh_of_occ + hspsqn(ith,occ(isps))%m
      enddo
      return
      end function jzh_of_occ

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  function par_of_occ
!
!  computes the parity (1 + +, 2 = -) of a haiku, using an occupation array
!
!  input: ith = species/handness
!         nh = # of particles
!      occ(:) = occupation of s.p. states
!
!
      function par_of_occ(ith,nh,occ)

      use sporbit
      use spstate
      use haiku_info

      implicit none
      integer :: ith
      integer :: nh
      integer :: occ(nhsps(abs(ith)))
      integer par_of_occ
      integer isps

      par_of_occ = 1
      if(allsameparity)return
      if(nh == 0)return

      do isps = 1,nh
         par_of_occ = par_of_occ*hspsqn(ith,occ(isps))%par
      enddo
      par_of_occ = (3-par_of_occ)/2
      if(par_of_occ < 1 .or. par_of_occ > 2)then
        print*,' problem with parity ',par_of_occ,nh
        stop
      endif
      return
  end       function par_of_occ

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!      subroutine test_haiku_qns
!
!   test the quantum numbers jz, wh, parh for a list of haikus
!   assumes the haikus have already been segregated by nh, w
!
!  INPUT:
!    ith  = species/handness
!    nh = # of particles, iwh = which W value
!    occ(:)  an array used repeatedly
!
!  SUBROUTINES CALLED:
!    convert_haiku : converts a haiku (hsd) to occupation array (occ)
!
!  FUNCTIONS CALLED:
!    jzh_of_occ, par_of_occ
!
      subroutine test_haiku_qns(ith,nh,iwh,occ)

      use haiku_info
      implicit none

      integer ith
      integer nh
      integer iwh
      integer :: occ(:)
!............................................
      integer nhaiku
      integer ih
      integer,pointer :: hsd(:)
!      integer jzh_of_occ,par_of_occ

      nhaiku = haiku(ith)%of_nh(nh)%of_w(iwh)%nhaiku

      do ih = 1,nhaiku
         hsd => haiku(ith)%of_nh(nh)%of_w(iwh)%hsd(:,ih)
         call convert_haiku(ith,nh,hsd,occ)
         if(haiku(ith)%of_nh(nh)%of_w(iwh)%jz(ih) /=  jzh_of_occ(ith,nh,occ))then

             print*,' problem with sorted jz '
             print*,ith,nh,iwh
             print*,ih
             print*,haiku(ith)%of_nh(nh)%of_w(iwh)%jz(ih)
             print*,jzh_of_occ(ith,nh,occ)
             print*,' '
             print*,nhaiku
             print*,haiku(ith)%of_nh(nh)%of_w(iwh)%jz(:)

            stop
         endif
         if(haiku(ith)%of_nh(nh)%of_w(iwh)%par(ih) /= par_of_occ(ith,nh,occ))then

            print*,' problem with sorted par '
             print*,ith,nh,iwh
             print*,ih
            stop
         endif
      enddo
      return
      end subroutine test_haiku_qns

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

end module haikus

max_stars_repo_name: noctildon/pyBigstick

max_issues_repo_name: noctildon/pyBigstick

max_forks_repo_name: noctildon/pyBigstick


text:
Ednas Playground, for All Dogs Big and Small was created and opened in 2006, thanks due to donation by the Gorland family.  Located just outside the CCAH building of the VMTH, Ednas Playground is a very small Dog Parks dog park and sitting area.  An outside sitting area has several benches on top of carefully laid bricks, with some of them dedicated to certain donaters and even lain as memorials for loved pets.  Unfortunately, its also located near a large field that often contains several horses and a donkey, and a corresponding smell is noticable with a light breeze.

The main focus is of course the agility park!  Agility is a common dog sport seen both on TV and showcased at notable UC Davis Events Board events such as Picnic Day. Theres a pause table, weave poles, a couple of jumps, and contact obstacles such as a dog walk, a seesaw, an Aframe, and a tunnel.

Also see Dogs for more information and on other local parks

20090301 14:12:39 nbsp It may be small but its the only one thats free in the greater Sacramento area.   Anyhow, my dogs love the place and its good place to work out a high energy dog. Users/atwong

20110117 19:48:07 nbsp Oh, this is free? I always meant to ask, and now I know!  Users/ChristyMarsden

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C    @(#)ldptimdc.f	20.7 1/15/98
C****************************************************************
C
C     File: ldptimdc.f
C
C     Purpose: Routine to load PTI multi-terminal d-c data from raw
C              data file
C
c     Return code:  n = 0 : Success
c                   n = 1 : Error
c
C     Author: Walt Powell  Date: 21 May 1996
C     Called by: load_pti.f
C
C****************************************************************
      integer function ldptimdc (xbuf, tempfile, options, numnver,
     &                           numrec, error)
      integer tempfile, numrec, options(*), numver, error
      character *(*) xbuf

      include 'ipfinc/parametr.inc'

      include 'ipfinc/prt.inc'
      include 'ipfinc/pti_data.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/branch.inc'
      include 'ipfinc/lfiles.inc'
      include 'ipfinc/wsccbase.inc'

      common /old_branch/ old_k1, old_k2, old_sect, old_type, old_id
      integer old_k1, old_k2, old_sect, old_type
      character old_id * 1

      character word(50)*10, busname*8, type*1, code*10, bus1*8
      logical first, found
      integer find_bus, ftn_atoi, bptr1, bptr2, add_bus,
     &        array(6,12), count, ptr, gtptinam, status

      ldptimdc = 0
      error = 0
      numrec = 1
      nconv = 0
      ndcbs = 0
      ndcln = 0
c
c     Process first record.
c
      last = index (xbuf, '/') - 1
      if (last .le. 0) last = lastch(xbuf)
      call ptiscan (xbuf(1:last), word, nwrd, ',',  ' ')

      do i = nwrd+1, 8
        word(i) = ' '
      enddo

      nptidc = ftn_atoi(word(1))
      if (nptidc .eq. 0) then
        ldptimdc = 1
      else
        nconv = ftn_atoi(word(2))
        ndcbs = ftn_atof(word(3))
        ndcln = ftn_atof(word(4))
        mdc = ftn_atof(word(5))
        nvconv = ftn_atoi(word(6))
        if (nconv .eq. 0) then
          write (errbuf(1), 100) nptidc
  100     format (' D-C system (', i2, ') has no specified number of a-c
     & commutator buses.')
          errbuf(2) = ' (' // xbuf(1:60) // ')'
          call prterx ('W',2)
          error = 1
        endif
        if (ndcbs .eq. 0) then
          write (errbuf(1), 110) nptidc
  110     format (' D-C system (', i2, ') has no specified number of d-c
     & converter buses.')
          errbuf(2) = ' (' // xbuf(1:60) // ')'
          call prterx ('W',2)
          error = 1
        endif
        if (ndcln .eq. 0) then
          write (errbuf(1), 120) nptidc
  120     format (' D-C system (', i2, ') has no specified number of d-c
     & lines.')
          errbuf(2) = ' (' // xbuf(1:60) // ')'
          call prterx ('W',2)
          error = 1
        endif
        if (nvconv .eq. 0) then
          write (errbuf(1), 130) nptidc
  130     format (' D-C system (', i2, ') has no specified voltage-contr
     &olled d-c bus.')
          errbuf(2) = ' (' // xbuf(1:60) // ')'
          call prterx ('W',2)
          error = 1
        endif
        if (error .ne. 0) go to 304
c
c       Process converter a-c records.
c
        pwr_rect = 0.0
        pwr_invt = 0.0
        kvc = 0
        do ix = 1, nconv
          numrec = numrec + 1
          read (tempfile, fmt='(a)', end=910) xbuf
          last = index (xbuf, '/') - 1
          if (last .le. 0) last = lastch(xbuf)
          call ptiscan (xbuf(1:last), word, nwrd, ',',  ' ')

          do i = nwrd+1, 16
            word(i) = ' '
          enddo
c
c         Set defaults
c
          if (word(10) .eq. ' ') word(10) = '1.50'
          if (word(11) .eq. ' ') word(11) = '0.51'
          if (word(12) .eq. ' ') word(12) = '0.00625'
c
c	  Converter bus
c
          status1 = gtptinam (word(1), npti1, kp1, isgn1, bus1, basekv1)
          if (bus1 .eq. ' ') then
            write (errbuf(1), 160) nptidc, npti1
  160       format (' D-C system (', i2, ') Converter bus (', i6,
     &        ') is not in system.')
            call prterx ('W',1)
            error = 1
            go to 302
          endif
          array(1,ix) = npti1
          array(2,ix) = kp1
          k1 = find_bus (bus1, basekv1)
          if (k1 .le. 0) then
            write (errbuf(1), 170) nptidc, npti1, bus1, basekv1
  170       format (' D-C system (', i2, ') Converter bus (', i6,
     &        1x, a8, f7.1, ') is not in system.')
            call prterx ('W',1)
            error = 1
            go to 302
          endif
          array(3,ix) = k1
          kbsdta(1,k1) = 12          ! Converter bus changed to type BM
c
c         Find converter commutating bus.
c
          ptr = kbsdta(16,k1)
          k1c = 0
          found = .false.
          do while (ptr .gt. 0 .and. .not. found)
            if (brtype(ptr) .eq. 3 .or. brtype(ptr) .eq. 4 .or.
     &          brtype(ptr) .eq. 5) then
              found = .true.
              k1c = ky(ptr)
            else
              ptr = brnch_nxt(ptr)
            endif
          enddo
          if (k1c .le. 0) then
            write (errbuf(1), 180) nptidc, npti1, bus1, basekv1
  180       format (' D-C system (', i2,
     &        ') Commutating bus for Converter (', i6, 1x, a8, f7.1,
     &        ') is not in system.')
            call prterx ('W',1)
            error = 1
            go to 302
          endif
          array(4,ix) = k1c
          br = ftn_atof (word(2))
          if (br .le. 0.0) then
            write (errbuf(1), 220) nptidc, npti1, bus(k1), base(k1)
  220       format (' D-C system (', i2,
     &          ') Converter bus (', i6, 1x, a8, f7.1,
     &          ') has zero bridges')
            call prterx ('W',1)
            error = 1
            go to 302
          endif
          busdta(3,k1) = br
          busdta(4,k1) = 0.0       ! Smoothing reactor

          alpha = ftn_atof (word(3))
          if (alpha .le. 0.0) then
            write (errbuf(1), 230) nptidc, npti1, bus(k1), base(k1)
  230       format (' D-C system (', i2,
     &        ') Converter bus (', i6, 1x, a8, f7.1,
     &        ') has zero Alpha_nom')
            call prterx ('W',1)
            error = 1
            go to 302
          endif
          busdta(6,k1) = 100.0
          busdta(10,k1) = alpha
          alphamin = ftn_atof (word(4))
          if (alphamin .le. 0.0) then
            write (errbuf(1), 240) nptidc, npti1, bus(k1), base(k1)
  240       format (' D-C system (', i2,
     &        ') Converter bus (', i6, 1x, a8, f7.1,
     &        ') has zero Alpha_min')
            call prterx ('W',1)
            error = 1
            go to 302
          endif
          busdta(5,k1) = alphamin
          busdta(11,k1) = alphamin

          busdta(7,k1) = 0.0       ! Voltage drop / bridge
          busdta(8,k1) = 0.0       ! Rating
          kbsdta(9,k1) = k1c

          sched = ftn_atof (word(13))
          if (npti1 .eq. nvconv) then
            if (sched .le. 0.0) then
              write (errbuf(1), 242) nptidc, npti1, bus(k1), base(k1),
     &          sched
  242         format (' D-C system (', i2,
     &        ') converter bus (', i6, 1x, a8, f7.1,
     &        ') has illegal scheduled voltage (', f8.1, ')')
              call prterx ('W',1)
              error = 1
            endif
            kvc = k1
            call putchr(1, ' ', kbsdta(12,k1))
            busdta(13,k1) = 0.0
            busdta(14,k1) = sched
          else
            pwr_rect = pwr_rect + amax1 (0.0, sched)
            pwr_invt = pwr_invt + amin1 (0.0, sched)
            busdta(13,k1) = sched
            busdta(14,k1) = 0.0
            if (sched .gt. 0.0) then
              call putchr(1, 'R', kbsdta(12,k1))
            else if (alpha .eq. alphamin) then
              call putchr(1, 'M', kbsdta(12,k1))
            else
              call putchr(1, 'I', kbsdta(12,k1))
            endif
          endif
c
c         The next test is trivial, since
c         baseratio = base(k1) / base(k1c) and both are known.
c
          baseratio = ftn_atof (word(8))
          if (baseratio .le. 0.0) then
            write (errbuf(1), 250) nptidc, npti1, bus(k1), base(k1)
  250       format (' D-C system (', i2,
     &        ') Converter bus (', i6, 1x, a8, f7.1,
     &        ') has zero base ratio')
            call prterx ('W',1)
            error = 1
            go to 302
          endif
          tol = (base(k1) / base(k1c)) / baseratio - 1.0
          if (abs(tol) .gt. 0.15) then
            write (errbuf(1), 260) nptidc, npti1, bus(k1), base(k1),
     &        baseratio
  260       format (' D-C system (', i2,
     &        ') Converter bus (', i6, 1x, a8, f7.1,
     &        ') has an unacceptable base ratio (', f8.4, ')')
            call prterx ('W',1)
            error = 1
            go to 302
          endif
          npolar = ftn_atoi (word(16))
          if (npolar .lt. 0) then
            write (errbuf(1), 270) nptidc
  270       format (' D-C system (', i2, ') has specified negative polar
     &ity. BPA model accepts only positive polarity')
            errbuf(2) = ' (' // xbuf(1:60) // ')'
            call prterx ('W',2)
            error = 1
            go to 302
          endif
c
c         Replace existing converter-commutating branch with an LTC
c
          ptr = kbsdta(16,k1)
          found = .false.
          do while (ptr .gt. 0 .and. .not. found)
            if (ky(ptr) .eq. k1c .and. brtype(ptr) .eq. 4) then
              found = .true.
            else
              ptr = brnch_nxt(ptr)
            endif
          enddo
          if (brtype(ptr) .ne. 4 .or. ky(ptr) .ne. k1c) then
c
c           Add LTC commuating record and link up branch

            if (ltot+1 .ge. MAXBRN) then
              write (errbuf(1), 280) MAXBRN
  280         format ('More than ',i5, ' branch records.')
              write (errbuf(2), 290) nptidc, bus(k1), base(k1),
     &          bus(k1c), base(k1c)
  290         format (' D-C system (', i2,
     &          ') added record (R ', a8, f7.1, 1x, a8, f7.1, ')')
              call prterx ('W',2)
              ltot = 1
              error = 1
              go to 302
            endif
            ltot = ltot + 1
            kbrnch(1,ltot) = 4
            call putchr(3, '   ', kbrnch(3,ltot))
            kbrnch(4,ltot) = k1
            brnch(6,ltot) = ftn_atof(word(10)) * base(k1c) ! Tmax (kV)
            brnch(7,ltot) = ftn_atof(word(11)) * base(k1c) ! Tmin (kV)
            taps = ftn_atof(word(12))
            if (taps .gt. 0.0) then
              numtaps = (brnch(6,ltot) - brnch(7,ltot))
     &                / (taps * base(k1c)) + 1.5
              if (numtaps .gt. 50) numtaps = 0
              brnch(8,ltot) = numtaps
            else
              brnch(8,ltot) = 0.0
            endif

            rateln(1,ltot) = 0.0        ! Rate A
            rateln(2,ltot) = 0.0        ! Rate B
            rateln(3,ltot) = 0.0        ! Rate C

            call lkbrdata (ltot, k1, k1c, bptr1, bptr2, error)
            kx(bptr1) = k1
            ky(bptr1) = k1c
            brid(bptr1) = ' '
            brsect(bptr1) = 0
            brtype(bptr1) = kbrnch(1,ltot)

            kx(bptr2) = k1c
            ky(bptr2) = k1
            brid(bptr2) = ' '
            brsect(bptr2) = 0
            brtype(bptr2) = kbrnch(1,ltot)
          endif
c
c         Replace existing bus tie with commutating transformer
c
          ptr = kbsdta(16,k1)
          found = .false.
          do while (ptr .gt. 0 .and. .not. found)
            if (ky(ptr) .eq. k1c .and.
     &         (brtype(ptr) .eq. 3 .or. brtype(ptr) .eq. 5)) then
              found = .true.
            else
              ptr = brnch_nxt(ptr)
            endif
          enddo
          if (ky(ptr) .ne. k1c) then
            write (errbuf(1), 300) nptidc, npti1, bus(k1c), base(k1c),
     &        bus(k1), base(k1)
  300       format (' D-C system (', i2,
     &        ') Converter commutating bus (', i6, 1x, a8, f7.1,
     &        ') is not adjacent to Converter.')
            call prterx ('W',1)
            error = 1
            go to 302
          endif
          nbr = iabs (brnch_ptr(ptr))
          rt = ftn_atof (word(5))
          xt = ftn_atof (word(6))
          cx = base(k1)**2 / bmva * br

          brtype(ptr) = 5
          kbrnch(1,nbr) = brtype(ptr)
          brnch(5,nbr) = rt / cx
          brnch(6,nbr) = xt / cx
          if (brnch_ptr(ptr) .gt. 0) then
            brnch(9,nbr) = base(k1)
            brnch(10,nbr) = base(k1c) * ftn_atof (word(9))
          else
            brnch(9,nbr) = base(k1c) * ftn_atof (word(9))
            brnch(10,nbr) = base(k1)
          endif
          brnch(16,nbr) = busdta(3,k1)
c
c         Redefine brtype() for transpose
c
          ptr = kbsdta(16,k1c)
          found = .false.
          do while (ptr .gt. 0 .and. .not. found)
            if (ky(ptr) .eq. k1 .and.
     &         (brtype(ptr) .eq. 3 .or. brtype(ptr) .eq. 5)) then
              found = .true.
            else
              ptr = brnch_nxt(ptr)
            endif
          enddo
          if (ptr .gt. 0) brtype(ptr) = kbrnch(1,nbr)

  302     continue
c
c         End of converter loop
        enddo
  304   continue
c
c       Determine code ("I" or "R") of voltage-controlled converter
c
        if (kvc .eq. 0) then
          write (errbuf(1), 306) nptidc, nvconv
  306     format (' D-C system (', i2,
     &      ') specified voltage controlled bus (', i6,
     &      ') has no assigned schedule')
          call prterx ('W',1)
          error = 1
        else
          if (0.90 * pwr_rect + pwr_invt .lt. 0.0) then
            call putchr(1, 'R', kbsdta(12,kvc))
          else
            call putchr(1, 'I', kbsdta(12,kvc))
          endif
        endif
c
c       Process converter d-c records.
c
        count = nconv
        do ix = 1, ndcbs
          numrec = numrec + 1
          read (tempfile, fmt='(a)', end=910) xbuf
          last = index (xbuf, '/') - 1
          if (last .le. 0) last = lastch(xbuf)
          call ptiscan (xbuf(1:last), word, nwrd, ',',  ' ')

          do i = nwrd+1, 7
            word(i) = ' '
          enddo
c
c         Set defaults
c
          if (word(3) .eq. ' ') word(3) = '1'
          if (word(4) .eq. ' ') word(4) = '1'
          if (word(5) .eq. ' ') then
            write (word(5), 310) ix
  310       format ('D-C # ', i2)
          endif

          nptiln = ftn_atoi(word(1))
          status = gtptinam (word(2), npti1, kp1, isgn1, bus1, basekv1)
          if (npti1 .eq. 0 .and. bus1 .eq. ' ') then
c
c           If a converter bus is not specified, it has no a-c entity
c           and a pseudo bus must be created.
c
            busname = word(5)
            basekv = base(array(1,1))
            indx = add_bus (busname, basekv, ntot+1)
C
C           This hashing function checks for duplicates, array
C           overflow, and automatically inserts the new bus in
C           position ntot+1.
c
            if (indx .lt. 0) then
              write (errbuf(1), 320) busname, basekv
  320         format (' Cannot generate a unique d-c bus name (',
     &          a8, f7.1, ')')
              errbuf(2) = ' (' // xbuf(1:60) // ')'
              call prterx ('W',2)
              error = 1
              go to 382
            else if (indx .eq. 0) then
              write (errbuf(1), 330) MAXBUS
  330         format (' More than ',i4,' buses in base system ')
              errbuf(2) = ' (' // xbuf(1:60) // ')'
              call prterx ('W',2)
              error = 1
              go to 382
            endif
            ntot = indx
            count = count + 1
            array(1,count) = 0
            array(2,count) = 0
            array(3,count) = ntot
            array(4,count) = 0
            array(5,count) = nptiln
            if (wsccflag) then
               wsccbase(ntot) = code (basekv, 4, 0)
            endif
          else
            if (bus1 .eq. ' ') then
              write (errbuf(1), 340) nptiln, npti1
  340         format (' D-C system (', i2, ') Converter bus (', i6,
     &        ') is not raw data file.')
              call prterx ('W',1)
              error = 1
              go to 382
            endif
            k1 = find_bus (bus1, basekv1)
            if (k1 .le. 0) then
              write (errbuf(1), 350) nptiln, npti1, bus1, basekv1
  350         format (' D-C system (', i2, ') converter bus (',
     &          i6, 1x, a8, f7.1, ') is not in system.')
              call prterx ('W',1)
              error = 1
              go to 382
            endif
            found = .false.
            i = 1
            do while (i .le. count .and. .not. found)
              if (array(3,i) .eq. k1) then
                found = .true.
                array(5,i) = nptiln
              else
                i = i + 1
              endif
            enddo
            if (.not. found) then
              write (errbuf(1), 360) nptiln, npti1, bus(k1), base(k1)
  360         format (' D-C system (', i2, ') converter bus (', i6,
     &          1x, a8, f7.1, ') has no d-c entity.')
              call prterx ('W',1)
              error = 1
              go to 382
            endif
          endif
          idc2 = ftn_atoi (word(6))
          if (idc2 .gt. 0) then
            write (errbuf(1), 370) nptiln, npti1, bus(k1), base(k1)
  370       format (' D-C system  (', i2, ') converter bus (', i6,
     &        1x, a8, f7.1, ') is series-connected. BPA model accepts on
     &ly parallel-connected.')
            call prterx ('W',1)
            error = 1
            go to 382
          endif
          rgnd = ftn_atof (word(7))
          if (rgnd .gt. 0) then
            write (errbuf(1), 380) nptiln, npti1, bus(k1), base(k1)
  380       format (' D-C system (', i2, ') Converter bus (', i6,
     &        1x, a8, f7.1, ') will ignore non-zero resistance to ground
     &')
            call prterx ('W',1)
            error = 1
          endif
  382     continue
        enddo
c
c       Process d-c line records.
c
        do ix = 1, ndcln
          numrec = numrec + 1
          read (tempfile, fmt='(a)', end=910) xbuf
          last = index (xbuf, '/') - 1
          if (last .le. 0) last = lastch(xbuf)
          call ptiscan (xbuf(1:last), word, nwrd, ',',  ' ')
c
c         Squeeze and blank fill comma fields (Note: first comma
c         is false)
c
          first = .true.
          i = 1
          do while (i .le. nwrd)
            if (word(i) .eq. ',') then
              if (first) then
                do j = i, nwrd-1
                  word(j) = word(j+1)
                enddo
                first = (word(i) .ne. ',')
                if (word(i) .eq. ',') word(i) = ' '
                nwrd = nwrd - 1
              else
                word(i) = ' '
              endif
            else
              first = .true.
            endif
            i = i + 1
          enddo
          do i = nwrd+1, 5
            word(i) = ' '
          enddo
c
c         Set defaults
c
          if (word(3) .eq. ' ') word(3) = '1'
          if (word(5) .eq. ' ') word(4) = '0.0'

          npti1 = ftn_atoi(word(1))
          isgn1 = isign(1, npti1)
          npti2 = ftn_atoi(word(2))
          isgn2 = isign(1, npti2)
          npti2 = iabs (npti2)

          found = .false.
          i = 1
          do while (i .le. count .and. .not. found)
            if (array(5,i) .eq. npti1) then
              found = .true.
              k1 = array(3,i)
            else
              i = i + 1
            endif
          enddo
          if (.not. found) then
            write (errbuf(1), 390) nptidc, npti1, npti2
  390       format (' D-C system (', i2, ') line (', i2, ' - ', i6,
     &         ') is missing terminal bus 1 d-c bus entity.')
            call prterx ('W',1)
            error = 1
            go to 460
          endif
          found = .false.
          i = 1
          do while (i .le. count .and. .not. found)
            if (array(5,i) .eq. npti2) then
              found = .true.
              k2 = array(3,i)
            else
              i = i + 1
            endif
          enddo
          if (.not. found) then
            write (errbuf(1), 400) nptidc, npti1, npti2
  400       format (' D-C system (', i2, ') line (', i2, ' - ', i6,
     &         ') is missing terminal bus 2 d-c bus entity.')
            call prterx ('W',1)
            error = 1
            go to 460
          endif
          if (kbsdta(1,k1) .ne. 12) then
            call typno (type, kbsdta(1,k1))
            write (errbuf(1), 410) nptidc, npti1, bus(k1), base(k1),
     &        npti2, bus(k2), base(k2), type
  410       format (' D-C system (', i2, ') line (', i2, 1x, a8, f7.1,
     &         ' - ', i6, 1x, a8, f7.1,
     &         ') terminal bus 1 is incorrect type (', a, ')')
            call prterx ('W',1)
            error = 1
            go to 460
          endif
          if (kbsdta(1,k2) .ne. 12) then
            call typno (type, kbsdta(1,k2))
            write (errbuf(1), 420) nptidc, npti1, bus(k1), base(k1),
     &        npti2, bus(k2), base(k2), type
  420       format (' D-C system (', i2, ') line (', i2, 1x, a8, f7.1,
     &         ' - ', i6, 1x, a8, f7.1,
     &         ') terminal bus 2 is incorrect type (', a, ')')
            call prterx ('W',1)
            error = 1
            go to 460
          endif
          rdc = ftn_atof (word(4))
          if (rdc .le. 0.0) then
            write (errbuf(1), 430) nptidc, npti1, bus(k1), base(k1),
     &        npti2, bus(k2), base(k2)
  430       format (' D-C system (', i2, ') line (', i2, 1x, a8, f7.1,
     &         ' - ', i6, 1x, a8, f7.1,
     &         ') has zero d-c resistance ')
            call prterx ('W',1)
            error = 1
            go to 460
          endif
          if (ltot+1 .ge. MAXBRN) then
            write (errbuf(1), 440) MAXBRN
  440       format ('More than ',i5, ' branch records. ')
            write (errbuf(2), 450) nptidc, npti1, bus(k1), base(k1),
     &        npti2, bus(k2), base(k2)
  450       format (' Adding d-c system (', i2, ') line (', i2, 1x,
     &         a8, f7.1, ' - ', i6, 1x, a8, f7.1)
            call prterx ('W',2)
            ltot = 1
            error = 1
            go to 460
          endif
          ltot = ltot + 1
          kbrnch(1,ltot) = 2
          call putchr(3, '   ', kbrnch(3,ltot))
          brnch(4,ltot) = 0.0              ! Rating
          brnch(5,ltot) = rdc              ! Rdc
          do i = 6, 18
            brnch(i,ltot) = 0.0
          enddo
          if (isgn1 .lt. 0 .or. isgn2 .gt. 0) then
            kbrnch(15,ltot) = 1            ! Metered end
          else
            kbrnch(15,ltot) = 2            ! Controlled end
          endif

          call lkbrdata (ltot, k1, k2, bptr1, bptr2, error)
          kx(bptr1) = k1
          ky(bptr1) = k2
          brid(bptr1) = word(3)
          brsect(bptr1) = 0
          brtype(bptr1) = kbrnch(1,ltot)

          kx(bptr2) = k2
          ky(bptr2) = k1
          brid(bptr2) = word(3)
          brsect(bptr2) = 0
          brtype(bptr2) = kbrnch(1,ltot)

  460     continue
        enddo

  900   continue
c
c       If any errors are encountered, dummy read any remaining records
c       in multi-terminal d-c set.
c
        do while (numrec .lt. 1+nconv+ndcbs+ndcln)
          read (tempfile, fmt='(a)', end=910) xbuf
          numrec = numrec + 1
        enddo

        go to 930

  910   write (errbuf(1), 920)
  920   format ('E-O-F encountered processing multi-terminal d-c data re
     &cords in raw data file ')
        call prterx ('W', 1)
        error = 1

  930   continue
      endif
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


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
!

module mod
logical expect(14), rslt(14)
type :: objects(k1,k2,l1,l2)
integer, kind :: k2 = selected_char_kind("ASCII")
integer, kind :: k1 = selected_int_kind(4)
integer,len :: l1 = 10
integer,len :: l2 = 10
integer(k1) :: p(l1,l2)
end type
contains
subroutine foo(n)
integer n
type(objects) :: x

rslt(5) = x%k1 .eq. selected_int_kind(4)
rslt(6) = x%k2 .eq. selected_char_kind("ASCII")
rslt(7) = (x%l1 .eq. n)
rslt(8) = size(x%p,dim=1) .eq. x%l1
rslt(9) = size(x%p,dim=2) .eq. x%l1

k=0
do i=1,size(x%p,dim=1)
do j=1, size(x%p,dim=2)
  x%p(i,j) = k
  k = k + 1
enddo
enddo

k=0
rslt(10) = .true.
do i=1,size(x%p,dim=1)
do j=1, size(x%p,dim=2)
  if (x%p(i,j) .ne. k) then
      print *, 'element ',i,',',j,' equals ',x%p(i,j),' but expecting ',k
      rslt(10) = .false.
  endif
  k = k + 1
enddo
enddo
rslt(11) = i .ne. size(x%p,dim=1)
rslt(12) = j .ne. size(x%p,dim=2)
rslt(13) = (x%l2 .eq. n)
rslt(14) = kind(x%p) .eq. x%k1

end subroutine

end module

program p
use mod
integer y
type(objects(1,20,30)) :: x

rslt = .false.
rslt(1) = x%k1 .eq. 1
rslt(2) = x%k2 .eq. 20
rslt(3) = x%l1 .eq. 30
rslt(4) = size(x%p, dim=1) .eq. x%l1

call foo(10)

expect = .true.
call check(rslt,expect,14)

end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text: One of the local radio stations. KSFM plays mostly hiphop/rap/r&b top 40. They also have a wiki:WikiPedia:KSFM wikipedia page.


text:
Jump to Timeline #Navigation Navigation

## If you are adding something, go ahead and remove the notice below!
  Sword and Sandals founded in December

## ################## Leave this navigation here:


text:
****************************************************************************
*                                                                          *
*     GENERAL SUBROUTINES FOR TROPUS                                       *
*     (C)Steve King, ISIS Facility, Rutherford Appleton Laboratory, UK     *
*                                                                          *
****************************************************************************

	SUBROUTINE LOAD_DATA(FNAME,DATA_CHANNELS,DATA_ARRAY)

	character*80	fname,prompt(2)
	integer*4		data_channels
	real*4		data_array(4096)

      prompt(1)='ERROR: Error opening Tropus input file: FATAL'
      prompt(2)='ERROR: Error reading Tropus input file: FATAL'

	data_channels=1
	open(12,file=fname,form='formatted',status='old',err=8001)
10	   read(12,*,err=8002,end=1000)data_array(data_channels)
	   data_channels=data_channels+1
	   goto 10
1000	close(12)
	data_channels=data_channels-1
	goto 9999

8001  call showprompt(prompt(1))
      stop

8002  call showprompt(prompt(2))
      stop

9999	return
	end

	SUBROUTINE SAVE_DATA(FNAME,DATA_CHANNELS,X_DATA_ARRAY,Y_DATA_ARRAY)

	character*80	fname,prompt(2)
	integer*4		data_channels
	real*4		x_data_array(4096),y_data_array(4096)

      prompt(1)='ERROR: Error opening Tropus output file: FATAL'
      prompt(2)='ERROR: Error writing Tropus output file: FATAL'

1000	format(f12.5,' ,',e12.5)

	open(12,file=fname,form='formatted',status='unknown',err=8001)
	   do i=1,data_channels
	      write(12,1000,err=8002)x_data_array(i),y_data_array(i)
	   end do
	close(12)
	goto 9999

8001  call showprompt(prompt(1))
      stop

8002  close(12)
      call showprompt(prompt(2))
      stop

9999	return
	end


text:
c     this is model YRH

c     Aug 9, 1999
c     derivs included
c

********************************************

      subroutine prepot
      integer g_p_
      parameter(g_p_=3)
      write(6,*)'Model YRH, prepot called'
      write(6,*)'Derivs included'
      write(6,*)'Max coupling = 0.2'
      call g_prepot11(g_p_)
      call g_prepot12(g_p_)
      call g_prepot22(g_p_)
      return
      end

********************************************

      subroutine pot(ro,e,de,nt,nsurf)
c      subroutine pot(ro,nt,nsurf)

      implicit double precision(a-h,o-z)
      dimension r(nt,3),De(3,nt),e(nt),ro(nt,3)
      dimension deo(3,nt)

      integer g_p_
      parameter (g_p_=3,ldg_r=3,ldg_e=3)
      dimension g_r(ldg_r,nt,3),g_e(ldg_e,nt)

c     permute the bonds
c     want r1 = YR
c          r2 = RH
c          r3 = YH

c     so arr #1 = Y + RH
c            #2 = R + YH
c            #3 = H + YR

c     in code:  r1 = RH
c               r2 = YH
c               r3 = YR

c     so...

      do i=1,nt
       r(i,1) = ro(i,2)
       r(i,2) = ro(i,3)
       r(i,3) = ro(i,1)
      enddo

c     seed matrix

        do i=1,nt
        do j=1,3
        do k=1,3
         if (j .eq. k) g_r(j,i,k)=1.0d0
         if (j .ne. k) g_r(j,i,k)=0.0d0
        enddo
        enddo
        enddo

      if (nsurf.eq.1) then
          call g_pot11(g_p_, r, g_r, ldg_r, e, Deo, ldg_e, nt)
      else if (nsurf.eq.2) then
          call g_pot12(g_p_, r, g_r, ldg_r, e, Deo, ldg_e, nt)
      else
          call g_pot22(g_p_, r, g_r, ldg_r, e, Deo, ldg_e, nt)
      endif

      do i=1,nt
       de(1,i) = Deo(3,i)
       de(2,i) = Deo(1,i)
       de(3,i) = Deo(2,i)
      enddo

      return
      end

*******************************************************************
*******************************************************************

C                           DISCLAIMER
C
C   This file was generated on 08/09/99 by the version of
C   ADIFOR compiled on June, 1998.
C
C   ADIFOR was prepared as an account of work sponsored by an
C   agency of the United States Government, Rice University, and
C   the University of Chicago.  NEITHER THE AUTHOR(S), THE UNITED
C   STATES GOVERNMENT NOR ANY AGENCY THEREOF, NOR RICE UNIVERSITY,
C   NOR THE UNIVERSITY OF CHICAGO, INCLUDING ANY OF THEIR EMPLOYEES
C   OR OFFICERS, MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES
C   ANY LEGAL LIABILITY OR RESPONSIBILITY FOR THE ACCURACY, COMPLETE-
C   NESS, OR USEFULNESS OF ANY INFORMATION OR PROCESS DISCLOSED, OR
C   REPRESENTS THAT ITS USE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS.
C
C*******************************************
C U11- lower diabat
C
      subroutine g_prepot11(g_p_)
C
        implicit double precision (a-h, o-z)
        dimension r(nt, 3), e(nt)
C
C     data is in bohr, and eV units                                    x
C
        save
        integer g_pmax_
        parameter (g_pmax_ = 3)
        integer g_i_, g_p_, ldg_r, ldg_e
        double precision d3_p, d7_b, d10_b, d7_v, d2_p, d5_b, d1_w, d10_
     *v, d2_v, d6_b
        double precision d1_p, d2_b, d3_b, d4_v, g_r1(g_pmax_), g_r(ldg_
     *r, nt, 3), g_r2(g_pmax_), g_r3(g_pmax_), g_d1_w(g_pmax_), g_s1(g_p
     *max_)
        double precision g_s2(g_pmax_), g_s3(g_pmax_), g_t1(g_pmax_), g_
     *t2(g_pmax_), g_t3(g_pmax_), g_coul1(g_pmax_), g_coul2(g_pmax_), g_
     *coul3(g_pmax_), g_exch1(g_pmax_), g_exch2(g_pmax_)
        double precision g_exch3(g_pmax_), g_w(g_pmax_), g_c1(g_pmax_),
     *g_rootw(g_pmax_), g_e(ldg_e, nt)
        integer g_ehfid
        save g_coul1, g_coul2, g_coul3, g_exch1, g_exch2, g_exch3, g_w,
     *g_c1, g_rootw
        save g_r1, g_r2, g_r3, g_d1_w, g_s1, g_s2, g_s3, g_t1, g_t2, g_t
     *3
        data de1_11 /3.9d0/, de2_11 /4.3d0/, de3_11 /0.4d0/
        data re1_11 /2.1d0/, re2_11 /2.1d0/, re3_11 /2.5d0/
        data beta1_11 /1.0d0/, beta2_11 /1.0d0/, beta3_11 /1.5d0/
        data ca_11 /1.5d0/, cb_11 /1.0d0/, cc_11 /0.15d0/
        data sp1_11 /0.9d0/, sp2_11 /0.9d0/, sp3_11 /0.2d0/
        data g_ehfid /0/
C
C
        if (g_p_ .gt. g_pmax_) then
          print *, 'Parameter g_p_ is greater than g_pmax_'
          stop
        endif
        cevau = 1.d0 / 27.2113961d0
C
        return
C
        entry g_pot11(g_p_, r, g_r, ldg_r, e, g_e, ldg_e, nt)
C
        do i = 1, nt
C
          do g_i_ = 1, g_p_
            g_r1(g_i_) = g_r(g_i_, i, 1)
          enddo
          r1 = r(i, 1)
C--------
          do g_i_ = 1, g_p_
            g_r2(g_i_) = g_r(g_i_, i, 2)
          enddo
          r2 = r(i, 2)
C--------
          do g_i_ = 1, g_p_
            g_r3(g_i_) = g_r(g_i_, i, 3)
          enddo
          r3 = r(i, 3)
C--------
C
C   SINGLETS (s1, s2, s3)
          do g_i_ = 1, g_p_
            g_d1_w(g_i_) = (-beta1_11) * g_r1(g_i_)
          enddo
          d1_w = (-beta1_11) * (r1 - re1_11)
          d2_v = exp(d1_w)
          d2_p =  d2_v
          d4_v = (1.d0 - d2_v) * (1.d0 - d2_v)
          d1_p = 2.0d0 * (1.d0 - d2_v)
          d6_b = (-(de1_11 * d1_p)) * d2_p
          do g_i_ = 1, g_p_
            g_s1(g_i_) = d6_b * g_d1_w(g_i_)
          enddo
          s1 = de1_11 * d4_v - de1_11
C--------
          do g_i_ = 1, g_p_
            g_d1_w(g_i_) = (-beta2_11) * g_r2(g_i_)
          enddo
          d1_w = (-beta2_11) * (r2 - re2_11)
          d2_v = exp(d1_w)
          d2_p =  d2_v
          d4_v = (1.d0 - d2_v) * (1.d0 - d2_v)
          d1_p = 2.0d0 * (1.d0 - d2_v)
          d6_b = (-(de2_11 * d1_p)) * d2_p
          do g_i_ = 1, g_p_
            g_s2(g_i_) = d6_b * g_d1_w(g_i_)
          enddo
          s2 = de2_11 * d4_v - de2_11
C--------
          do g_i_ = 1, g_p_
            g_d1_w(g_i_) = (-beta3_11) * g_r3(g_i_)
          enddo
          d1_w = (-beta3_11) * (r3 - re3_11)
          d2_v = exp(d1_w)
          d2_p =  d2_v
          d4_v = (1.d0 - d2_v) * (1.d0 - d2_v)
          d1_p = 2.0d0 * (1.d0 - d2_v)
          d6_b = (-(de3_11 * d1_p)) * d2_p
          do g_i_ = 1, g_p_
            g_s3(g_i_) = d6_b * g_d1_w(g_i_)
          enddo
          s3 = de3_11 * d4_v - de3_11
C--------
C
C   TRIPLETS (t1, t2, t3)
          do g_i_ = 1, g_p_
            g_d1_w(g_i_) = (-beta1_11) * g_r1(g_i_)
          enddo
          d1_w = (-beta1_11) * (r1 - re1_11)
          d2_v = exp(d1_w)
          d2_p =  d2_v
          d4_v = (1.d0 + d2_v) * (1.d0 + d2_v)
          d1_p = 2.0d0 * (1.d0 + d2_v)
          d7_b = 0.5d0 * de1_11 * d1_p * d2_p
          do g_i_ = 1, g_p_
            g_t1(g_i_) = d7_b * g_d1_w(g_i_)
          enddo
          t1 = 0.5d0 * (de1_11 * d4_v - de1_11)
C--------
          do g_i_ = 1, g_p_
            g_d1_w(g_i_) = (-beta2_11) * g_r2(g_i_)
          enddo
          d1_w = (-beta2_11) * (r2 - re2_11)
          d2_v = exp(d1_w)
          d2_p =  d2_v
          d4_v = (1.d0 + d2_v) * (1.d0 + d2_v)
          d1_p = 2.0d0 * (1.d0 + d2_v)
          d7_b = 0.5d0 * de2_11 * d1_p * d2_p
          do g_i_ = 1, g_p_
            g_t2(g_i_) = d7_b * g_d1_w(g_i_)
          enddo
          t2 = 0.5d0 * (de2_11 * d4_v - de2_11)
C--------
          do g_i_ = 1, g_p_
            g_d1_w(g_i_) = (-beta3_11) * g_r3(g_i_)
          enddo
          d1_w = (-beta3_11) * (r3 - re3_11)
          d2_v = exp(d1_w)
          d2_p =  d2_v
          d4_v = (1.d0 + d2_v) * (1.d0 + d2_v)
          d1_p = 2.0d0 * (1.d0 + d2_v)
          d7_b = 0.5d0 * de3_11 * d1_p * d2_p
          do g_i_ = 1, g_p_
            g_t3(g_i_) = d7_b * g_d1_w(g_i_)
          enddo
          t3 = 0.5d0 * (de3_11 * d4_v - de3_11)
C--------
C
          do g_i_ = 1, g_p_
            g_t1(g_i_) = sp1_11 * g_t1(g_i_)
          enddo
          t1 = t1 * sp1_11
C--------
          do g_i_ = 1, g_p_
            g_t2(g_i_) = sp2_11 * g_t2(g_i_)
          enddo
          t2 = t2 * sp2_11
C--------
          do g_i_ = 1, g_p_
            g_t3(g_i_) = sp3_11 * g_t3(g_i_)
          enddo
          t3 = t3 * sp3_11
C--------
C
C  EXCHANGE AND COULOMB  (exch1, coul1, etc.)
          do g_i_ = 1, g_p_
            g_coul1(g_i_) = 0.5d0 * g_t1(g_i_) + 0.5d0 * g_s1(g_i_)
          enddo
          coul1 = 0.5d0 * (s1 + t1)
C--------
          do g_i_ = 1, g_p_
            g_coul2(g_i_) = 0.5d0 * g_t2(g_i_) + 0.5d0 * g_s2(g_i_)
          enddo
          coul2 = 0.5d0 * (s2 + t2)
C--------
          do g_i_ = 1, g_p_
            g_coul3(g_i_) = 0.5d0 * g_t3(g_i_) + 0.5d0 * g_s3(g_i_)
          enddo
          coul3 = 0.5d0 * (s3 + t3)
C--------
          do g_i_ = 1, g_p_
            g_exch1(g_i_) = (-0.5d0) * g_t1(g_i_) + 0.5d0 * g_s1(g_i_)
          enddo
          exch1 = 0.5d0 * (s1 - t1)
C--------
          do g_i_ = 1, g_p_
            g_exch2(g_i_) = (-0.5d0) * g_t2(g_i_) + 0.5d0 * g_s2(g_i_)
          enddo
          exch2 = 0.5d0 * (s2 - t2)
C--------
          do g_i_ = 1, g_p_
            g_exch3(g_i_) = (-0.5d0) * g_t3(g_i_) + 0.5d0 * g_s3(g_i_)
          enddo
          exch3 = 0.5d0 * (s3 - t3)
C--------
C
          d4_v = (exch1 - exch2) * (exch1 - exch2)
          d3_p = 2.0d0 * (exch1 - exch2)
          d7_v = (exch2 - exch3) * (exch2 - exch3)
          d2_p = 2.0d0 * (exch2 - exch3)
          d10_v = (exch3 - exch1) * (exch3 - exch1)
          d1_p = 2.0d0 * (exch3 - exch1)
          d5_b = d1_p + (-d2_p)
          d6_b = -d1_p + d3_p
          d10_b = d2_p + (-d3_p)
          do g_i_ = 1, g_p_
            g_w(g_i_) = d5_b * g_exch3(g_i_) + d10_b * g_exch2(g_i_) + d
     *6_b * g_exch1(g_i_)
          enddo
          w = d4_v + d7_v + d10_v
C--------
C
          do g_i_ = 1, g_p_
            g_d1_w(g_i_) = (-cc_11) * g_r3(g_i_) + (-cc_11) * g_r2(g_i_)
     * + (-cc_11) * g_r1(g_i_) + (-cb_11) * g_w(g_i_)
          enddo
          d1_w = (-cb_11) * w - cc_11 * (r1 + r2 + r3)
          d2_v = exp(d1_w)
          d1_p =  d2_v
          d3_b = ca_11 * d1_p
          do g_i_ = 1, g_p_
            g_c1(g_i_) = d3_b * g_d1_w(g_i_)
          enddo
          c1 = ca_11 * d2_v
C--------
C
          d2_b = 1.0d0 / 2.d0
          d5_b = d2_b * c1 + d2_b * c1
          do g_i_ = 1, g_p_
            g_d1_w(g_i_) = d5_b * g_c1(g_i_) + d2_b * g_w(g_i_)
          enddo
          d1_w = (w + c1 * c1) / 2.d0
          d2_v = sqrt(d1_w)

          if ( d1_w .gt. 0.0d0 ) then
             d1_p = 1.0d0 / (2.0d0 *  d2_v)
          else
             d1_p = 0.d0
          endif
          do g_i_ = 1, g_p_
            g_rootw(g_i_) = d1_p * g_d1_w(g_i_)
          enddo
          rootw = d2_v
C--------
C
          do g_i_ = 1, g_p_
            g_e(g_i_, i) = -g_rootw(g_i_) + g_coul3(g_i_) + g_coul2(g_i_
     *) + g_coul1(g_i_)
          enddo
          e(i) = coul1 + coul2 + coul3 - rootw + de2_11
C--------
          do g_i_ = 1, g_p_
            g_e(g_i_, i) = cevau * g_e(g_i_, i)
          enddo
          e(i) = e(i) * cevau
C--------
C
C  end the nt loop here:
        enddo
C
        return
      end
C
C**************************************************

C**************************************************

C  U22- upper diabat
C
      subroutine g_prepot22(g_p_)
C
        implicit double precision (a-h, o-z)
        dimension r(nt, 3), e(nt)
C
C     data is in bohr and eV units
C
        save
        integer g_pmax_
        parameter (g_pmax_ = 3)
        integer g_i_, g_p_, ldg_r, ldg_e
        double precision d6_b, d2_v, d7_b, d4_v, d2_p, d1_p, d1_w, g_r1(
     *g_pmax_), g_r(ldg_r, nt, 3), g_r2(g_pmax_)
        double precision g_r3(g_pmax_), g_d1_w(g_pmax_), g_s1(g_pmax_),
     *g_t2(g_pmax_), g_t3(g_pmax_), g_e(ldg_e, nt)
        integer g_ehfid
        save g_r1, g_r2, g_r3, g_d1_w, g_s1, g_t2, g_t3
        data de2_11 /4.3d0/, excite /0.36d0/, wgt /0.2d0/
        data de1_22 /3.9d0/, de2_22 /4.3d0/, de3_22 /0.4d0/
        data re1_22 /2.1d0/, re2_22 /2.1d0/, re3_22 /2.5d0/
        data beta1_22 /1.d0/, beta2_22 /1.0d0/, beta3_22 /1.5d0/
        data g_ehfid /0/
C
C
        if (g_p_ .gt. g_pmax_) then
          print *, 'Parameter g_p_ is greater than g_pmax_'
          stop
        endif
        cevau = 1.d0 / 27.2113961d0
C
        return
C
        entry g_pot22(g_p_, r, g_r, ldg_r, e, g_e, ldg_e, nt)
C
        do i = 1, nt
C
          do g_i_ = 1, g_p_
            g_r1(g_i_) = g_r(g_i_, i, 1)
          enddo
          r1 = r(i, 1)
C--------
          do g_i_ = 1, g_p_
            g_r2(g_i_) = g_r(g_i_, i, 2)
          enddo
          r2 = r(i, 2)
C--------
          do g_i_ = 1, g_p_
            g_r3(g_i_) = g_r(g_i_, i, 3)
          enddo
          r3 = r(i, 3)
C--------
C
          do g_i_ = 1, g_p_
            g_d1_w(g_i_) = (-beta1_22) * g_r1(g_i_)
          enddo
          d1_w = (-beta1_22) * (r1 - re1_22)
          d2_v = exp(d1_w)
          d2_p =  d2_v
          d4_v = (1.d0 - d2_v) * (1.d0 - d2_v)
          d1_p = 2.0d0 * (1.d0 - d2_v)
          d6_b = (-(de1_22 * d1_p)) * d2_p
          do g_i_ = 1, g_p_
            g_s1(g_i_) = d6_b * g_d1_w(g_i_)
          enddo
          s1 = de1_22 * d4_v - de1_22
C--------
          do g_i_ = 1, g_p_
            g_d1_w(g_i_) = (-beta2_22) * g_r2(g_i_)
          enddo
          d1_w = (-beta2_22) * (r2 - re2_22)
          d2_v = exp(d1_w)
          d2_p =  d2_v
          d4_v = (1.d0 + d2_v) * (1.d0 + d2_v)
          d1_p = 2.0d0 * (1.d0 + d2_v)
          d7_b = 0.5d0 * de2_22 * d1_p * d2_p
          do g_i_ = 1, g_p_
            g_t2(g_i_) = d7_b * g_d1_w(g_i_)
          enddo
          t2 = 0.5d0 * (de2_22 * d4_v - de2_22)
C--------
          do g_i_ = 1, g_p_
            g_d1_w(g_i_) = (-beta3_22) * g_r3(g_i_)
          enddo
          d1_w = (-beta3_22) * (r3 - re3_22)
          d2_v = exp(d1_w)
          d2_p =  d2_v
          d4_v = (1.d0 + d2_v) * (1.d0 + d2_v)
          d1_p = 2.0d0 * (1.d0 + d2_v)
          d7_b = 0.5d0 * de3_22 * d1_p * d2_p
          do g_i_ = 1, g_p_
            g_t3(g_i_) = d7_b * g_d1_w(g_i_)
          enddo
          t3 = 0.5d0 * (de3_22 * d4_v - de3_22)
C--------
C
          do g_i_ = 1, g_p_
            g_e(g_i_, i) = wgt * g_t3(g_i_) + wgt * g_t2(g_i_) + g_s1(g_
     *i_)
          enddo
          e(i) = s1 + (t2 + t3) * wgt + excite + de2_11
C--------
          do g_i_ = 1, g_p_
            g_e(g_i_, i) = cevau * g_e(g_i_, i)
          enddo
          e(i) = e(i) * cevau
C--------
C
C  end the nt loop here:
        enddo
C
        return
      end
C
C*****************************************

C  U12- diabatic coupling
C
      subroutine g_prepot12(g_p_)
C
        implicit double precision (a-h, o-z)
C      dimension r(nt,3),De(3,nt),e(nt)
        dimension r(nt, 3), e(nt)
C
        save
        integer g_pmax_
        parameter (g_pmax_ = 3)
        integer g_i_, g_p_, ldg_r, ldg_e
        double precision d1_p, d1_w, d4_b, d3_b, d2_v, g_r1(g_pmax_), g_
     *r(ldg_r, nt, 3), g_r2(g_pmax_), g_s1(g_pmax_), g_s2(g_pmax_)
        double precision g_d1_w(g_pmax_), g_tmp1(g_pmax_), g_tmp2(g_pmax
     *_), g_e(ldg_e, nt)
        integer g_ehfid
        save g_r1, g_r2, g_s1, g_s2, g_d1_w, g_tmp1, g_tmp2
        data a1 /0.8d0/, a2 /1.0d0/
        data re1_12 /2.2d0/, re2_12 /4.2d0/
        data thshift /45.d0/
        data scale /0.2d0/
        data g_ehfid /0/
C
C
        if (g_p_ .gt. g_pmax_) then
          print *, 'Parameter g_p_ is greater than g_pmax_'
          stop
        endif
        cevau = 1.d0 / 27.2113961d0
C
        return
C
        entry g_pot12(g_p_, r, g_r, ldg_r, e, g_e, ldg_e, nt)
C
        do i = 1, nt
C
          do g_i_ = 1, g_p_
            g_r1(g_i_) = g_r(g_i_, i, 1)
          enddo
          r1 = r(i, 1)
C--------
          do g_i_ = 1, g_p_
            g_r2(g_i_) = g_r(g_i_, i, 2)
          enddo
          r2 = r(i, 2)
C--------
          r3 = r(i, 3)
C
          thrad = dacos(-1.d0) / 180.d0 * thshift
C
          csth = dcos(thrad)
          sith = dsin(thrad)
C
          do g_i_ = 1, g_p_
            g_s1(g_i_) = sith * g_r2(g_i_) + csth * g_r1(g_i_)
          enddo
          s1 = (r1 - re1_12) * csth + (r2 - re2_12) * sith
C--------
          do g_i_ = 1, g_p_
            g_s2(g_i_) = (-csth) * g_r2(g_i_) + sith * g_r1(g_i_)
          enddo
          s2 = (r1 - re1_12) * sith - (r2 - re2_12) * csth
C--------
C
          d2_v = s1 * s1
          d1_p = 2.0d0 * s1
          d3_b = (-a1) * d1_p
          do g_i_ = 1, g_p_
            g_d1_w(g_i_) = d3_b * g_s1(g_i_)
          enddo
          d1_w = (-a1) * d2_v
          d2_v = exp(d1_w)
          d1_p =  d2_v
          do g_i_ = 1, g_p_
            g_tmp1(g_i_) = d1_p * g_d1_w(g_i_)
          enddo
          tmp1 = d2_v
C--------
          d2_v = s2 * s2
          d1_p = 2.0d0 * s2
          d3_b = (-a2) * d1_p
          do g_i_ = 1, g_p_
            g_d1_w(g_i_) = d3_b * g_s2(g_i_)
          enddo
          d1_w = (-a2) * d2_v
          d2_v = exp(d1_w)
          d1_p =  d2_v
          do g_i_ = 1, g_p_
            g_tmp2(g_i_) = d1_p * g_d1_w(g_i_)
          enddo
          tmp2 = d2_v
C--------
C
          d3_b = scale * tmp2
          d4_b = scale * tmp1
          do g_i_ = 1, g_p_
            g_e(g_i_, i) = d4_b * g_tmp2(g_i_) + d3_b * g_tmp1(g_i_)
          enddo
          e(i) = tmp1 * tmp2 * scale
C--------
C
          do g_i_ = 1, g_p_
            g_e(g_i_, i) = cevau * g_e(g_i_, i)
          enddo
          e(i) = e(i) * cevau
C--------
C
C   end the nt loop here:
        enddo
C
        return
      end
C*****************************************
C
C

max_stars_repo_name: kevinmooreiii/DiNT

max_issues_repo_name: kevinmooreiii/DiNT

max_forks_repo_name: kevinmooreiii/DiNT


text:
       subroutine renum
      common/gflmn/sz,ks,kp,met
       met=0
       return
       end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


text:
subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)

! Riemann solver for the elasticity equations in 2d, with varying
! material properties rho, lambda, and mu
!
! waves: 4
! equations: 5
! aux fields: 5

! Conserved quantities:
!       1 sigma_11
!       2 sigma_22
!       3 sigma_12
!       4 u
!       5 v

! Auxiliary variables:
!       1  density
!       2  lamda
!       3  mu
!       4  cp
!       5  cs

! Note that although there are 5 eigenvectors, one eigenvalue
! is always zero and so we only need to compute 4 waves.
!
! solve Riemann problems along one slice of data.
!
! On input, ql contains the state vector at the left edge of each cell
!           qr contains the state vector at the right edge of each cell
!
! Note that the i'th Riemann problem has left state qr(:,i-1)
!                                    and right state ql(:,i)
! From the basic clawpack routines, this routine is called with ql = qr
!
! This data is along a slice in the x-direction if ixy=1
!                            or the y-direction if ixy=2.
!
! Here it is assumed that auxl=auxr gives the cell values
! for this slice.
!
! On output, wave contains the waves,
!            s the speeds,
!            amdq the  left-going flux difference  A^- \Delta q
!            apdq the right-going flux difference  A^+ \Delta q
!
! Note that the waves are *not* in order of increasing lambda.
! Instead the 1- and 2-waves are the P-waves and the 3- and 4-waves
! are the S-waves.   (The 5th wave would have speed zero and is not computed.)

    implicit none

    integer, intent(in) :: ixy, maxm, meqn, mwaves, mbc, mx, maux
    double precision, intent(in) :: ql, qr, auxl, auxr
    double precision, intent(out) :: wave, s, amdq, apdq

    dimension wave( meqn, mwaves, 1-mbc:maxm+mbc)
    dimension    s(mwaves, 1-mbc:maxm+mbc)
    dimension   ql(meqn, 1-mbc:maxm+mbc)
    dimension   qr(meqn, 1-mbc:maxm+mbc)
    dimension apdq(meqn, 1-mbc:maxm+mbc)
    dimension amdq(meqn, 1-mbc:maxm+mbc)
    dimension auxl(maux, 1-mbc:maxm+mbc)
    dimension auxr(maux, 1-mbc:maxm+mbc)

    integer :: ksig11, ksig22, ku, kv, i, m
    double precision :: dsig11, dsig22, dsig12, du, dv
    double precision :: alamr, amur, bulkr, cpr, csr
    double precision :: alaml, amul, bulkl, cpl, csl
    double precision :: det, a1, a2, a3, a4

    ! set ku to point to  the component of the system that corresponds
    ! to velocity in the direction of this slice, kv to the orthogonal
    ! velocity.  Similarly ksig11 and ksig22 point to normal stresses.
    ! 3rd component is always shear stress sig12.

    if (ixy.eq.1) then
        ksig11 = 1
        ksig22 = 2
        ku = 4
        kv = 5
    else
        ksig11 = 2
        ksig22 = 1
        ku = 5
        kv = 4
    endif

    ! note that notation for u and v reflects assumption that the
    ! Riemann problems are in the x-direction with u in the normal
    ! direciton and v in the orthogonal direcion, but with the above
    ! definitions of ku and kv the routine also works with ixy=2

    ! split the jump in q at each interface into waves
    ! The jump is split into leftgoing waves traveling at speeds -cp, -cs
    ! relative to the material properties to the left of the interface,
    ! and rightgoing waves traveling at speeds +cp, +cs
    ! relative to the material properties to the right of the interface,

    do i = 2-mbc, mx+mbc
        dsig11 = ql(ksig11,i) - qr(ksig11,i-1)
        dsig22 = ql(ksig22,i) - qr(ksig22,i-1)
        dsig12 = ql(3,i) - qr(3,i-1)
        du = ql(ku,i) - qr(ku,i-1)
        dv = ql(kv,i) - qr(kv,i-1)

        ! material properties in cells i (on right) and i-1 (on left):

        alamr = auxl(2,i)
        amur = auxl(3,i)
        bulkr = alamr + 2.d0*amur
        cpr = auxl(4,i)
        csr = auxl(5,i)

        alaml = auxr(2,i-1)
        amul = auxr(3,i-1)
        bulkl = alaml + 2.d0*amul
        cpl = auxr(4,i-1)
        csl = auxr(5,i-1)

        ! P-wave strengths:
        det = bulkl*cpr + bulkr*cpl
        if (det.eq.0.d0) then
            write(6,*) 'det=0 in rpn2'
            stop
        endif
        a1 = (cpr*dsig11 + bulkr*du) / det
        a2 = (cpl*dsig11 - bulkl*du) / det

        ! S-wave strengths:
        det = amul*csr + amur*csl
        if (det.eq.0.d0) then
            ! no s-waves
            a3 = 0.d0
            a4 = 0.d0
        else
            a3 = (csr*dsig12 + amur*dv) / det
            a4 = (csl*dsig12 - amul*dv) / det
        endif

        ! 5th wave has velocity 0 so is not computed or propagated.

        ! Compute the waves.

        wave(ksig11,1,i) = a1 * bulkl
        wave(ksig22,1,i) = a1 * alaml
        wave(3,1,i)  = 0.d0
        wave(ku,1,i) = a1 * cpl
        wave(kv,1,i) = 0.d0
        s(1,i) = -cpl

        wave(ksig11,2,i) = a2 * bulkr
        wave(ksig22,2,i) = a2 * alamr
        wave(3,2,i)  = 0.d0
        wave(ku,2,i) = -a2 * cpr
        wave(kv,2,i) = 0.d0
        s(2,i) = cpr

        wave(ksig11,3,i) = 0.d0
        wave(ksig22,3,i) = 0.d0
        wave(3,3,i)  = a3*amul
        wave(ku,3,i) = 0.d0
        wave(kv,3,i) = a3*csl
        s(3,i) = -csl

        wave(ksig11,4,i) = 0.d0
        wave(ksig22,4,i) = 0.d0
        wave(3,4,i)  = a4*amur
        wave(ku,4,i) = 0.d0
        wave(kv,4,i) = -a4*csr
        s(4,i) = csr

        ! compute the leftgoing and rightgoing flux differences:
        ! Note s(i,1),s(i,3) < 0   and   s(i,2),s(i,4) > 0.
        do m=1,meqn
            amdq(m,i) = s(1,i)*wave(m,1,i) + s(3,i)*wave(m,3,i)
            apdq(m,i) = s(2,i)*wave(m,2,i) + s(4,i)*wave(m,4,i)
        enddo
    enddo

    return
end subroutine rpn2

max_stars_repo_name: ranocha/riemann

max_issues_repo_name: ranocha/riemann

max_forks_repo_name: ranocha/riemann


text:
c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine x_solve(rho_i, qs, square, u, rhs, nx, nxmax, ny, nz)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c
c     Performs line solves in X direction by first factoring
c     the block-tridiagonal matrix into an upper triangular matrix,
c     and then performing back substitution to solve for the unknow
c     vectors of each line.
c
c     Make sure we treat elements zero to cell_size in the direction
c     of the sweep.
c
c---------------------------------------------------------------------
#ifdef READEX_LEVEL_1
#include "readex.inc"
#endif
      include 'header.h'
      include 'work_lhs.h'

      integer nx, nxmax, ny, nz
      double precision rho_i (  0:nxmax-1,0:ny-1,0:nz-1),
     $                 qs    (  0:nxmax-1,0:ny-1,0:nz-1),
     $                 square(  0:nxmax-1,0:ny-1,0:nz-1),
     $                 u     (5,0:nxmax-1,0:ny-1,0:nz-1),
     $                 rhs   (5,0:nxmax-1,0:ny-1,0:nz-1)

      integer i,j,k,m,n,isize
#ifdef READEX_LEVEL_1

      READEX_REGION_DEFINE( x_slv )
#endif
c---------------------------------------------------------------------
c---------------------------------------------------------------------

      if (timeron) call timer_start(t_xsolve)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     This function computes the left hand side in the xi-direction
c---------------------------------------------------------------------
#ifdef READEX_LEVEL_1
      READEX_REGION_START( x_slv , "x_solve" , SCOREP_USER_REGION_TYPE_COMMON )
#endif

!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(n,m,i,j,k,isize)
!$OMP&  SHARED(dx5,dx4,dx3,dx2,dx1,tx2,tx1,dt,c1345,c3c4,con43,c1,c2,
!$OMP&         nx,ny,nz)
      isize = nx-1

c---------------------------------------------------------------------
c     determine a (labeled f) and n jacobians
c---------------------------------------------------------------------
!$OMP DO
      do k = 1, nz-2
         do j = 1, ny-2
            do i = 0, isize

               tmp1 = rho_i(i,j,k)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2
c---------------------------------------------------------------------
c
c---------------------------------------------------------------------
               fjac(1,1,i) = 0.d0
               fjac(1,2,i) = 1.d0
               fjac(1,3,i) = 0.d0
               fjac(1,4,i) = 0.d0
               fjac(1,5,i) = 0.d0

               fjac(2,1,i) = -(u(2,i,j,k) * tmp2 *
     >              u(2,i,j,k))
     >              + c2 * qs(i,j,k)
               fjac(2,2,i) = ( 2.d0 - c2 )
     >              * ( u(2,i,j,k) / u(1,i,j,k) )
               fjac(2,3,i) = - c2 * ( u(3,i,j,k) * tmp1 )
               fjac(2,4,i) = - c2 * ( u(4,i,j,k) * tmp1 )
               fjac(2,5,i) = c2

               fjac(3,1,i) = - ( u(2,i,j,k)*u(3,i,j,k) ) * tmp2
               fjac(3,2,i) = u(3,i,j,k) * tmp1
               fjac(3,3,i) = u(2,i,j,k) * tmp1
               fjac(3,4,i) = 0.d0
               fjac(3,5,i) = 0.d0

               fjac(4,1,i) = - ( u(2,i,j,k)*u(4,i,j,k) ) * tmp2
               fjac(4,2,i) = u(4,i,j,k) * tmp1
               fjac(4,3,i) = 0.d0
               fjac(4,4,i) = u(2,i,j,k) * tmp1
               fjac(4,5,i) = 0.d0

               fjac(5,1,i) = ( c2 * 2.0d0 * square(i,j,k)
     >              - c1 * u(5,i,j,k) )
     >              * ( u(2,i,j,k) * tmp2 )
               fjac(5,2,i) = c1 *  u(5,i,j,k) * tmp1
     >              - c2
     >              * ( u(2,i,j,k)*u(2,i,j,k) * tmp2
     >              + qs(i,j,k) )
               fjac(5,3,i) = - c2 * ( u(3,i,j,k)*u(2,i,j,k) )
     >              * tmp2
               fjac(5,4,i) = - c2 * ( u(4,i,j,k)*u(2,i,j,k) )
     >              * tmp2
               fjac(5,5,i) = c1 * ( u(2,i,j,k) * tmp1 )

               njac(1,1,i) = 0.d0
               njac(1,2,i) = 0.d0
               njac(1,3,i) = 0.d0
               njac(1,4,i) = 0.d0
               njac(1,5,i) = 0.d0

               njac(2,1,i) = - con43 * c3c4 * tmp2 * u(2,i,j,k)
               njac(2,2,i) =   con43 * c3c4 * tmp1
               njac(2,3,i) =   0.d0
               njac(2,4,i) =   0.d0
               njac(2,5,i) =   0.d0

               njac(3,1,i) = - c3c4 * tmp2 * u(3,i,j,k)
               njac(3,2,i) =   0.d0
               njac(3,3,i) =   c3c4 * tmp1
               njac(3,4,i) =   0.d0
               njac(3,5,i) =   0.d0

               njac(4,1,i) = - c3c4 * tmp2 * u(4,i,j,k)
               njac(4,2,i) =   0.d0
               njac(4,3,i) =   0.d0
               njac(4,4,i) =   c3c4 * tmp1
               njac(4,5,i) =   0.d0

               njac(5,1,i) = - ( con43 * c3c4
     >              - c1345 ) * tmp3 * (u(2,i,j,k)**2)
     >              - ( c3c4 - c1345 ) * tmp3 * (u(3,i,j,k)**2)
     >              - ( c3c4 - c1345 ) * tmp3 * (u(4,i,j,k)**2)
     >              - c1345 * tmp2 * u(5,i,j,k)

               njac(5,2,i) = ( con43 * c3c4
     >              - c1345 ) * tmp2 * u(2,i,j,k)
               njac(5,3,i) = ( c3c4 - c1345 ) * tmp2 * u(3,i,j,k)
               njac(5,4,i) = ( c3c4 - c1345 ) * tmp2 * u(4,i,j,k)
               njac(5,5,i) = ( c1345 ) * tmp1

            enddo
c---------------------------------------------------------------------
c     now jacobians set, so form left hand side in x direction
c---------------------------------------------------------------------
            call lhsinit(lhs, isize)
            do i = 1, isize-1

               tmp1 = dt * tx1
               tmp2 = dt * tx2

               lhs(1,1,aa,i) = - tmp2 * fjac(1,1,i-1)
     >              - tmp1 * njac(1,1,i-1)
     >              - tmp1 * dx1
               lhs(1,2,aa,i) = - tmp2 * fjac(1,2,i-1)
     >              - tmp1 * njac(1,2,i-1)
               lhs(1,3,aa,i) = - tmp2 * fjac(1,3,i-1)
     >              - tmp1 * njac(1,3,i-1)
               lhs(1,4,aa,i) = - tmp2 * fjac(1,4,i-1)
     >              - tmp1 * njac(1,4,i-1)
               lhs(1,5,aa,i) = - tmp2 * fjac(1,5,i-1)
     >              - tmp1 * njac(1,5,i-1)

               lhs(2,1,aa,i) = - tmp2 * fjac(2,1,i-1)
     >              - tmp1 * njac(2,1,i-1)
               lhs(2,2,aa,i) = - tmp2 * fjac(2,2,i-1)
     >              - tmp1 * njac(2,2,i-1)
     >              - tmp1 * dx2
               lhs(2,3,aa,i) = - tmp2 * fjac(2,3,i-1)
     >              - tmp1 * njac(2,3,i-1)
               lhs(2,4,aa,i) = - tmp2 * fjac(2,4,i-1)
     >              - tmp1 * njac(2,4,i-1)
               lhs(2,5,aa,i) = - tmp2 * fjac(2,5,i-1)
     >              - tmp1 * njac(2,5,i-1)

               lhs(3,1,aa,i) = - tmp2 * fjac(3,1,i-1)
     >              - tmp1 * njac(3,1,i-1)
               lhs(3,2,aa,i) = - tmp2 * fjac(3,2,i-1)
     >              - tmp1 * njac(3,2,i-1)
               lhs(3,3,aa,i) = - tmp2 * fjac(3,3,i-1)
     >              - tmp1 * njac(3,3,i-1)
     >              - tmp1 * dx3
               lhs(3,4,aa,i) = - tmp2 * fjac(3,4,i-1)
     >              - tmp1 * njac(3,4,i-1)
               lhs(3,5,aa,i) = - tmp2 * fjac(3,5,i-1)
     >              - tmp1 * njac(3,5,i-1)

               lhs(4,1,aa,i) = - tmp2 * fjac(4,1,i-1)
     >              - tmp1 * njac(4,1,i-1)
               lhs(4,2,aa,i) = - tmp2 * fjac(4,2,i-1)
     >              - tmp1 * njac(4,2,i-1)
               lhs(4,3,aa,i) = - tmp2 * fjac(4,3,i-1)
     >              - tmp1 * njac(4,3,i-1)
               lhs(4,4,aa,i) = - tmp2 * fjac(4,4,i-1)
     >              - tmp1 * njac(4,4,i-1)
     >              - tmp1 * dx4
               lhs(4,5,aa,i) = - tmp2 * fjac(4,5,i-1)
     >              - tmp1 * njac(4,5,i-1)

               lhs(5,1,aa,i) = - tmp2 * fjac(5,1,i-1)
     >              - tmp1 * njac(5,1,i-1)
               lhs(5,2,aa,i) = - tmp2 * fjac(5,2,i-1)
     >              - tmp1 * njac(5,2,i-1)
               lhs(5,3,aa,i) = - tmp2 * fjac(5,3,i-1)
     >              - tmp1 * njac(5,3,i-1)
               lhs(5,4,aa,i) = - tmp2 * fjac(5,4,i-1)
     >              - tmp1 * njac(5,4,i-1)
               lhs(5,5,aa,i) = - tmp2 * fjac(5,5,i-1)
     >              - tmp1 * njac(5,5,i-1)
     >              - tmp1 * dx5

               lhs(1,1,bb,i) = 1.d0
     >              + tmp1 * 2.d0 * njac(1,1,i)
     >              + tmp1 * 2.d0 * dx1
               lhs(1,2,bb,i) = tmp1 * 2.d0 * njac(1,2,i)
               lhs(1,3,bb,i) = tmp1 * 2.d0 * njac(1,3,i)
               lhs(1,4,bb,i) = tmp1 * 2.d0 * njac(1,4,i)
               lhs(1,5,bb,i) = tmp1 * 2.d0 * njac(1,5,i)

               lhs(2,1,bb,i) = tmp1 * 2.d0 * njac(2,1,i)
               lhs(2,2,bb,i) = 1.d0
     >              + tmp1 * 2.d0 * njac(2,2,i)
     >              + tmp1 * 2.d0 * dx2
               lhs(2,3,bb,i) = tmp1 * 2.d0 * njac(2,3,i)
               lhs(2,4,bb,i) = tmp1 * 2.d0 * njac(2,4,i)
               lhs(2,5,bb,i) = tmp1 * 2.d0 * njac(2,5,i)

               lhs(3,1,bb,i) = tmp1 * 2.d0 * njac(3,1,i)
               lhs(3,2,bb,i) = tmp1 * 2.d0 * njac(3,2,i)
               lhs(3,3,bb,i) = 1.d0
     >              + tmp1 * 2.d0 * njac(3,3,i)
     >              + tmp1 * 2.d0 * dx3
               lhs(3,4,bb,i) = tmp1 * 2.d0 * njac(3,4,i)
               lhs(3,5,bb,i) = tmp1 * 2.d0 * njac(3,5,i)

               lhs(4,1,bb,i) = tmp1 * 2.d0 * njac(4,1,i)
               lhs(4,2,bb,i) = tmp1 * 2.d0 * njac(4,2,i)
               lhs(4,3,bb,i) = tmp1 * 2.d0 * njac(4,3,i)
               lhs(4,4,bb,i) = 1.d0
     >              + tmp1 * 2.d0 * njac(4,4,i)
     >              + tmp1 * 2.d0 * dx4
               lhs(4,5,bb,i) = tmp1 * 2.d0 * njac(4,5,i)

               lhs(5,1,bb,i) = tmp1 * 2.d0 * njac(5,1,i)
               lhs(5,2,bb,i) = tmp1 * 2.d0 * njac(5,2,i)
               lhs(5,3,bb,i) = tmp1 * 2.d0 * njac(5,3,i)
               lhs(5,4,bb,i) = tmp1 * 2.d0 * njac(5,4,i)
               lhs(5,5,bb,i) = 1.d0
     >              + tmp1 * 2.d0 * njac(5,5,i)
     >              + tmp1 * 2.d0 * dx5

               lhs(1,1,cc,i) =  tmp2 * fjac(1,1,i+1)
     >              - tmp1 * njac(1,1,i+1)
     >              - tmp1 * dx1
               lhs(1,2,cc,i) =  tmp2 * fjac(1,2,i+1)
     >              - tmp1 * njac(1,2,i+1)
               lhs(1,3,cc,i) =  tmp2 * fjac(1,3,i+1)
     >              - tmp1 * njac(1,3,i+1)
               lhs(1,4,cc,i) =  tmp2 * fjac(1,4,i+1)
     >              - tmp1 * njac(1,4,i+1)
               lhs(1,5,cc,i) =  tmp2 * fjac(1,5,i+1)
     >              - tmp1 * njac(1,5,i+1)

               lhs(2,1,cc,i) =  tmp2 * fjac(2,1,i+1)
     >              - tmp1 * njac(2,1,i+1)
               lhs(2,2,cc,i) =  tmp2 * fjac(2,2,i+1)
     >              - tmp1 * njac(2,2,i+1)
     >              - tmp1 * dx2
               lhs(2,3,cc,i) =  tmp2 * fjac(2,3,i+1)
     >              - tmp1 * njac(2,3,i+1)
               lhs(2,4,cc,i) =  tmp2 * fjac(2,4,i+1)
     >              - tmp1 * njac(2,4,i+1)
               lhs(2,5,cc,i) =  tmp2 * fjac(2,5,i+1)
     >              - tmp1 * njac(2,5,i+1)

               lhs(3,1,cc,i) =  tmp2 * fjac(3,1,i+1)
     >              - tmp1 * njac(3,1,i+1)
               lhs(3,2,cc,i) =  tmp2 * fjac(3,2,i+1)
     >              - tmp1 * njac(3,2,i+1)
               lhs(3,3,cc,i) =  tmp2 * fjac(3,3,i+1)
     >              - tmp1 * njac(3,3,i+1)
     >              - tmp1 * dx3
               lhs(3,4,cc,i) =  tmp2 * fjac(3,4,i+1)
     >              - tmp1 * njac(3,4,i+1)
               lhs(3,5,cc,i) =  tmp2 * fjac(3,5,i+1)
     >              - tmp1 * njac(3,5,i+1)

               lhs(4,1,cc,i) =  tmp2 * fjac(4,1,i+1)
     >              - tmp1 * njac(4,1,i+1)
               lhs(4,2,cc,i) =  tmp2 * fjac(4,2,i+1)
     >              - tmp1 * njac(4,2,i+1)
               lhs(4,3,cc,i) =  tmp2 * fjac(4,3,i+1)
     >              - tmp1 * njac(4,3,i+1)
               lhs(4,4,cc,i) =  tmp2 * fjac(4,4,i+1)
     >              - tmp1 * njac(4,4,i+1)
     >              - tmp1 * dx4
               lhs(4,5,cc,i) =  tmp2 * fjac(4,5,i+1)
     >              - tmp1 * njac(4,5,i+1)

               lhs(5,1,cc,i) =  tmp2 * fjac(5,1,i+1)
     >              - tmp1 * njac(5,1,i+1)
               lhs(5,2,cc,i) =  tmp2 * fjac(5,2,i+1)
     >              - tmp1 * njac(5,2,i+1)
               lhs(5,3,cc,i) =  tmp2 * fjac(5,3,i+1)
     >              - tmp1 * njac(5,3,i+1)
               lhs(5,4,cc,i) =  tmp2 * fjac(5,4,i+1)
     >              - tmp1 * njac(5,4,i+1)
               lhs(5,5,cc,i) =  tmp2 * fjac(5,5,i+1)
     >              - tmp1 * njac(5,5,i+1)
     >              - tmp1 * dx5

            enddo

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     performs gaussian elimination on this cell.
c
c     assumes that unpacking routines for non-first cells
c     preload C' and rhs' from previous cell.
c
c     assumed send happens outside this routine, but that
c     c'(IMAX) and rhs'(IMAX) will be sent to next cell
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     outer most do loops - sweeping in i direction
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     multiply c(0,j,k) by b_inverse and copy back to c
c     multiply rhs(0) by b_inverse(0) and copy to rhs
c---------------------------------------------------------------------
            call binvcrhs( lhs(1,1,bb,0),
     >                        lhs(1,1,cc,0),
     >                        rhs(1,0,j,k) )

c---------------------------------------------------------------------
c     begin inner most do loop
c     do all the elements of the cell unless last
c---------------------------------------------------------------------
            do i=1,isize-1

c---------------------------------------------------------------------
c     rhs(i) = rhs(i) - A*rhs(i-1)
c---------------------------------------------------------------------
               call matvec_sub(lhs(1,1,aa,i),
     >                         rhs(1,i-1,j,k),rhs(1,i,j,k))

c---------------------------------------------------------------------
c     B(i) = B(i) - C(i-1)*A(i)
c---------------------------------------------------------------------
               call matmul_sub(lhs(1,1,aa,i),
     >                         lhs(1,1,cc,i-1),
     >                         lhs(1,1,bb,i))

c---------------------------------------------------------------------
c     multiply c(i,j,k) by b_inverse and copy back to c
c     multiply rhs(1,j,k) by b_inverse(1,j,k) and copy to rhs
c---------------------------------------------------------------------
               call binvcrhs( lhs(1,1,bb,i),
     >                        lhs(1,1,cc,i),
     >                        rhs(1,i,j,k) )

            enddo

c---------------------------------------------------------------------
c     rhs(isize) = rhs(isize) - A*rhs(isize-1)
c---------------------------------------------------------------------
            call matvec_sub(lhs(1,1,aa,isize),
     >                         rhs(1,isize-1,j,k),rhs(1,isize,j,k))

c---------------------------------------------------------------------
c     B(isize) = B(isize) - C(isize-1)*A(isize)
c---------------------------------------------------------------------
            call matmul_sub(lhs(1,1,aa,isize),
     >                         lhs(1,1,cc,isize-1),
     >                         lhs(1,1,bb,isize))

c---------------------------------------------------------------------
c     multiply rhs() by b_inverse() and copy to rhs
c---------------------------------------------------------------------
            call binvrhs( lhs(1,1,bb,isize),
     >                       rhs(1,isize,j,k) )

c---------------------------------------------------------------------
c     back solve: if last cell, then generate U(isize)=rhs(isize)
c     else assume U(isize) is loaded in un pack backsub_info
c     so just use it
c     after call u(istart) will be sent to next cell
c---------------------------------------------------------------------

            do i=isize-1,0,-1
               do m=1,BLOCK_SIZE
                  do n=1,BLOCK_SIZE
                     rhs(m,i,j,k) = rhs(m,i,j,k)
     >                    - lhs(m,n,cc,i)*rhs(n,i+1,j,k)
                  enddo
               enddo
            enddo

         enddo
      enddo
!$OMP END DO nowait
!$OMP END PARALLEL
#ifdef READEX_LEVEL_1
      READEX_REGION_STOP( x_slv )

#endif
      if (timeron) call timer_stop(t_xsolve)

      return
      end

max_stars_repo_name: readex-eu/readex-apps

max_issues_repo_name: readex-eu/readex-apps

max_forks_repo_name: readex-eu/readex-apps


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

!   REAL with complex arguments and with/without KIND argument

program accuracy
  implicit none
  complex(kind=8)   :: z8
  complex(kind=4)   :: z4
  real(kind=8)      :: result(8), expect(8)
  data expect/ 0.5555555555555556d0,0.5555555555555556d0,0.5555555555555556d0, &
	       0.5555555820465088d0,0.5555555820465088d0,0.5555555820465088d0, &
	       0.5555555820465088d0,0.5555555820465088d0/

  z8 = (0.5555555555555555555555555d0,0.55555555555555555555555d0)
  result(1) = real(z8)
  result(2) = dreal(z8)
  result(3) = real(z8,kind=8)
  result(4) = real(z8,kind=4)
  z4 = (0.5555555555555555555555555e0,0.55555555555555555555555e0)
  result(5) = real(z4)
  result(6) = dreal(z4)
  result(7) = real(z4,kind=8)
  result(8) = real(z4,kind=4)
!  print *,result
!  print *,expect
  call check(result,expect,16)
end program accuracy


!---------------------------------------------------!
! Copyright (c) 2017 Shunsuke A. Sato               !
! Released under the MIT license                    !
! https://opensource.org/licenses/mit-license.php   !
!---------------------------------------------------!
subroutine propagation
  use global_variables
  implicit none
  integer :: it

  it = 0
  call write_field(it)

  do it = 1,Nt

    call dt_evolve

    if(mod(it,100) == 0)call write_field(it)

  end do

end subroutine propagation