text:
!==========================================================================2
! MODULES FOR BIGSTICK
!===========================================================================
!
!  MODULES:
!       ------- BASICS AND SINGLE-PARTICLE SPACE ----

!    system_parameters: many-body space defined by Z,N, Jz, parity...
!    verbosity : flags for writing out information, file numbers, etc
!    bitstuff : # of bits per word assumed
!    sporbit:  single particle ORBITs & quantum numbers
!        key derived variable: orbqn
!    spstate:  single particle STATEs & quantum numbers
!        key derived variable: spsqn
!        [ NB: an orbit labeled by (nlj) while state by (nljm) ]
!      -------  HAIKUS ----------
!    W_info   : truncations on weighting W
!    spsections : partition of the s.p. space by W; foundation of templates
!    templates : partition of the many-body space by W; like configurations,
!                  only not as fine-grained
!    builder : derived types used in constructing templates
!    haiku_info: information on haikus, mostly limits on W
!    haikus  : half-slater determinants used to construct...full slater determinants
!    blocks : organization of haikus
!     ---------- BASIS -------------
!    sdlimits: q# limits on SDs, used in constructing basis
!    basis : basic information on basis, including dim and start arrays
!    sectors: organization of basis
!    ----------- HAIKU HOPS ---------------
!    bugs   : derived types for hops
!        --- CONSTRUCTING JUMPS --------
!    descendents:
!    geneologies:
!    chains
!       ----- JUMPS ------------
!    jumpdef
!    jumpNbody
!      ------- INTERACTION -----------------
!    pairdef
!    interaction
!    diagh
!      ------- APPLYING THE HAMILTONIAN
!    lanczos_info
!      ------- OUTPUT ---------------------
!===========================================================================
!
!------  Modules for MPI node information ----------------------------------

!-------------- OVERALL QUANTUM NUMBERS OF THE SYSTEMS----------------------
module system_parameters
  implicit none
  integer :: np(2) ! Number of particles
  integer :: jz
  integer :: iparity
  character(len=1) :: cparity
  integer :: Acore    ! a of the core

!...... INFORMATION FOR PARTICLE-HOLE CONJUGATION........
  logical :: phconj(2)
  integer :: npmax(2)
  integer :: npeff(2)  ! effective number of particles regarding particle-hole conjugation;
                       ! added in 7.5.6

end module system_parameters

!...........................................................................
module menu_choices
  implicit none
  character(len=3)  :: menu_char
  logical :: menu_dx_omathematica ! density in mathematica format
end module menu_choices

!------------------------------------------
! MODULE PROGRAM_INFO
!
! contains information about code
!

  module program_info
     implicit none

     character(6) :: version = '7.9.12'
     character(9) :: lastmodified = ' Mar 2021'

  end module program_info

!-------------------------------------------------------------------
!  MODULE VERBOSITY
!-------------------------------------------------------------------
!  Various flags for printing out information.
!  This allows one to easily print out or suppress information.
!
!  NOTE: verbose_X  means write out information to terminal
!        chatty_X  means write out information as being computed (rare)
!        print_X    means write to file
!        X_file     is unit number for file
!-----------------------------------------------------------------
module verbosity
  implicit none

  logical :: print4modelinfo       ! prints out information needed for modeling
  integer :: modelinfo_file    =  44
  logical :: verbose_winnow    = .false.
  logical :: print_groups      = .false.
  integer :: group_file        = 37
  logical :: print_hsps        = .false.
  integer :: hsps_file         = 37
  logical :: verbose_templates =  .false.
  logical :: verbose_haikus    = .false.
  logical :: print_haikus      = .false.
  integer :: haiku_file        =  37
  logical :: verbose_blocks    = .false.
  logical :: verbose_hops      = .false.
  logical :: print_hops        = .false.
  integer :: hop_file          =  38
  logical :: verbose_sectors   = .false.
  logical :: print_sectors     = .false.
  integer :: sector_file       = 37
  logical :: chatty_descent    = .false.
  logical :: chatty_links      = .false.

  logical :: print_genes       = .false.
  integer :: gene_file         = 38
  logical :: print_basis       = .false.
  integer :: basis_file        = 31

  integer :: xsd_file          = 32  ! added 4/2011 CWJ
  logical :: verbose_uncouple  = .false.

  logical :: verbose_jumps  = .false.

  logical :: print_matrix      = .false.
  integer :: matrix_file       = 55
  logical :: verbose_orthog    = .false.

  logical :: print_block_ops   = .false.
  logical :: print_node_load   = .false.

  logical :: print_sectorjumps = .false.    ! added 5/2010 CWJ
  integer :: sectorjump_file   = 39

  logical :: print_jumps = .false.   ! added 4/2011 CWJ
  integer :: jump_file = 46

  logical :: verbose_3body_readin = .false.
  logical :: verbose_3bodyjumps = .false.

  logical :: verbose_distrib = .false.  ! added 8/2011 CWJ
  logical :: verbose_bundles  = .true.  ! added 10/2011 CWJ

  logical :: verbose_fragments = .false.
  logical :: verbose_distro    =.false.

end module verbosity

!...........................................................................
module bitstuff
  implicit none
  integer max_bit_word         ! the max. # of usable bits in a
  parameter (max_bit_word=30)  ! given word to describe an SD
end module bitstuff

!---------------------------------------------------------------------------
!  declaring size of variables
!---------------------------------------------------------------------------
module precisions
  integer,parameter :: lanc_prec =    4
  integer,parameter :: basis_prec   = 8
  integer,parameter :: obs_prec     = 8
  integer,parameter :: egv_prec     = 8
  integer(4) :: MPI_lanc
end module precisions

!-----------------------------------------------------
!  INFORMATION ON TIME PER OPERATION
!  USEFUL FOR LOAD BALANCING
!  DEFAULTS CAN BE MODIFIED FOR A GIVEN MACHINE
!-------------------------------------------------------
module timeinfo

  implicit none
!.............. TIMING PER OPERATION (in ns)............
  integer :: timefile = 59

  logical :: modifytimewts = .true.   ! enables modification of timing
                                     ! by reading in timing from file "timebigstick.info"
                                     ! (if it exists; otherwise use defaults)
!   DEFAULT VALUES IN ABSENCE OF OTHER INFORMATION

  real(8),parameter :: tpoSPEdefault = 2.0
  real(8),parameter :: tpoPPdefault = 6.0
  real(8),parameter :: tpoNNdefault = 6.0
  real(8),parameter :: tpoPNdefault = 8.0
  real(8),parameter :: tpoPPPdefault = 6.0
  real(8),parameter :: tpoNNNdefault = 6.0
  real(8),parameter :: tpoPPNdefault = 8.0
  real(8),parameter :: tpoPNNdefault = 8.0

! added by shan for backward PN
  real(8),parameter :: tpoPNbdefault = 8.0
  real(8),parameter :: tpoPPbdefault = 8.0

! added by CWJ for backward PPN, PNN

real(8),parameter :: tpoPPNbdefault = 8.0
real(8),parameter :: tpoPNNbdefault = 8.0

!  INFORMATION FROM PAST, CURRENT RUNS
  real(8) :: tpoSPEprior, tpoSPE
  real(8) :: tpoPPprior, tpoPP
  real(8) :: tpoNNprior, tpoNN
  real(8) :: tpoPNprior, tpoPN
  real(8) :: tpoPPPprior, tpoPPP
  real(8) :: tpoNNNprior, tpoNNN
  real(8) :: tpoPPNprior, tpoPPN
  real(8) :: tpoPNNprior, tpoPNN
  real(8) :: tpoPPNbprior, tpoPPNb
  real(8) :: tpoPNNbprior, tpoPNNb

  real(8) :: tpoPPbprior, tpoPPb
  real(8) :: tpoPNbprior, tpoPNb

  logical :: priortimeinfo

end module timeinfo
!...........................................................................
!
!  GATHERS MISCELLANEOUS INFORMATION FOR END-OF-RUN REPORT
!
module reporter
   implicit none
   character(len=40) :: spfilename
   character(len=50) :: tbmefilename(50)    ! assumes no more than 50 input files
   real              :: tbmescale(50),spescale(50)      ! scaling of 2-body matrix elements,spes
   character(1)      :: tbmetype(50)        ! type of TBME
   integer           :: ntbmefiles          ! # of tbme files; needs to be initialized
   character(len=50) :: wfn_in_filename

end module reporter
!...........................................................................
module io
  use precisions
  implicit none

  character(55) :: outfile              ! filename for outputs
  logical       :: writeout
  integer       :: autoinputfile   = 8
  integer       :: resultfile      = 12
  integer       :: denresultfile   = 32
  integer       :: wfnfile         = 13
  integer       :: oldwfnfile      = 14
  integer       :: mdensfile       = 15 ! mathematica density file
  integer       :: logfile         = 36 ! human-readable log file
  integer       :: binlogfile      = 22 ! binary log file, used for speeding up creations
  integer       :: fragfile        = 46 ! file for fragment info
  integer       :: configfile      = 91  ! file for configuration occupation
  logical       :: auto_readin         ! flag to read in information from .wfn file
  logical       :: auto_input          ! flag to read from autoinput.bigstick file
  logical       :: ham_readin          ! flag to read in hamiltonian files
  logical       :: op_readin           ! flag to read in (nonscalar) one-body operator
  logical       :: strengthflag        ! flag to compute strength functions starting from a pivot
  logical       :: blockstrengthflag   ! flag to compute strength functions starting from a pivot block ADDED
  logical       :: strengthnorm        !  if TRUE then normalize output wfns for strength function ADDED 7.9.8
  logical      :: dotflag               ! if true, compute dot produce, else compute relative entropy adding 7.9.9

  logical       :: densityflag         ! to compute density matrices inline
  logical       :: trdensout           ! a flag to write to file suitable for Navratil's TRDENS code
  logical       :: greenflag           ! flag to apply resolve/Green's function to a pivot ADDED 7.8.8
  logical       :: restart_enabled = .true.  ! if set = .false., makes for easier restart
                                       ! else lanczos vectors automatically erased
  logical       :: write_wfn           ! flag to write wavefunctions to file (.wfn) (added by P.G.K)
  logical       :: get_JT               ! flag to enable calculation of J, T; added 7.6.8
  logical       :: baseASCII           ! write out .bas file in human-readable ASCII   added 7.9.1

  integer(kind=basis_prec) :: dimbasischeck

  ! KSM:  Added to support writing large files to scratch directories
  !       on HPC machines
  character (len=64) :: nersc_host
  character (len=1024) :: base_scratch_dir  ! like /tmp or $SCRATCH
  character (len=1024) :: scratch_dir       ! add on resultfile
  logical printouthamflag  ! flag for printing out nonzeroes of Hamiltonian
  logical printouttrans1flag  ! flag for printing out nonzeroes  of 1- body transition
  logical :: modeldensities  ! added 7.9.6, May2020

end module io

!--------------------------------------------------
!
! module for controlling hermiticity
!
!module hermit

!  logical hermitian

!end module hermit

!--------------------------------------------------------------------------
!  single-particle ORBIT information
!  key derived variable: orbqn
!---------------------------------------------------------------------------
module sporbit
   implicit none
   integer :: numorb(2)  ! # of s.p. orbits
   integer :: numorbmax	        ! max of numorb
   logical :: isoflag           ! in isospin formalism
   logical :: pnwtflag          ! flag to allow p,n to have different weights even if in isospin formalism
                                ! added in version 7.2.8 by CWJ
   logical :: spinless          ! for "spinless" fermions
   logical :: allsameparity     ! flag for all sp orbits same parity
   logical :: allsameW          ! flag for all sp orbits same W
!------------ CREATE A DEFINED TYPE-----------------------------------------
   type orb
      integer :: nr            ! radial quantum number
      integer :: j             ! 2 x j
      integer :: l             ! L
      integer :: par           ! parity = +/- 1
      integer :: w             ! excitation
   end type orb
   type (orb),allocatable :: orbqn(:,:) ! orbqn(species,iorb)

   character*200 :: sps_path    ! environment variable for searching for directory containing single particle orbits
   integer :: length_sps_path
   integer :: maxorblabel       ! added 7.4.1; used for error traps in reading in matrix element files

contains
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
!  function parmult
!
! "multiplies" parities where 1 = + and 2 = -
!
!  INPUT: par1,par2
!
      integer function parmult(par1,par2)

      implicit none
      integer par1,par2  ! input
!..........................................
      integer i1,i2,i3

      if(par1 < 1 .or. par2 > 2)then
        print*,' wrong parity ',i1,i2
        stop
      endif

      if(par2 < 1 .or. par2 > 2)then
        print*,' wrong parity (2)',i1,i2
        stop
      endif

      i1 = 3-2*par1
      i2 = 3-2*par2
      i3 = i1*i2
      parmult = (3-i3)/2

      return
      end function parmult
end module sporbit

!---------------------------------------------------------------------------
!  single-particle STATE information
!  key dervied variable: spsqn
!---------------------------------------------------------------------------
module spstate
  implicit none
  integer :: nsps(2) ! # of s.p. states
  integer :: nspsmax
!------------ CREATE A DEFINED TYPE-----------------------------------------
  type spst
     integer :: nr       ! radial quantum number
     integer :: j        ! 2 x j
     integer :: m        ! 2 x jz
     integer :: l        ! L
     integer :: w        ! excitation
     integer :: par      ! parity
     integer :: orb      ! orbit label
     integer :: group    ! which group -- labels by jz, w, par
     integer :: tr       ! time-reversed state   (added Aug 2010 by CWJ @ SDSU)
  end type spst
  type (spst),allocatable :: spsqn(:,:)  ! spsqn(ispecies,istate)

!--------- INFO ON HAIKU S.P. STATES:---------------------------------------
!  haikus are divided up into "left" built from s.p. states with jz < 0
!                         and "right" built from s.p. states with jz >= 0
!------------------ CONVENTION for hspsqn (and elsewhere)-------------------
!                   if ispecies < 0, then "left" haiku  jz < 0
!                               >0        "right" haiku  jz >= 0
!...........................................................................
  integer :: nhsps0(-2:2)
  integer :: nhspsmax
  type (spst),allocatable :: hspsqn(:,:) ! hspsqn(ispecies,istate)

!---------------------- GROUPS -- group together hsps by par, jz, par-------
  integer :: ngroups(-2:2)

  type grouplist
     integer, pointer :: jz(:)
     integer, pointer :: w(:)
     integer, pointer :: par(:)
     integer, pointer :: start(:)
     integer, pointer :: fin(:)
     logical, pointer :: active(:)
  end type grouplist

  type (grouplist) :: group(-2:2)
  logical :: groupwinnow = .true.  ! flag to winnow down groups
end module spstate

!---------------------------------------------------------------------------
!  Information on W (excitation weighting)
! (equivalent to Nhw in REDSTICK, t in ANTOINE)
!---------------------------------------------------------------------------
module W_info
  implicit none
  integer :: maxW(2)      ! max excitation allowed for p, n
  integer :: minW(2)      ! min excitation allowed for p, n
  integer :: maxWtot,minWtot     ! max,min total excitation
  integer :: wceiling(2)  ! ceiling to limit creating too many haikus
end module W_info

! NOTE: in 7.6.9 module spsections MOVED to file bsections.f90

!---------------------------------------------------------------------------
!  used to construct the basis
!---------------------------------------------------------------------------
!
! NOTE: in 7.6.9 module templates MOVED to file bbasis_template.f90
!
!...........................................................................
module builder
!
!  USEFUL (?) DERIVED TYPE  used in basislib1.f90 routines

  implicit none
  type alist
     integer,pointer :: loc(:)
  end type alist

  type blist
     type (alist), pointer :: list(:)
  end type blist

end module builder

!---------------------------------------------------------------------------
!  Information about haikus
!  a HAIKU is a slater determinant, constructed from s.p. states
!  "left haikus" built from s.p. states with jz < 0
!  "right haikus" built from s.p. states with Jz >= 0
!---------------------------------------------------------------------------
module haiku_info
  implicit none

                          ! nhsps is set in subroutine limit_hspspace in bspstates.f90
  integer             :: nhsps(-2:2) !these needs to be made consistent
                                                   ! and so does this:
  integer             :: Nword(2)           ! # of words needed to construct a haiku

  integer             :: minNh(-2:2) ! min number of particles in any haiku
  integer             :: maxNh(-2:2) ! max number of particles in any haiku
  integer,allocatable :: minWh(:,:)                ! min w found in haiku w/ n particles
  integer,allocatable :: maxWh(:,:)                ! maxWh(ispecies, nh); ispecies = +/- 1, 2
                                                   ! if ispecies -1,-2  = "left" haikus  jz < 0
                                                   !  if ispecies +1,+2 = "right" haikus jz >= 0
end module haiku_info

!  MODULE haikus moved to file bbasis_haikus.f90 (as of 7.6.9)

!  MODULE blocks moved to file bbasis_blocks.f90 (as of 7.6.9)

!  MODULE slaterlimits moved to file bbasis_limits.f90 in 7.6.9

!  MODULE basis moved to file bbasislib.f90 in 7.6.9

!---------------------------------------------------------------------------
!  sectors of the basis are divided by (first) jz, then parity, finally W
!  sectors are constructed from blocks of haikus
!---------------------------------------------------------------------------
module sectors

  implicit none
  integer     :: nsectors(2)                ! # of sections
  type sectorinfo
     integer(8)  :: xsdstart                          ! where, in the list of X-species sds, does this start
     integer(8)  :: xsdend
     integer(8)  :: nxsd                              ! # of SDs in this sector
	 integer(8)  :: ncxsd
     integer  :: jzX
     integer  :: parX
     integer  :: wX
     integer  :: ncsectors                         ! # of opposite-species sections can couple to
     integer, pointer :: csector(:)                ! list of opposite-species sections
     integer  :: nhblocks                          ! # of pairs of haiku blocks used to construct this section
     integer, pointer :: rhblock(:),lhblock(:)     ! list of right- and left-haiku blocks
     integer(8), pointer :: blockstart(:),blockend(:) ! where, for a given block, this starts and ends
     integer,pointer :: cut(:)                     ! ADDED 7.4.1 to be used to cut proton sectors for fragments
     integer         :: ncuts
     integer(kind=8) :: basisstart,basisend        ! ADDED 7.4.9: for protons, where basis starts and stops
     logical :: ibelong                            ! ADDED 7.5.9: denotes if a sector belongs on a given MPI process
                                                   ! especially when vectors are broken into fragments
  end type sectorinfo
  type mastersect
     type (sectorinfo), pointer :: sector(:)
  end type mastersect

  type (mastersect) :: xsd(2)

end module sectors

!.... MODULE bugs MOVED TO bhopslib.f90....(and renamed hoppy)

!.... MODULE nutple_info moved to bjumplib_ntuples.f90

!=============================================================
!  THE FOLLOWING MODULES ARE USED TO CONSTRUCT THE JUMPS
!==============================================================

!==============================================================

module geneologies
implicit none

type genlink
   type (genlink), pointer :: next
   integer, pointer :: gen(:)   ! stores list of block hops; if > 0 then right blockhop, if < 0 then left
   integer :: iblock,fblock
!-------------- NOTE: I MIGHT BE ABLE TO GET RID OF THESE THROUGH SOME REVISION
   integer :: irblock,ilblock   ! indices of initial, final blocks
   integer :: frblock,flblock   ! indices of initial, final blocks

end type genlink

type (genlink), target :: geneology

end module geneologies

!==============================================================
!
!  definition of derived types for sector jumps which overlay jumps
!
module jumpdef

implicit none

type conjump  ! information on conjugate jumps
   integer ncjmps
   integer, pointer :: cjump(:)
end type conjump

type jumper
   integer(8) :: njumps
   integer(8) :: nstart
!   integer(8) :: stride
   logical diag
!............. ADDED IN 7.5.3....  FOR INTRON DELETION SHIFTS.........
   logical :: containsexons     ! if "exons" (active jumps) for current MPI proc are here
   integer(8) :: exstart,exstop  ! stop/stop of exons
   integer(8) :: exshift        ! shift of exons for storage
!............. ADDED in 7.9.0.... FOR SPLITTING SECTOR JUMPS.........
!           these are filled out in bjumplib_master.f90
!           and then used to reconstruct the sector jumps
!
   integer(4) :: nsplit   ! number of splits
   integer(8),allocatable :: splitstart(:)   ! new starts for splits

end type jumper

type jumpsect
   integer :: nsectjumps
   integer, pointer :: isector(:)
   integer, pointer :: fsector(:)
   type (jumper), pointer :: sjmp(:)
   type (conjump), pointer :: csjmp(:)
end type jumpsect

end module jumpdef
!=====================================================

!
!  contains info on 1- and 2-body jumps
!
!  "JUMPS" are the fundamental constructs that define operations
!  a jump is always within a species and contains the following information:
!  -- index of initial and final (proton or neutron) Slater determinants
!  -- index/ices of operators that take us from initial to final Slater det
!     here either a 1-body (a^dagger a) or 2-body (a^dagger a^dagger a a) operator
!  -- phase information from anticommutation
!
! the master routines for creating jumps are found in bjumplib1.f90
!         master1bodyjumps  and master2bodyjumps
!  however jumps themselves are constructed (from "chains", an intermediate construct)
!  in subroutine weld in bjumplib3.f90
!

module jumpNbody
   use precisions
   use jumpdef
   implicit none

!---- MADE EXPLICIT IN 7.7.4 -- memory storage per jump
   integer(4) :: bytesper2Bjump = 25    ! 2 x 8 isd/fsd + 8 m.e. label + 1 phase
   integer(4) :: bytesper1Bjump = 21    ! 2 x 8 isd/fsd + 2 x 2 indexs + 1 phase

   logical :: MPIjumps  = .true.
   logical :: skipzeros =  .false.  !  flag to skip jumps with "known" zero matrix elements
                          ! this is straightforward for PP, NN jumps; not easy for others

!---------- SECTOR JUMPS define clumps of jumps between proton and neutron sectors
!           they are created in routines
!            set1bsectorjumps and set2bsectorjumps in bjumplib1.f90
!           and set3bsectorjumps in bjumplib5.f90
!
!           each sector jumps is associated with "conjugate" sector jumps
!           for the opposite species
!       "conjugate" sector jumps are computed in
!           masterfindconjumps1b and masterfindconjumps2b in bjumplib3.f90
!           and masterfindconjumps3b  in bjumplib5.f90
!

   type (jumpsect), target :: x1bjump(2)
   type (jumpsect),target :: x2bjump(2)
!-------------------
   ! one-body operators stay within species and include kinetic energy + external
   ! potential.  The external potential could be a sum over a core.
   ! We also have components of H_{xy} which have the same form.
   integer(8) :: totn1bjumps
   ! n1b_isd:  map neutron jump id to initial proton slater determinant
   ! n1b_fsd:  map neutron jump id to final proton slater determinant
   integer(kind=basis_prec),allocatable, target :: n1b_isd(:),n1b_fsd(:)
   integer(kind=basis_prec),allocatable, target :: n1b_isd0(:),n1b_fsd0(:)   ! used only when printing out jumps

   ! n1b_cop:  map neutron jump id to to creation op single particle label
   ! n1b_dop:  map neutron jump id to to destruction op single particle label
   ! n1b_phase maps proton jump id to phase to apply to matrix element
   integer(kind=2),allocatable, target :: n1b_cop(:),n1b_dop(:)   ! used for 1-body jumps
   integer(kind=1), allocatable, target :: n1b_phase(:)

   integer(8) :: totp1bjumps
   ! p1b_isd:  map proton jump id to initial proton slater determinant
   ! p1b_fsd:  map proton jump id to final proton slater determinant
   integer(kind=basis_prec),allocatable, target :: p1b_isd(:),p1b_fsd(:)
   integer(kind=basis_prec),allocatable, target :: p1b_isd0(:),p1b_fsd0(:)  ! used only when printing out jumps

   ! p1b_cop:  map proton jump id to to creation op single particle label
   ! p1b_dop:  map proton jump id to to destruction op single particle label
   !   they are used together    a^+_c a_d for hamiltonian entry
   ! p1b_phase maps proton jump id to phase to apply to matrix element
   integer(kind=2),allocatable, target :: p1b_cop(:),p1b_dop(:)   ! used for 1-body jumps
   integer(kind=1), allocatable, target :: p1b_phase(:)

   ! two body operators
   ! H_{xx} = (1/4) \sum {V_{ijkl}^{(xx)} {a^+_i(x) a^+_j(x) a_l(x) a_k(x)}   paper eqn 15
   ! must be diagonal in quantum numbers because species y is not involved.
   integer(8):: totn2bjumps,totp2bjumps
   integer(kind=basis_prec),allocatable, target :: n2b_isd(:),n2b_fsd(:),p2b_isd(:),p2b_fsd(:)
   integer(kind=basis_prec),allocatable, target :: n2b_isd0(:),n2b_fsd0(:),p2b_isd0(:),p2b_fsd0(:) ! used only when printing out jumps

   integer,allocatable, target :: n2b_op(:),p2b_op(:)
   real(4), allocatable, target :: n2b_me(:), p2b_me(:)
   integer(kind=1), allocatable, target :: n2b_phase(:),p2b_phase(:)

   logical :: onebodyonly    ! flag to signal only one-body; added 7.6.8
   logical :: alt_PNb_loop   ! flag to allow for alternate loop for PPb in mat-vec mult added 7.8.0

end module jumpNbody
!==============================================

!
!  INFORMATION ON COUPLED two-body pairs
!
!
      module pairdefcoupled

!
!  information on COUPLED two-body pairs
!   restricting on W (truncation #) for storage of coupled TBMEs
!

      type coupair_qn
        integer :: par
        integer :: indx
        integer :: ia,ib
      end type coupair_qn

      type coupairinfo
        type (coupair_qn), allocatable :: pairc(:)
        integer :: meref(2), mestart(2)  ! i parity
      end type coupairinfo

!--- ADDED 7.9.3--- ALLOWS FOR SORTING INTO BLOCKS WITH W----
      type coupairinfo_adv
        type (coupair_qn), allocatable :: pairc(:)
		integer :: nw,wstart,wstop
		integer,allocatable :: wlist(:)
        integer :: meref(3), mestart(3)  ! i parity
      end type coupairinfo_adv

      end module pairdefcoupled

!==============================================================
!
!  INFORMATION ON COUPLED MATRIX ELEMENTS (i.e., single -particle energies,
!   couple 2-body matrix elements, etc);
!
!  ENCODING for Hamiltonian matrix elements
!    Vxx(ab,cd):
!     a >= b, c >= d,  (ab) >= (cd)
!       pair1 = a*(a-1)/2 + b  -> XXcouplemap(pair1)
!      pair2 = c*(c-1)/2 + d   -> XXcouplemap(pair2)
!    pair1 >= pair2
!   iref = XXcouples(it)%meref(ipar)
!   istart = XXcouples(it)%mestart(ipar)
!   indx = (pair1-iref)*(pair1-1-iref)/2+pair2-iref+istart
!
!   Vpn(ab,cd)
!   a >= c; if a==c, b >= d
!      pair1 = numorb(2)*(a-1) + b  -> PNcouplemap(pair1)
!     pair2 = numorb(2)*(c-1) + d -> PNcouplemap(pair2)
!     iref = PNcouples%meref(ipar)
!     istart = PNcouples%mestart(ipar)
!    indx = (pair1-iref)*(pair1-1-iref)/2+pair2-iref+istart
!
!  ENCODING for two-body density matrix elements  ADDED 7.9.1
!    Vxx(ab,cd):
!     a >= b, c >= d
!       pair1 = a*(a-1)/2 + b  -> XXcouplemap(pair1)
!      pair2 = c*(c-1)/2 + d   -> XXcouplemap(pair2)
!    pair1 >= pair2
!   iref = XXcouples(it)%meref(ipar)
!   istart = XXcouples(it)%mestart(ipar)
!   indx = (pair1-iref)*numXXpairs(ipar) + pair2-iref+istart
!
!   Vpn(ab,cd)
!   a >= c; if a==c, b >= d
!      pair1 = numorb(2)*(a-1) + b  -> PNcouplemap(pair1)
!     pair2 = numorb(2)*(c-1) + d -> PNcouplemap(pair2)
!     iref = PNcouples%meref(ipar)
!     istart = PNcouples%mestart(ipar)
!    indx = (pair1-iref)*(pair1-1-iref)/2+pair2-iref+istart
!
!  REQUIRES MODULE pairdefcoupled

   module coupledmatrixelements
	   use pairdefcoupled
	   implicit none
!               restrict on total W

	    integer ncouplesXX(2)
	    type (coupairinfo) :: XXcouples(2)
	    integer, allocatable, target :: PPcouplemap(:), NNcouplemap(:)
	    integer ncouplesPN
	    type (coupairinfo) :: PNcouples
     integer, allocatable :: PNcouplemap(:)

     logical dens2bflag     ! flag for which encoding to use--Hamiltonian or 2-body density ADDED 7.9.1
	 logical diag_den2b     ! flag for writing out only 2-body densities for initial=final, only Jt=0  ADDED 7.9.5
!	 logical iso_den2b    ! flag to write out diagonal 2-body densities with good isospin (may leave out options)
!...........................................................................
!-----------   Logical variables to control how interaction matrix elements
!-----------   are put into pp, nn, and pn components when using isospin
!-----------   formalism
	     logical pp_int,nn_int,pn_int,iso_int,isov_int,isot_int, coul_int, mixed_T

!--------- INFORMATION ABOUT FOLDERS TO FIND MATRIX ELEMENTS

	      character*200 :: int_path
	      integer :: length_int_path

		  logical :: emptyHam     ! no interaction read
		  real  :: hw   ! for c.m.
		  real :: ephshift   ! overall shift in energy from particle-hole tranformation

		  character*55 intfilename
		  real,allocatable, target :: pspe(:),nspe(:)
		  real,allocatable, target :: psppot(:,:),nsppot(:,:)  ! added in 7.7.9; off-diagonal single-particle
		                                                       ! potentials
		  logical :: call_spe  ! flag whether or not to call SPE routines

		  integer :: nv2bmedim(0:2)  ! 0 = pn, 1 = pp, 2 = nn
		  integer :: norb_allow(2)  ! # of allowed orbits
		  integer :: nv2bme_tot

!----------- coupled  TBMES -----------------

		  type vjs
		     integer  :: Jmin,Jmax
		     real,pointer :: v(:)
		  end type vjs
		  type (vjs), allocatable,target :: ppme(:),pnme(:),nnme(:)
!----------- COUPLE two-body density matrices
          !
		  type vjs2b
		     integer  :: Jabmin,Jabmax,Jcdmin,Jcdmax,Jmin,Jmax
		     real,pointer :: v(:,:,:)
		  end type vjs2b
		  type (vjs2b), allocatable,target :: pp2bden(:),pn2bden(:),nn2bden(:)

   end module coupledmatrixelements

!==============================================================
!
!  INFORMATION ON UNCOUPLED two-body pairs
!
! REQUIRES MODULE pairdef
!
      module pairdef

      implicit none
!
!  information on UNCOUPLED two-body pairs | ia, ib : JM >
!
      type pair_qn
        integer :: M
        integer :: par
        integer :: W
        integer :: indx
        integer :: ia,ib  ! important: if ia, ib > 0, then m is >= 0
                          ! else if ia < 0, then ma < 0
        integer :: tr        ! index of time-reversed pair
        integer(1) :: trphase   ! phase from time-reversal
		integer :: wb   ! added 7.9.4 for ordering by W of destruction operator in a^+ b
      end type pair_qn

      type pairinfo
        type (pair_qn), pointer :: pair(:)
        integer(kind=8), pointer       :: meref(:,:), mestart(:,:)  ! indices are M, parity
		integer(kind=8),pointer :: block(:,:)  ! ADDED in 7.9.1, size of blocks
      end type pairinfo
      end module pairdef

!==============================================================
!  storage of uncoupled two-body matrix elements
!    key derived  variables :vtbme; hmatXX,hmatpp,pairXX,pairpn
module interaction

  use pairdef
  use pairdefcoupled
  use precisions
  implicit none

!--- # of uncoupled TBMEs ---------------------
!  integer(8) :: nmatpp,nmatnn  NOT USED ANY MORE
  integer(8) :: nmatpn
  integer(8) :: nmatXX(2)

!---------- uncoupled TBMEs ------------------

  real, allocatable, target :: hmatpp(:),hmatnn(:)
  real, allocatable, target :: obsmatpp(:),obsmatnn(:)
  real, allocatable,target :: hmatpn(:),obsmatpn(:)
  real, pointer      :: vmatpn(:)

!--------- for 2-body density matrices, need double precision because adding many small numbers
  real(kind=obs_prec), allocatable, target :: dmatpp(:),dmatnn(:),dmatpn(:)

!------------ coding of the pairs -------------
!  take states i,j,k,l  from list in hspsqn
!
!  We encode a_i a_j, or a^+_i a^+_j, as follows:
!
!  assume i >= j
!
!  then  (i,j) =>  i(i-1)/2 +j

  integer npairXX(2)
  type (pairinfo) :: XX2(2)
  integer :: npairpn
  type (pairinfo) :: PN2
  integer :: nrhophX(2)
  type (pairinfo) :: rhophX(2)  ! added 7.9.4  -- one-body particle-hole operator for protons, neutrons

  integer(kind=8), allocatable :: cpnpair(:,:),dpnpair(:,:)
  integer(kind=8), allocatable,target :: mappairPP(:),mappairNN(:),mappairpn(:)
  logical :: useTR = .false.  ! use time-reversal to restrict some matrix elements
			      ! that is, Vijkl = V-i -j -k -l
  integer(1), allocatable :: phasepnpair(:,:)

!C--------------- ENCODE particle-hole pairs

  type (pair_qn),allocatable,target :: php(:),phn(:)  !particle-hole

  integer nphp,nphn           ! # of particle-hole operators
  integer,allocatable,target :: mapphp(:),mapphn(:)

!------ SPECIALIZED FLAG-- FOR RANDOMIZED UNCOUPLED MATRIX ELEMENTS----
!  SET FALSE EXCEPT FOR SPECIAL CASES

  logical :: random_mscheme = .false.

end module interaction
!===========================================================================
!
!  information for deformed single-particle spaces; only used in special cases
!
module interaction_deformed
	implicit none
	logical :: deformed = .false.
	real, allocatable, target :: speXm(:,:)

end module interaction_deformed
!===========================================================================
! module diagh moved to bdiagh.f90

!===========================================================================
! Observables, densities, applied operators, etc.
!===========================================================================
module obs
  use precisions
  implicit none
  real(kind=obs_prec) :: xj2,xt2
  real(kind=4),allocatable :: energy(:),xjlist(:), xtlist(:)
  real(kind=obs_prec) :: xj2tmp,xt2tmp ! for global sum /allreduce/ (added by P.G.K.)
  logical :: twoobsflag

end module obs

!=======================================================================
!module densities moved to bdensities.f90
!
!===================================================================
!
! coupled matrix elements for transition operators   ADDED 7.5.2
! FOR application to a wavefunction
! (corresponding uncoupled matrix elements found in module onebodypot)

module opmatrixelements
	implicit none
	logical :: pnoperators    ! flag to signal if operators are in pn format or isospin format
	  real, allocatable :: op1bod(:,:)  ! reduced matrix elements for one-body operartor; to be deleted

	real,allocatable,target :: pop1bod(:,:),nop1bod(:,:)  ! singly-reduced matrix elements for proton, neutron 1-body operators
	integer :: jop                         ! J of transition operator
	integer :: top                         ! (optional) T of transition operator

    logical :: subtract_enabled=.false.       ! for MKGK application of r^2 Y00
    logical :: subtract

end module opmatrixelements

!===================================================================
!
! application of one-body (uncoupled) potentials, including self-consistent mean-field for starting
! the original coupled matrix elements found in module opmatrixelements
!
      module onebodypot

      implicit none

      real, pointer :: pPOT(:,:),nPOT(:,:)

      real, allocatable, target :: ppot_obs(:,:),npot_obs(:,:)
      real, allocatable, target :: ppot_h(:,:),npot_h(:,:)

      end module onebodypot

!===========================================================================
! LANCZOS_INFO
! NOTE some info found in module flagger in bmodule_flags.f90
!
!===========================================================================
module lanczos_info
  use precisions
  implicit none

  integer(4):: niter  ! # of iteractions
  integer nkeep  ! # iterations to keep
  integer :: actual_iterations   ! counts how many actual iterations have taken place
                                 ! added by CWJ 8/2010

  logical :: reorthog = .true.   ! logical flag for reorthogonalization
  logical :: storevectorsenabled = .false.  ! if false, store Lanczos vectors in RAM if possible
  logical :: writetodisk

  integer(4) :: lvec_file = 30
  integer :: coef_file = 21  ! file for alpha, beta coefficients
  real (kind = lanc_prec),allocatable :: lvec(:,:)
  character(2) :: lanczchar
  integer :: kpow  ! exponent for polynomial filtered Lanczos; otherwise == 1
!------------- CONVERGENCE CONTROL ------
  logical :: fixed_iter     ! whether or not to used fixed # of iterations
!--------------THICK RESTART ---------------------
  integer :: nthick_keep     ! # of vectors to keep when restarting
  logical :: thick_restart
  integer :: max_num_thick_restart  ! max # of times to thick-restart
  real(8) :: Etarget        ! target energy for finding excited states
  logical :: targetX
!----------- CONSTRUCTING INITIAL PIVOT --------------  ADDED in 7.6.3 ----
  logical :: initializing_enabled=.false.
  integer :: initial_maxiter     ! max # of initial iterations
  logical :: initializing        ! flag to signal initializing, so don't do everything
  logical :: applypntrace         ! apply  average pn diagonal
  logical :: applyXXonly         ! only apply PP,NN, not PN
  logical :: diagonalsectorsonly ! only apply to "diagonal" sector jumps
  integer :: nsave               ! # of vectors to keep for starting, typically 1 but could be more

!---------- MISC (MOSTLY OBSOLETE)
  logical :: coef_only    ! used to get just coefficients
  integer :: time_loop_max = 6400  ! # of loops for timing

end module lanczos_info

!===========================================================================
! TIMING
!===========================================================================
module timing
  implicit none

  real(8) :: startall, endall
  real(8) :: startbasis,endbasis
  real(8) :: startham,endham
  real(8) :: startlanczos,endlanczos
  real(8) :: startobs,endobs,startwriteobs,endwriteobs
  real(8) :: startdens,enddens
  real(8) :: starts1b,ends1b
  real(8) :: startp2b,endp2b
  real(8) :: startn2b,endn2b

  real(8) :: time1body
  real(8) :: time2body
  real(8) :: timereorthog

! Added by PGK ( extra timing ).............................................
  real(8) :: time_applobs,time_writeobs
  real(8) :: time_reduce
  real(8) :: time_gather
  real(8) :: time_distr
  real(8) :: time_writeeigvec
  real(8) :: time_read_orthog
  real(8) :: time_pn
  real(8) :: time_pp
  real(8) :: time_nn
  real(8) :: time_ppb,time_pnb  ! 'backwards'
  real(8) :: time_pro
  real(8) :: time_dot
  real(8) :: time_swp
  real(8) :: time_spe
  real(8) :: time_rre     ! reduce (reorthog. only)
  real(8) :: time_ror     ! read in reorthog.
  real(8) :: time_sort1b     ! read in reorthog.
  real(8) :: time_sort2b     ! read in reorthog.
  real(8) :: time_hops    ! time to generate hops
  real(8) :: time_ppp, time_nnn, time_pnn, time_ppn
  real(8) :: time_ppnf,time_ppnb, time_pnnf,time_pnnb  ! added  7.9.2
  real(8) :: time_jmpcnt, time_jmpcnt_last   ! for counting jumps
  real(8) :: time_meset,  time_meset_last    ! set up for matrix elements
  real(8) :: time_munc,  time_munc_last    ! uncouple matrix elements
  real(8) :: time_desc,  time_desc_last     ! compute descendents setting up for jumps
  real(8) :: time_intron,  time_intron_last    ! map and delete introns
  real(8) :: time_pivot, time_pivot_last       ! time to set up pivot
  real(8) :: time_p1b,time_p1b_last, time_n1b,time_n1b_last  ! time for one-body densities

  real(8) :: time_hmult
  real(8) :: time_eig   ! time to get eigenvalues
  real(kind=8) :: time_eigvec
  real(kind=8) :: time_trdens

! Separate "timelast" variables are necessary to prevent "time-mixing"......
  real(8) :: timelast1    ! 1-body
  real(8) :: timelast2    ! 2-body
  real(8) :: timelast3    ! reorthog
  real(8) :: timelast4    ! applobs
  real(8) :: timelast5    ! reduce (applyh only)
  real(8) :: timelast6    ! distribution
  real(8) :: timelast7    ! writeeigvec
  real(8) :: timelast8    ! read_orthog
  real(8) :: timelast9    ! gather
  real(8) :: time_pn_last
  real(8) :: time_pp_last ! 2-body (pp)
  real(8) :: time_nn_last ! 2-body (nn)
  real(8) :: time_ppb_last  ! backwards pp
  real(8) :: time_pnb_last ! backwards pn
  real(8) :: time_ppp_last, time_ppn_last, time_pnn_last, time_nnn_last
  real(8) :: time_ppnf_last, time_ppnb_last,time_pnnf_last,time_pnnb_last ! added 7.9.2
  real(8) :: time_hmult_last
  real(8) :: time_pro_last! project Lanczos vector
  real(8) :: time_dot_last! dot product
  real(8) :: time_swp_last! "swap" v <---> w
  real(8) :: time_spe_last! find eigen-values
  real(8) :: time_rre_last! reduce (reorthog. only)
  real(8) :: time_ror_last! read in reorthog.
  real(8) :: time_eigvec_last   ! build eigenvectors
  real(8) :: time_trd_last   ! build trdens eigenvectors
  real(8) :: time_eig_last

! Temporary timing .........................................................
  real(8) :: t1start,t1end
  real(8) :: t8start,t8end
  real(8) :: t2,t2last
  real(8) :: t3,t3last
  real(8) :: t4,t4last
  real(8) :: t5,t5last
  real(8) :: t6,t6last
  real(8) :: t7,t7last

  real(8) :: jumptime,jumptimestart   ! used to check on timing across MPI processes
  real(8) ::time_hmult_wait,time_wait_last    ! USED ONLY FOR DEBUGGING, TIME WAITING TO REDUCE
!...........................................................................
end module timing

!===========================================================================
!
! for computing trace
!
module tracy
  implicit none
  integer(8) :: noperation, ioperation

  integer, allocatable :: indx1(:), indx2(:)
  real, allocatable :: me(:)

  type mex
     integer(kind=8) :: ncol
     real, allocatable :: v(:)
     integer(kind=8), allocatable :: col(:)
  end type mex

  type (mex), allocatable :: metrace(:)

  logical :: tracecount

!.... ADDED IN 7.6.3: COMPUTE THE PN TRACE SECTOR BY SECTOR
!           IN 7.6.6: TRACE (actually: average) OF ALL SECTOR BY SECTOR
   real, allocatable :: sectortrace(:)
!....... ADDED IN 7.9.0.... so one only has to compute the diagonal...
   logical :: getvariance = .false.
   real(kind=8) :: tr1,tr2

end module tracy

!===================================================================

max_stars_repo_name: noctildon/pyBigstick

max_issues_repo_name: noctildon/pyBigstick

max_forks_repo_name: noctildon/pyBigstick


text:
SUBROUTINE contac(task,iwrit,dtcal,ttime,velnp,maxve)

  !     main contac routine

  USE ctrl_db, ONLY: mconts, nconp, numct, ctype, cactive, top, bottom
  USE outp_db, ONLY: nreqc,nprqc
  IMPLICIT NONE

  !        Dummy arguments
  CHARACTER(len=*),INTENT(IN) :: task  !task to perform
  INTEGER (kind=4),INTENT(IN) :: iwrit
  INTEGER (kind=4),INTENT(IN), OPTIONAL :: maxve
  REAL (kind=8),INTENT(IN), OPTIONAL :: dtcal,ttime,velnp(:,:)
  !        Local variables
  INTEGER (kind=4) :: i

  INTERFACE
    INCLUDE 'conta1.h'
    INCLUDE 'conta2.h'
    INCLUDE 'conta3.h'
    INCLUDE 'dbea3d.h'
  END INTERFACE

  IF( TRIM(task) == 'DUMPIN' .OR.  &
      TRIM(task) == 'RESTAR' .OR.  &
      TRIM(task) == 'UPDLON')THEN

    IF(cactive(1)) CALL conta1(task,dtcal,ttime,iwrit,velnp,maxve)
    IF(cactive(2)) CALL conta2(task,dtcal,ttime,iwrit,velnp,maxve)
    IF(cactive(3)) CALL conta3(task,dtcal,ttime,iwrit,velnp,maxve)
    IF(cactive(4)) CALL dbea3d(task,ttime,iwrit)

    IF(TRIM(task) == 'UPDLON')THEN
      !  check
      numct = 0                 !initializes
      DO i=1,mconts                 !for each contact algorithm
        cactive(i) = nconp(i) > 0
        IF( cactive(i) ) THEN       !if there are contact pairs using this alg.
          numct = numct + 1         !increase the number of cont algors used
          ctype (numct) = i         !keep for this set the element type
        END IF
      END DO
      IF( numct == 0 ) THEN
        nreqc = 0
        IF( ASSOCIATED(nprqc) )DEALLOCATE(nprqc)
        bottom = .FALSE. ; top = .FALSE.
      END IF
    END IF

  ELSE

    DO i = 1, numct

      SELECT CASE (ctype(i))

      CASE (1)
        CALL conta1(task,dtcal,ttime,iwrit,velnp,maxve)
      CASE (2)
        CALL conta2(task,dtcal,ttime,iwrit,velnp,maxve)
      CASE (3)
        CALL conta3(task,dtcal,ttime,iwrit,velnp,maxve)
      CASE (4)
        CALL dbea3d(task,ttime,iwrit)

      END SELECT
    END DO
  END IF
  RETURN
END SUBROUTINE contac

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


text:
program preprocessor11
! #if with >, <, ==
implicit none

#define X 1
#define Y 5
#define Z 1+Y*3

print *, Z

#if X == 1
print *, "X is 1"
#else
print *, "X is not 1"
#endif

#if X == 2
print *, "X is 2"
#else
print *, "X is not 2"
#endif

#if X != 1
print *, "X is not 1"
#else
print *, "Not X is not 1"
#endif

#if X == 1 && Y == 5
print *, "X is 1 and Y is 5"
#else
print *, "Not (X is 1 and Y is 5)"
#endif

#if X == 2 && Y == 5
print *, "X is 2 and Y is 5"
#else
print *, "Not (X is 2 and Y is 5)"
#endif

#if X < 3 && Y > 3
print *, "X<3 and Y>3"
#else
print *, "Not (X<3 and Y>3)"
#endif

#if X > 3 && Y < 3
print *, "X>3 and Y<3"
#else
print *, "Not (X>3 and Y<3)"
#endif

#if X == 3-2
print *, "X is 3-2"
#else
print *, "X is not 3-2"
#endif

#if X == 3-1
print *, "X is 3-1"
#else
print *, "X is not 3-1"
#endif

#if Y == (6*5-5)/5
print *, "Y is (6*5-5)/5"
#else
print *, "Y is not (6*5-5)/5"
#endif

#if Y == (6*5-6)/5
print *, "Y is (6*5-6)/5"
#else
print *, "Y is not (6*5-6)/5"
#endif

#if Y == (-5+6*5)/5
print *, "Y is (-5+6*5)/5"
#else
print *, "Y is not (-5+6*5)/5"
#endif

#if Y == (-6+6*5)/5
print *, "Y is (-6+6*5)/5"
#else
print *, "Y is not (-6+6*5)/5"
#endif

#if Y == (+5+6*5)/7
print *, "Y is (+5+6*5)/7"
#else
print *, "Y is not (+5+6*5)/7"
#endif

#if Y == (+4+6*5)/7
print *, "Y is (+4+6*5)/7"
#else
print *, "Y is not (+4+6*5)/7"
#endif

#if Y == (+5+6*Y)/7
print *, "Y is (+5+6*Y)/7"
#else
print *, "Y is not (+5+6*Y)/7"
#endif

#if Z == 16
print *, "Z is 16"
#else
print *, "Z is not 16"
#endif

#if defined(Y) && (3 <= Y)
print *, "Y defined and 1 <= Y"
#else
print *, "Y not defined or not 1 <= Y"
#endif

#if defined(A) && (3 <= A)
print *, "A defined and 1 <= A"
#else
print *, "A not defined or not 1 <= A"
#endif

end program

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


text:
FUNCTION FUNC3(FUNC1, FUNC2, x, y)
  REAL, EXTERNAL :: FUNC1, FUNC2
  REAL :: FUNC3
  REAL :: x, y

  FUNC3 = FUNC1(x) * FUNC2(y)
END FUNCTION FUNC3

ext: f

max_stars_repo_name: mullikine/RosettaCodeData

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
          program demo_sinh
          use, intrinsic :: iso_fortran_env, only : real_kinds, real32, real64, real128
          implicit none
            real(kind=real64) :: x = - 1.0_real64
            x = sinh(x)
          end program demo_sinh

max_stars_repo_name: urbanjost/fortran-intrinsic-manpages

max_issues_repo_name: urbanjost/fortran-intrinsic-manpages

max_forks_repo_name: urbanjost/fortran-intrinsic-manpages


text:
      SUBROUTINE SW_NEXT
     &          (pflag,nsgmnt0,nsgmnt)

c     Locates the next segment which matches the one just done

c pflag          profile index flag
c nsgmnt0        index of starting segment
c nsgmnt         index of current segment

      parameter     (mxsgmnt=201)

      common/mf_sw_1/dst(mxsgmnt),xla(mxsgmnt),xlo(mxsgmnt),
     &               azm(mxsgmnt),xdp(mxsgmnt),fld(mxsgmnt),
     &               sgm(mxsgmnt),eps(mxsgmnt),ncd(mxsgmnt),
     &               bta(mxsgmnt),hpr(mxsgmnt),npr(mxsgmnt),
     &               num(mxsgmnt),nrsgmnt

      logical        end loop
      integer        pflag

      end loop=.false.
      do while (.not.end loop)
         nsgmnt=nsgmnt+1
         if (nsgmnt .gt. nrsgmnt) then

            end loop=.true.
         else

            if (sgm(nsgmnt) .eq. sgm(nsgmnt0) .and.
     &          num(nsgmnt) .eq. 0) then

               if (pflag .gt. 1) then

                  if (pflag .lt. 5) then

                     if (hpr(nsgmnt) .eq. hpr(nsgmnt0)) end loop=.true.
                  else

                     if (npr(nsgmnt) .eq. npr(nsgmnt0)) end loop=.true.
                  end if
               else

                  end loop=.true.
               end if
            end if
         end if
      end do

      RETURN
      END      ! SW_NEXT

ext: for

max_stars_repo_name: spinorkit/LWPC

max_issues_repo_name: spinorkit/LWPC

max_forks_repo_name: spinorkit/LWPC


! module DoLoop  ! It is recommended to always place functions and subroutines within modules.
!     implicit none
!     private

!     public pi, radius, height, volume

! contains
!     subroutine CalculateVolume(radius, height)
!         implicit none

!         real(dp),intent(in):: pi, radius, height, volume

!         volume = pi * radius**2 * height
!     end subroutine CalculateVolume

! end module DoLoop

! !First, every module starts with comments documenting the purpose and content of the module.
! !https://fortran-lang.org/learn/best_practices/modules_programs

! !Second, Similarly, every procedure starts with a comment briefly describing its purpose and the intent of the dummy arguments.

! !3rd,imports (use) and exports (public) are explicitly given.


text:
I created this account to let everyone know how unsatisfied I am with Clarity Massage on G Street in Davis.  Below is the posting I wrote for the Clarity Massage page:

20080610 15:49:20   I had a terrible experience at Clarity Massage. I have had many massages in my life with different masseuses and I have never been so dissatisfied as I was this time. I felt extremely uncomfortable and entirely exposed the entire time I was there. I asked Stuart to stop early and I left the office as quickly as possible. His inappropriate caresses and toocloseforcomfort rubbing left me feeling dirty and takenadvantageof. I was never adequately covered and the placement of my body as well as his was unprofessional and sickening. I DO NOT recommend Clarity Massage to anyone, especially if you are unfamiliar with massages as you may not know what is appropriate and what is not while you are there. The massage I received was entirely inappropriate and I do not want anyone else to risk having the same experience I had. I have not tried any other massage services in Davis yet and I have never posted on Davis Wiki, but my experience at Clarity Massage was horrific enough for me to create a Davis Wiki account to let everyone know about my experience there. Please be careful when choosing a masseuse and always speak up if you feel uncomfortable. —UnsatisfiedWithClarityMassage

It took me a few weeks to get over my experience and post something on DavisWiki (which is where I learned of Clarity Massage in the first place).  There is still only one other comment posted on his page, and it was the comment that made me think I would have a great experience at Clarity Massage.  Because I was so very misinformed, I felt I needed to post something as well to let others know about my experience there.  If there is anyone else who has had a similar experience at Clarity Massage, I hope you were able to leave the office with some dignity.  If not, Im so sorry for what you went through, and I hope you will have the courage to let someone know and post something here about your experience so others dont have to go through what we did.

20080610 23:14:56 nbsp If you were touched inappropriately, the best thing to do would be to report this to the police. Its very hard to give much credence to an anonymous poster making unverifiable claims who might have a grudge against the person in question. The police are better equipped to deal with this kind of thing, anyway. Just a thought.  Users/WilliamLewis

20080626 16:12:04 nbsp I was considering sending my girlfriend to Clarity Massage, but your comment has me spooked. Users/KellyM

20080723 09:57:14 nbsp KellyM, it helps so much to know that my experience may have helped you and your friend make the most informed decision about where to get a massage.  Other action has been taken against Stuart Steinman, and my comments on this page were meant only to alert other potential customers of the experience I had so they could better make their own decision about going to Clarity Massage.  To WilliamLewis, your response to my comment being untruthful and unverifiable because I prefer to remain anonymous was hurtful and I can only hope that someone else in a similar situation would never attempt to ask you for help about it.  You dont know what other action I have already taken and I prefer to stay anonymous out of precaution for my own safety and to avoid persecution from users like you. Users/UnsatisfiedWithClarityMassage
 I never said you were untruthful, I only said the manner in which youre doing this is not the way to go. Also, where am I persecuting you? Users/WilliamLewis

20080723 10:38:03 nbsp I removed your comment from the Clarity Massage page because it alleged grossly illegal activity in a way that is essentially unverifiable. Such a thing would be ok if you stood behind your real name. It would also be ok if you were a long time editor of this site with a history of good edits. However, you came in here, created an anonymous account for a single purpose, and then trashed someone in a way that could have serious legal consequences to someone. That looks rather shady. Weve had a major problem with anonymous users coming on here and bashing various businesses. Weve often been able to find out that the anonymous users were associated with competing businesses. If you stood behind your real name, you would be providing a great measure of assurance that your claims are honest ones. For if they were not, you would be risking a libel lawsuit. Users/WilliamLewis

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C**************************************************************************
C      FFTPACK: Available at http://www.scd.ucar.edu/softlib/mathlib.html
C
C   ***NOTE (C.T.) *** this is a subset of FFTPACK, which only includes
C            the routines for the complex FFT (forward & inverse).
C
C Modified: November 1999 by Arjan van Dijk to include IMPLICIT NONE and
C           to convert all routines to DOUBLE precision.
C**************************************************************************

C      * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C      *                                          *
C      *                        FFTPACK             *
C      *                                          *
C      *     A package of Fortran subprograms for calculating          *
C      *     fast Fourier transforms for both complex and real         *
C      *      periodic sequences and certain other symmetric           *
C      *             sequences that are listed below            *
C      *               (Version 4.1 November 1988)                *
C      *                          by                    *
C      *                     Paul Swarztrauber                   *
C      *                          of                    *
C      *         The National Center for Atmospheric Research         *
C      *               Boulder, Colorado  (80307)  U.S.A.           *
C      *                  which is sponsored by                  *
C      *              the National Science Foundation          *
C      *                                          *
C      * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C  Any source code available through this distribution interface at NCAR
C  is free of charge, but there is no guarantee.
C
C FFTPACK breaks the FORTRAN 77 ANSI Standard
C by passing REAL arrays to subroutines and using the arrays within
C the subroutines as DOUBLE PRECISION or other types.  This infraction
C may cause data alignment problems when the source code is compiled
C and loaded in an executable.

C****************************************************************************
C     SUBROUTINE CFFTI(N,WSAVE)
C
C     SUBROUTINE CFFTI INITIALIZES THE ARRAY WSAVE WHICH IS USED IN
C     BOTH CFFTF AND CFFTB. THE PRIME FACTORIZATION OF N TOGETHER WITH
C     A TABULATION OF THE TRIGONOMETRIC FUNCTIONS ARE COMPUTED AND
C     STORED IN WSAVE.
C
C     INPUT PARAMETER
C
C     N       THE LENGTH OF THE SEQUENCE TO BE TRANSFORMED
C
C     OUTPUT PARAMETER
C
C     WSAVE   A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 4*N+15
C             THE SAME WORK ARRAY CAN BE USED FOR BOTH CFFTF AND CFFTB
C             AS LONG AS N REMAINS UNCHANGED. DIFFERENT WSAVE ARRAYS
C             ARE REQUIRED FOR DIFFERENT VALUES OF N. THE CONTENTS OF
C             WSAVE MUST NOT BE CHANGED BETWEEN CALLS OF CFFTF OR CFFTB.
C
      SUBROUTINE CFFTI (N,WSAVE)
      IMPLICIT NONE
      INTEGER N,IW1,IW2
      DOUBLE PRECISION WSAVE
      DIMENSION       WSAVE(*)
C
      IF (N .EQ. 1) RETURN
      IW1 = N+N+1
      IW2 = IW1+N+N
      CALL CFFTI1 (N,WSAVE(IW1),WSAVE(IW2))
      RETURN
      END

C****************************************************************************
C     SUBROUTINE CFFTF(N,C,WSAVE)
C
C     SUBROUTINE CFFTF COMPUTES THE FORWARD COMPLEX DISCRETE FOURIER
C     TRANSFORM (THE FOURIER ANALYSIS). EQUIVALENTLY , CFFTF COMPUTES
C     THE FOURIER COEFFICIENTS OF A COMPLEX PERIODIC SEQUENCE.
C     THE TRANSFORM IS DEFINED BELOW AT OUTPUT PARAMETER C.
C
C     THE TRANSFORM IS NOT NORMALIZED. TO OBTAIN A NORMALIZED TRANSFORM
C     THE OUTPUT MUST BE DIVIDED BY N. OTHERWISE A CALL OF CFFTF
C     FOLLOWED BY A CALL OF CFFTB WILL MULTIPLY THE SEQUENCE BY N.
C
C     THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE CFFTF MUST BE
C     INITIALIZED BY CALLING SUBROUTINE CFFTI(N,WSAVE).
C
C     INPUT PARAMETERS
C
C
C     N      THE LENGTH OF THE COMPLEX SEQUENCE C. THE METHOD IS
C            MORE EFFICIENT WHEN N IS THE PRODUCT OF SMALL PRIMES. N
C
C     C      A COMPLEX ARRAY OF LENGTH N WHICH CONTAINS THE SEQUENCE
C
C     WSAVE   A REAL WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 4N+15
C             IN THE PROGRAM THAT CALLS CFFTF. THE WSAVE ARRAY MUST BE
C             INITIALIZED BY CALLING SUBROUTINE CFFTI(N,WSAVE) AND A
C             DIFFERENT WSAVE ARRAY MUST BE USED FOR EACH DIFFERENT
C             VALUE OF N. THIS INITIALIZATION DOES NOT HAVE TO BE
C             REPEATED SO LONG AS N REMAINS UNCHANGED THUS SUBSEQUENT
C             TRANSFORMS CAN BE OBTAINED FASTER THAN THE FIRST.
C             THE SAME WSAVE ARRAY CAN BE USED BY CFFTF AND CFFTB.
C
C     OUTPUT PARAMETERS
C
C     C      FOR J=1,...,N
C
C              C(J)=THE SUM FROM K=1,...,N OF
C
C                    C(K)*EXP(-I*(J-1)*(K-1)*2*PI/N)
C
C                       WHERE I=SQRT(-1)
C
C     WSAVE   CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT BE
C             DESTROYED BETWEEN CALLS OF SUBROUTINE CFFTF OR CFFTB
C
      SUBROUTINE CFFTF (N,C,WSAVE)
      IMPLICIT NONE
      INTEGER N,IW1,IW2
      DOUBLE PRECISION C,WSAVE
      DIMENSION       C(*)       ,WSAVE(*)
C
      IF (N .EQ. 1) RETURN
      IW1 = N+N+1
      IW2 = IW1+N+N
      CALL CFFTF1 (N,C,WSAVE,WSAVE(IW1),WSAVE(IW2))
      RETURN
      END

C****************************************************************************
C     SUBROUTINE CFFTB(N,C,WSAVE)
C
C     SUBROUTINE CFFTB COMPUTES THE BACKWARD COMPLEX DISCRETE FOURIER
C     TRANSFORM (THE FOURIER SYNTHESIS). EQUIVALENTLY , CFFTB COMPUTES
C     A COMPLEX PERIODIC SEQUENCE FROM ITS FOURIER COEFFICIENTS.
C     THE TRANSFORM IS DEFINED BELOW AT OUTPUT PARAMETER C.
C
C     A CALL OF CFFTF FOLLOWED BY A CALL OF CFFTB WILL MULTIPLY THE
C     SEQUENCE BY N.
C
C     THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE CFFTB MUST BE
C     INITIALIZED BY CALLING SUBROUTINE CFFTI(N,WSAVE).
C
C     INPUT PARAMETERS
C
C
C     N      THE LENGTH OF THE COMPLEX SEQUENCE C. THE METHOD IS
C            MORE EFFICIENT WHEN N IS THE PRODUCT OF SMALL PRIMES.
C
C     C      A COMPLEX ARRAY OF LENGTH N WHICH CONTAINS THE SEQUENCE
C
C     WSAVE   A REAL WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 4N+15
C             IN THE PROGRAM THAT CALLS CFFTB. THE WSAVE ARRAY MUST BE
C             INITIALIZED BY CALLING SUBROUTINE CFFTI(N,WSAVE) AND A
C             DIFFERENT WSAVE ARRAY MUST BE USED FOR EACH DIFFERENT
C             VALUE OF N. THIS INITIALIZATION DOES NOT HAVE TO BE
C             REPEATED SO LONG AS N REMAINS UNCHANGED THUS SUBSEQUENT
C             TRANSFORMS CAN BE OBTAINED FASTER THAN THE FIRST.
C             THE SAME WSAVE ARRAY CAN BE USED BY CFFTF AND CFFTB.
C
C     OUTPUT PARAMETERS
C
C     C      FOR J=1,...,N
C
C              C(J)=THE SUM FROM K=1,...,N OF
C
C                    C(K)*EXP(I*(J-1)*(K-1)*2*PI/N)
C
C                       WHERE I=SQRT(-1)
C
C     WSAVE   CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT BE
C             DESTROYED BETWEEN CALLS OF SUBROUTINE CFFTF OR CFFTB
C
      SUBROUTINE CFFTB (N,C,WSAVE)
      IMPLICIT NONE
      INTEGER N,IW1,IW2
      DOUBLE PRECISION C,WSAVE
      DIMENSION       C(*)       ,WSAVE(*)
C
      IF (N .EQ. 1) RETURN
      IW1 = N+N+1
      IW2 = IW1+N+N
      CALL CFFTB1 (N,C,WSAVE,WSAVE(IW1),WSAVE(IW2))
      RETURN
      END

C****************************************************************************
      SUBROUTINE CFFTI1 (N,WA,WIFAC)
      IMPLICIT NONE
      INTEGER N,NTRYH,NL,NF,J,NTRY,NQ,NR,IB,I,L1,K1,IP,LD,L2,
     &  IDO,IDOT,IPM,I1,II
      DOUBLE PRECISION WA,TPI,ARGH,FI,ARGLD,ARG,WIFAC
      DIMENSION       WA(*)      ,WIFAC(*)    ,NTRYH(4)
      DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/3,4,2,5/
      NL = N
      NF = 0
      J = 0
  101 J = J+1
      IF (J-4) 102,102,103
  102 NTRY = NTRYH(J)
      GO TO 104
  103 NTRY = NTRY+2
  104 NQ = NL/NTRY
      NR = NL-NTRY*NQ
      IF (NR) 101,105,101
  105 NF = NF+1
      WIFAC(NF+2) = DBLE(NTRY)
      NL = NQ
      IF (NTRY .NE. 2) GO TO 107
      IF (NF .EQ. 1) GO TO 107
      DO 106 I=2,NF
         IB = NF-I+2
         WIFAC(IB+2) = WIFAC(IB+1)
  106 CONTINUE
      WIFAC(3) = 2.D0
  107 IF (NL .NE. 1) GO TO 104
      WIFAC(1) = DBLE(N)
      WIFAC(2) = DBLE(NF)
      TPI = 2.D0*(4.D0*ATAN(1.D0))
      ARGH = TPI/DBLE(N)
      I = 2
      L1 = 1
      DO 110 K1=1,NF
         IP = NINT(WIFAC(K1+2))
         LD = 0
         L2 = L1*IP
         IDO = N/L2
         IDOT = IDO+IDO+2
         IPM = IP-1
         DO 109 J=1,IPM
            I1 = I
            WA(I-1) = 1.D0
            WA(I) = 0.D0
            LD = LD+L1
            FI = 0.D0
            ARGLD = DBLE(LD)*ARGH
            DO 108 II=4,IDOT,2
               I = I+2
               FI = FI+1.D0
               ARG = FI*ARGLD
               WA(I-1) = COS(ARG)
               WA(I) = SIN(ARG)
  108       CONTINUE
            IF (IP .LE. 5) GO TO 109
            WA(I1-1) = WA(I-1)
            WA(I1) = WA(I)
  109    CONTINUE
         L1 = L2
  110 CONTINUE
      RETURN
      END

C****************************************************************************
      SUBROUTINE CFFTF1 (N,C,CH,WA,WIFAC)
      IMPLICIT NONE
      INTEGER N,NF,NA,L1,IW,K1,IP,L2,IDO,IDOT,IDL1,IX2,IX3,IX4,N2,
     &  I,NAC
      DOUBLE PRECISION C,CH,WA,WIFAC
      DIMENSION       CH(*)      ,C(*)           ,WA(*)   ,WIFAC(*)
      NF = NINT(WIFAC(2))
      NA = 0
      L1 = 1
      IW = 1
      DO 116 K1=1,NF
         IP = NINT(WIFAC(K1+2))
         L2 = IP*L1
         IDO = N/L2
         IDOT = IDO+IDO
         IDL1 = IDOT*L1
         IF (IP .NE. 4) GO TO 103
         IX2 = IW+IDOT
         IX3 = IX2+IDOT
         IF (NA .NE. 0) GO TO 101
         CALL PASSF4 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
         GO TO 102
  101    CALL PASSF4 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
  102    NA = 1-NA
         GO TO 115
  103    IF (IP .NE. 2) GO TO 106
         IF (NA .NE. 0) GO TO 104
         CALL PASSF2 (IDOT,L1,C,CH,WA(IW))
         GO TO 105
  104    CALL PASSF2 (IDOT,L1,CH,C,WA(IW))
  105    NA = 1-NA
         GO TO 115
  106    IF (IP .NE. 3) GO TO 109
         IX2 = IW+IDOT
         IF (NA .NE. 0) GO TO 107
         CALL PASSF3 (IDOT,L1,C,CH,WA(IW),WA(IX2))
         GO TO 108
  107    CALL PASSF3 (IDOT,L1,CH,C,WA(IW),WA(IX2))
  108    NA = 1-NA
         GO TO 115
  109    IF (IP .NE. 5) GO TO 112
         IX2 = IW+IDOT
         IX3 = IX2+IDOT
         IX4 = IX3+IDOT
         IF (NA .NE. 0) GO TO 110
         CALL PASSF5 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
         GO TO 111
  110    CALL PASSF5 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
  111    NA = 1-NA
         GO TO 115
  112    IF (NA .NE. 0) GO TO 113
         CALL PASSF (NAC,IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
         GO TO 114
  113    CALL PASSF (NAC,IDOT,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
  114    IF (NAC .NE. 0) NA = 1-NA
  115    L1 = L2
         IW = IW+(IP-1)*IDOT
  116 CONTINUE
      IF (NA .EQ. 0) RETURN
      N2 = N+N
      DO 117 I=1,N2
         C(I) = CH(I)
  117 CONTINUE
      RETURN
      END

C****************************************************************************
      SUBROUTINE CFFTB1 (N,C,CH,WA,WIFAC)
      IMPLICIT NONE
      INTEGER N,NF,NA,L1,IW,K1,IP,L2,IDO,IDOT,IDL1,IX2,IX3,IX4,N2,
     &  I,NAC
      DOUBLE PRECISION C,CH,WA,WIFAC
      DIMENSION       CH(*)      ,C(*)           ,WA(*)   ,WIFAC(*)
      NF = NINT(WIFAC(2))
      NA = 0
      L1 = 1
      IW = 1
      DO 116 K1=1,NF
         IP = NINT(WIFAC(K1+2))
         L2 = IP*L1
         IDO = N/L2
         IDOT = IDO+IDO
         IDL1 = IDOT*L1
         IF (IP .NE. 4) GO TO 103
         IX2 = IW+IDOT
         IX3 = IX2+IDOT
         IF (NA .NE. 0) GO TO 101
         CALL PASSB4 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
         GO TO 102
  101    CALL PASSB4 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
  102    NA = 1-NA
         GO TO 115
  103    IF (IP .NE. 2) GO TO 106
         IF (NA .NE. 0) GO TO 104
         CALL PASSB2 (IDOT,L1,C,CH,WA(IW))
         GO TO 105
  104    CALL PASSB2 (IDOT,L1,CH,C,WA(IW))
  105    NA = 1-NA
         GO TO 115
  106    IF (IP .NE. 3) GO TO 109
         IX2 = IW+IDOT
         IF (NA .NE. 0) GO TO 107
         CALL PASSB3 (IDOT,L1,C,CH,WA(IW),WA(IX2))
         GO TO 108
  107    CALL PASSB3 (IDOT,L1,CH,C,WA(IW),WA(IX2))
  108    NA = 1-NA
         GO TO 115
  109    IF (IP .NE. 5) GO TO 112
         IX2 = IW+IDOT
         IX3 = IX2+IDOT
         IX4 = IX3+IDOT
         IF (NA .NE. 0) GO TO 110
         CALL PASSB5 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
         GO TO 111
  110    CALL PASSB5 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
  111    NA = 1-NA
         GO TO 115
  112    IF (NA .NE. 0) GO TO 113
         CALL PASSB (NAC,IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
         GO TO 114
  113    CALL PASSB (NAC,IDOT,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
  114    IF (NAC .NE. 0) NA = 1-NA
  115    L1 = L2
         IW = IW+(IP-1)*IDOT
  116 CONTINUE
      IF (NA .EQ. 0) RETURN
      N2 = N+N
      DO 117 I=1,N2
         C(I) = CH(I)
  117 CONTINUE
      RETURN
      END

C****************************************************************************
      SUBROUTINE PASSB (NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
      IMPLICIT NONE
      INTEGER NAC,IDO,IP,L1,IDL1,IDOT,IPP2,IPPH,IDP,J,K,I,IDL,LC,L,
     &  IK,IDLJ,JC,INC,IDIJ,IDJ
      DOUBLE PRECISION CC,C1,C2,CH,CH2,WA,WAR,WAI
      DIMENSION       CH(IDO,L1,IP)          ,CC(IDO,IP,L1)      ,
     1         C1(IDO,L1,IP)          ,WA(*)  ,C2(IDL1,IP),
     2         CH2(IDL1,IP)
      IDOT = IDO/2
      IPP2 = IP+2
      IPPH = (IP+1)/2
      IDP = IP*IDO
C
      IF (IDO .LT. L1) GO TO 106
      DO 103 J=2,IPPH
         JC = IPP2-J
         DO 102 K=1,L1
            DO 101 I=1,IDO
               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
  101       CONTINUE
  102    CONTINUE
  103 CONTINUE
      DO 105 K=1,L1
         DO 104 I=1,IDO
            CH(I,K,1) = CC(I,1,K)
  104    CONTINUE
  105 CONTINUE
      GO TO 112
  106 DO 109 J=2,IPPH
         JC = IPP2-J
         DO 108 I=1,IDO
            DO 107 K=1,L1
               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
  107       CONTINUE
  108    CONTINUE
  109 CONTINUE
      DO 111 I=1,IDO
         DO 110 K=1,L1
            CH(I,K,1) = CC(I,1,K)
  110    CONTINUE
  111 CONTINUE
  112 IDL = 2-IDO
      INC = 0
      DO 116 L=2,IPPH
         LC = IPP2-L
         IDL = IDL+IDO
         DO 113 IK=1,IDL1
            C2(IK,L) = CH2(IK,1)+WA(IDL-1)*CH2(IK,2)
            C2(IK,LC) = WA(IDL)*CH2(IK,IP)
  113    CONTINUE
         IDLJ = IDL
         INC = INC+IDO
         DO 115 J=3,IPPH
            JC = IPP2-J
            IDLJ = IDLJ+INC
            IF (IDLJ .GT. IDP) IDLJ = IDLJ-IDP
            WAR = WA(IDLJ-1)
            WAI = WA(IDLJ)
            DO 114 IK=1,IDL1
               C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J)
               C2(IK,LC) = C2(IK,LC)+WAI*CH2(IK,JC)
  114       CONTINUE
  115    CONTINUE
  116 CONTINUE
      DO 118 J=2,IPPH
         DO 117 IK=1,IDL1
            CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
  117    CONTINUE
  118 CONTINUE
      DO 120 J=2,IPPH
         JC = IPP2-J
         DO 119 IK=2,IDL1,2
            CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC)
            CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC)
            CH2(IK,J) = C2(IK,J)+C2(IK-1,JC)
            CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC)
  119    CONTINUE
  120 CONTINUE
      NAC = 1
      IF (IDO .EQ. 2) RETURN
      NAC = 0
      DO 121 IK=1,IDL1
         C2(IK,1) = CH2(IK,1)
  121 CONTINUE
      DO 123 J=2,IP
         DO 122 K=1,L1
            C1(1,K,J) = CH(1,K,J)
            C1(2,K,J) = CH(2,K,J)
  122    CONTINUE
  123 CONTINUE
      IF (IDOT .GT. L1) GO TO 127
      IDIJ = 0
      DO 126 J=2,IP
         IDIJ = IDIJ+2
         DO 125 I=4,IDO,2
            IDIJ = IDIJ+2
            DO 124 K=1,L1
               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
  124       CONTINUE
  125    CONTINUE
  126 CONTINUE
      RETURN
  127 IDJ = 2-IDO
      DO 130 J=2,IP
         IDJ = IDJ+IDO
         DO 129 K=1,L1
            IDIJ = IDJ
            DO 128 I=4,IDO,2
               IDIJ = IDIJ+2
               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
  128       CONTINUE
  129    CONTINUE
  130 CONTINUE
      RETURN
      END

C****************************************************************************
      SUBROUTINE PASSB2 (IDO,L1,CC,CH,WA1)
      IMPLICIT NONE
      INTEGER IDO,L1,K,I
      DOUBLE PRECISION CC,CH,WA1,TR2,TI2
      DIMENSION       CC(IDO,2,L1)           ,CH(IDO,L1,2)        ,
     1         WA1(1)
      IF (IDO .GT. 2) GO TO 102
      DO 101 K=1,L1
         CH(1,K,1) = CC(1,1,K)+CC(1,2,K)
         CH(1,K,2) = CC(1,1,K)-CC(1,2,K)
         CH(2,K,1) = CC(2,1,K)+CC(2,2,K)
         CH(2,K,2) = CC(2,1,K)-CC(2,2,K)
  101 CONTINUE
      RETURN
  102 DO 104 K=1,L1
         DO 103 I=2,IDO,2
            CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K)
            TR2 = CC(I-1,1,K)-CC(I-1,2,K)
            CH(I,K,1) = CC(I,1,K)+CC(I,2,K)
            TI2 = CC(I,1,K)-CC(I,2,K)
            CH(I,K,2) = WA1(I-1)*TI2+WA1(I)*TR2
            CH(I-1,K,2) = WA1(I-1)*TR2-WA1(I)*TI2
  103    CONTINUE
  104 CONTINUE
      RETURN
      END

C****************************************************************************
      SUBROUTINE PASSB3 (IDO,L1,CC,CH,WA1,WA2)
      IMPLICIT NONE
      INTEGER IDO,L1,K,I
      DOUBLE PRECISION CC,CH,WA1,WA2,TAUR,TAUI,TR2,CR2,TI2,CI2,CR3,CI3,
     &  DR2,DR3,DI2,DI3
      DIMENSION       CC(IDO,3,L1)           ,CH(IDO,L1,3)        ,
     1         WA1(*)     ,WA2(*)
      DATA TAUR,TAUI /-.5D0,.866025403784439D0/
      IF (IDO .NE. 2) GO TO 102
      DO 101 K=1,L1
         TR2 = CC(1,2,K)+CC(1,3,K)
         CR2 = CC(1,1,K)+TAUR*TR2
         CH(1,K,1) = CC(1,1,K)+TR2
         TI2 = CC(2,2,K)+CC(2,3,K)
         CI2 = CC(2,1,K)+TAUR*TI2
         CH(2,K,1) = CC(2,1,K)+TI2
         CR3 = TAUI*(CC(1,2,K)-CC(1,3,K))
         CI3 = TAUI*(CC(2,2,K)-CC(2,3,K))
         CH(1,K,2) = CR2-CI3
         CH(1,K,3) = CR2+CI3
         CH(2,K,2) = CI2+CR3
         CH(2,K,3) = CI2-CR3
  101 CONTINUE
      RETURN
  102 DO 104 K=1,L1
         DO 103 I=2,IDO,2
            TR2 = CC(I-1,2,K)+CC(I-1,3,K)
            CR2 = CC(I-1,1,K)+TAUR*TR2
            CH(I-1,K,1) = CC(I-1,1,K)+TR2
            TI2 = CC(I,2,K)+CC(I,3,K)
            CI2 = CC(I,1,K)+TAUR*TI2
            CH(I,K,1) = CC(I,1,K)+TI2
            CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K))
            CI3 = TAUI*(CC(I,2,K)-CC(I,3,K))
            DR2 = CR2-CI3
            DR3 = CR2+CI3
            DI2 = CI2+CR3
            DI3 = CI2-CR3
            CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2
            CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2
            CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3
            CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3
  103    CONTINUE
  104 CONTINUE
      RETURN
      END

C****************************************************************************
      SUBROUTINE PASSB4 (IDO,L1,CC,CH,WA1,WA2,WA3)
      IMPLICIT NONE
      INTEGER IDO,L1,K,I
      DOUBLE PRECISION CC,CH,WA1,WA2,WA3,TI1,TI2,TI3,TI4,TR1,TR2,TR3,
     &  TR4,CR3,CR2,CI3,CR4,CI4,CI2
      DIMENSION       CC(IDO,4,L1)           ,CH(IDO,L1,4)        ,
     1         WA1(*)     ,WA2(*)     ,WA3(*)
      IF (IDO .NE. 2) GO TO 102
      DO 101 K=1,L1
         TI1 = CC(2,1,K)-CC(2,3,K)
         TI2 = CC(2,1,K)+CC(2,3,K)
         TR4 = CC(2,4,K)-CC(2,2,K)
         TI3 = CC(2,2,K)+CC(2,4,K)
         TR1 = CC(1,1,K)-CC(1,3,K)
         TR2 = CC(1,1,K)+CC(1,3,K)
         TI4 = CC(1,2,K)-CC(1,4,K)
         TR3 = CC(1,2,K)+CC(1,4,K)
         CH(1,K,1) = TR2+TR3
         CH(1,K,3) = TR2-TR3
         CH(2,K,1) = TI2+TI3
         CH(2,K,3) = TI2-TI3
         CH(1,K,2) = TR1+TR4
         CH(1,K,4) = TR1-TR4
         CH(2,K,2) = TI1+TI4
         CH(2,K,4) = TI1-TI4
  101 CONTINUE
      RETURN
  102 DO 104 K=1,L1
         DO 103 I=2,IDO,2
            TI1 = CC(I,1,K)-CC(I,3,K)
            TI2 = CC(I,1,K)+CC(I,3,K)
            TI3 = CC(I,2,K)+CC(I,4,K)
            TR4 = CC(I,4,K)-CC(I,2,K)
            TR1 = CC(I-1,1,K)-CC(I-1,3,K)
            TR2 = CC(I-1,1,K)+CC(I-1,3,K)
            TI4 = CC(I-1,2,K)-CC(I-1,4,K)
            TR3 = CC(I-1,2,K)+CC(I-1,4,K)
            CH(I-1,K,1) = TR2+TR3
            CR3 = TR2-TR3
            CH(I,K,1) = TI2+TI3
            CI3 = TI2-TI3
            CR2 = TR1+TR4
            CR4 = TR1-TR4
            CI2 = TI1+TI4
            CI4 = TI1-TI4
            CH(I-1,K,2) = WA1(I-1)*CR2-WA1(I)*CI2
            CH(I,K,2) = WA1(I-1)*CI2+WA1(I)*CR2
            CH(I-1,K,3) = WA2(I-1)*CR3-WA2(I)*CI3
            CH(I,K,3) = WA2(I-1)*CI3+WA2(I)*CR3
            CH(I-1,K,4) = WA3(I-1)*CR4-WA3(I)*CI4
            CH(I,K,4) = WA3(I-1)*CI4+WA3(I)*CR4
  103    CONTINUE
  104 CONTINUE
      RETURN
      END

C****************************************************************************
      SUBROUTINE PASSB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
      IMPLICIT NONE
      INTEGER IDO,L1,K,I
      DOUBLE PRECISION CC,CH,WA1,WA2,WA3,WA4,TR11,TI11,TR12,TI12,
     &  TI2,TI3,TI4,TI5,TR2,TR3,TR4,TR5,CR2,CI2,CR3,CI3,CR4,CI4,CR5,CI5,
     &  DR3,DI3,DR4,DI4,DR5,DI5,DR2,DI2
      DIMENSION       CC(IDO,5,L1)           ,CH(IDO,L1,5)        ,
     1         WA1(*)     ,WA2(*)     ,WA3(*)        ,WA4(*)
      DATA TR11,TI11,TR12,TI12 /.309016994374947D0,.951056516295154D0,
     1-.809016994374947D0,.587785252292473D0/
      IF (IDO .NE. 2) GO TO 102
      DO 101 K=1,L1
         TI5 = CC(2,2,K)-CC(2,5,K)
         TI2 = CC(2,2,K)+CC(2,5,K)
         TI4 = CC(2,3,K)-CC(2,4,K)
         TI3 = CC(2,3,K)+CC(2,4,K)
         TR5 = CC(1,2,K)-CC(1,5,K)
         TR2 = CC(1,2,K)+CC(1,5,K)
         TR4 = CC(1,3,K)-CC(1,4,K)
         TR3 = CC(1,3,K)+CC(1,4,K)
         CH(1,K,1) = CC(1,1,K)+TR2+TR3
         CH(2,K,1) = CC(2,1,K)+TI2+TI3
         CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
         CI2 = CC(2,1,K)+TR11*TI2+TR12*TI3
         CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
         CI3 = CC(2,1,K)+TR12*TI2+TR11*TI3
         CR5 = TI11*TR5+TI12*TR4
         CI5 = TI11*TI5+TI12*TI4
         CR4 = TI12*TR5-TI11*TR4
         CI4 = TI12*TI5-TI11*TI4
         CH(1,K,2) = CR2-CI5
         CH(1,K,5) = CR2+CI5
         CH(2,K,2) = CI2+CR5
         CH(2,K,3) = CI3+CR4
         CH(1,K,3) = CR3-CI4
         CH(1,K,4) = CR3+CI4
         CH(2,K,4) = CI3-CR4
         CH(2,K,5) = CI2-CR5
  101 CONTINUE
      RETURN
  102 DO 104 K=1,L1
         DO 103 I=2,IDO,2
            TI5 = CC(I,2,K)-CC(I,5,K)
            TI2 = CC(I,2,K)+CC(I,5,K)
            TI4 = CC(I,3,K)-CC(I,4,K)
            TI3 = CC(I,3,K)+CC(I,4,K)
            TR5 = CC(I-1,2,K)-CC(I-1,5,K)
            TR2 = CC(I-1,2,K)+CC(I-1,5,K)
            TR4 = CC(I-1,3,K)-CC(I-1,4,K)
            TR3 = CC(I-1,3,K)+CC(I-1,4,K)
            CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
            CH(I,K,1) = CC(I,1,K)+TI2+TI3
            CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
            CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
            CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
            CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
            CR5 = TI11*TR5+TI12*TR4
            CI5 = TI11*TI5+TI12*TI4
            CR4 = TI12*TR5-TI11*TR4
            CI4 = TI12*TI5-TI11*TI4
            DR3 = CR3-CI4
            DR4 = CR3+CI4
            DI3 = CI3+CR4
            DI4 = CI3-CR4
            DR5 = CR2+CI5
            DR2 = CR2-CI5
            DI5 = CI2-CR5
            DI2 = CI2+CR5
            CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2
            CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2
            CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3
            CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3
            CH(I-1,K,4) = WA3(I-1)*DR4-WA3(I)*DI4
            CH(I,K,4) = WA3(I-1)*DI4+WA3(I)*DR4
            CH(I-1,K,5) = WA4(I-1)*DR5-WA4(I)*DI5
            CH(I,K,5) = WA4(I-1)*DI5+WA4(I)*DR5
  103    CONTINUE
  104 CONTINUE
      RETURN
      END

C****************************************************************************
      SUBROUTINE PASSF (NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
      IMPLICIT NONE
      INTEGER NAC,IDO,IP,L1,IDL1,IDOT,IPP2,IPPH,IDP,J,JC,K,I,IDL,
     &  INC,L,LC,IK,IDLJ,IDIJ,IDJ
      DOUBLE PRECISION CC,C1,C2,CH,CH2,WA,WAR,WAI
      DIMENSION       CH(IDO,L1,IP)          ,CC(IDO,IP,L1)      ,
     1         C1(IDO,L1,IP)          ,WA(*)  ,C2(IDL1,IP),
     2         CH2(IDL1,IP)
      IDOT = IDO/2
      IPP2 = IP+2
      IPPH = (IP+1)/2
      IDP = IP*IDO
C
      IF (IDO .LT. L1) GO TO 106
      DO 103 J=2,IPPH
         JC = IPP2-J
         DO 102 K=1,L1
            DO 101 I=1,IDO
               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
  101       CONTINUE
  102    CONTINUE
  103 CONTINUE
      DO 105 K=1,L1
         DO 104 I=1,IDO
            CH(I,K,1) = CC(I,1,K)
  104    CONTINUE
  105 CONTINUE
      GO TO 112
  106 DO 109 J=2,IPPH
         JC = IPP2-J
         DO 108 I=1,IDO
            DO 107 K=1,L1
               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
  107       CONTINUE
  108    CONTINUE
  109 CONTINUE
      DO 111 I=1,IDO
         DO 110 K=1,L1
            CH(I,K,1) = CC(I,1,K)
  110    CONTINUE
  111 CONTINUE
  112 IDL = 2-IDO
      INC = 0
      DO 116 L=2,IPPH
         LC = IPP2-L
         IDL = IDL+IDO
         DO 113 IK=1,IDL1
            C2(IK,L) = CH2(IK,1)+WA(IDL-1)*CH2(IK,2)
            C2(IK,LC) = -WA(IDL)*CH2(IK,IP)
  113    CONTINUE
         IDLJ = IDL
         INC = INC+IDO
         DO 115 J=3,IPPH
            JC = IPP2-J
            IDLJ = IDLJ+INC
            IF (IDLJ .GT. IDP) IDLJ = IDLJ-IDP
            WAR = WA(IDLJ-1)
            WAI = WA(IDLJ)
            DO 114 IK=1,IDL1
               C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J)
               C2(IK,LC) = C2(IK,LC)-WAI*CH2(IK,JC)
  114       CONTINUE
  115    CONTINUE
  116 CONTINUE
      DO 118 J=2,IPPH
         DO 117 IK=1,IDL1
            CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
  117    CONTINUE
  118 CONTINUE
      DO 120 J=2,IPPH
         JC = IPP2-J
         DO 119 IK=2,IDL1,2
            CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC)
            CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC)
            CH2(IK,J) = C2(IK,J)+C2(IK-1,JC)
            CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC)
  119    CONTINUE
  120 CONTINUE
      NAC = 1
      IF (IDO .EQ. 2) RETURN
      NAC = 0
      DO 121 IK=1,IDL1
         C2(IK,1) = CH2(IK,1)
  121 CONTINUE
      DO 123 J=2,IP
         DO 122 K=1,L1
            C1(1,K,J) = CH(1,K,J)
            C1(2,K,J) = CH(2,K,J)
  122    CONTINUE
  123 CONTINUE
      IF (IDOT .GT. L1) GO TO 127
      IDIJ = 0
      DO 126 J=2,IP
         IDIJ = IDIJ+2
         DO 125 I=4,IDO,2
            IDIJ = IDIJ+2
            DO 124 K=1,L1
               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J)
               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J)
  124       CONTINUE
  125    CONTINUE
  126 CONTINUE
      RETURN
  127 IDJ = 2-IDO
      DO 130 J=2,IP
         IDJ = IDJ+IDO
         DO 129 K=1,L1
            IDIJ = IDJ
            DO 128 I=4,IDO,2
               IDIJ = IDIJ+2
               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J)
               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J)
  128       CONTINUE
  129    CONTINUE
  130 CONTINUE
      RETURN
      END

C****************************************************************************
      SUBROUTINE PASSF2 (IDO,L1,CC,CH,WA1)
      IMPLICIT NONE
      INTEGER IDO,L1,K,I
      DOUBLE PRECISION CC,CH,WA1,TR2,TI2
      DIMENSION       CC(IDO,2,L1)           ,CH(IDO,L1,2)        ,
     1         WA1(*)
      IF (IDO .GT. 2) GO TO 102
      DO 101 K=1,L1
         CH(1,K,1) = CC(1,1,K)+CC(1,2,K)
         CH(1,K,2) = CC(1,1,K)-CC(1,2,K)
         CH(2,K,1) = CC(2,1,K)+CC(2,2,K)
         CH(2,K,2) = CC(2,1,K)-CC(2,2,K)
  101 CONTINUE
      RETURN
  102 DO 104 K=1,L1
         DO 103 I=2,IDO,2
            CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K)
            TR2 = CC(I-1,1,K)-CC(I-1,2,K)
            CH(I,K,1) = CC(I,1,K)+CC(I,2,K)
            TI2 = CC(I,1,K)-CC(I,2,K)
            CH(I,K,2) = WA1(I-1)*TI2-WA1(I)*TR2
            CH(I-1,K,2) = WA1(I-1)*TR2+WA1(I)*TI2
  103    CONTINUE
  104 CONTINUE
      RETURN
      END

C****************************************************************************
      SUBROUTINE PASSF3 (IDO,L1,CC,CH,WA1,WA2)
      IMPLICIT NONE
      INTEGER IDO,L1,K,I
      DOUBLE PRECISION CC,CH,WA1,WA2,TAUR,TAUI,TR2,CR2,TI2,CI2,CR3,CI3,
     &  DR2,DR3,DI2,DI3
      DIMENSION       CC(IDO,3,L1)           ,CH(IDO,L1,3)        ,
     1         WA1(*)     ,WA2(*)
      DATA TAUR,TAUI /-.5D0,-.866025403784439D0/
      IF (IDO .NE. 2) GO TO 102
      DO 101 K=1,L1
         TR2 = CC(1,2,K)+CC(1,3,K)
         CR2 = CC(1,1,K)+TAUR*TR2
         CH(1,K,1) = CC(1,1,K)+TR2
         TI2 = CC(2,2,K)+CC(2,3,K)
         CI2 = CC(2,1,K)+TAUR*TI2
         CH(2,K,1) = CC(2,1,K)+TI2
         CR3 = TAUI*(CC(1,2,K)-CC(1,3,K))
         CI3 = TAUI*(CC(2,2,K)-CC(2,3,K))
         CH(1,K,2) = CR2-CI3
         CH(1,K,3) = CR2+CI3
         CH(2,K,2) = CI2+CR3
         CH(2,K,3) = CI2-CR3
  101 CONTINUE
      RETURN
  102 DO 104 K=1,L1
         DO 103 I=2,IDO,2
            TR2 = CC(I-1,2,K)+CC(I-1,3,K)
            CR2 = CC(I-1,1,K)+TAUR*TR2
            CH(I-1,K,1) = CC(I-1,1,K)+TR2
            TI2 = CC(I,2,K)+CC(I,3,K)
            CI2 = CC(I,1,K)+TAUR*TI2
            CH(I,K,1) = CC(I,1,K)+TI2
            CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K))
            CI3 = TAUI*(CC(I,2,K)-CC(I,3,K))
            DR2 = CR2-CI3
            DR3 = CR2+CI3
            DI2 = CI2+CR3
            DI3 = CI2-CR3
            CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2
            CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2
            CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3
            CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3
  103    CONTINUE
  104 CONTINUE
      RETURN
      END

C****************************************************************************
      SUBROUTINE PASSF4 (IDO,L1,CC,CH,WA1,WA2,WA3)
      IMPLICIT NONE
      INTEGER IDO,L1,K,I
      DOUBLE PRECISION CC,CH,WA1,WA2,WA3,TI1,TI2,TI3,TI4,TR1,TR2,TR3,
     &  TR4,CR2,CR3,CR4,CI2,CI3,CI4
      DIMENSION       CC(IDO,4,L1)           ,CH(IDO,L1,4)        ,
     1         WA1(*)     ,WA2(*)     ,WA3(*)
      IF (IDO .NE. 2) GO TO 102
      DO 101 K=1,L1
         TI1 = CC(2,1,K)-CC(2,3,K)
         TI2 = CC(2,1,K)+CC(2,3,K)
         TR4 = CC(2,2,K)-CC(2,4,K)
         TI3 = CC(2,2,K)+CC(2,4,K)
         TR1 = CC(1,1,K)-CC(1,3,K)
         TR2 = CC(1,1,K)+CC(1,3,K)
         TI4 = CC(1,4,K)-CC(1,2,K)
         TR3 = CC(1,2,K)+CC(1,4,K)
         CH(1,K,1) = TR2+TR3
         CH(1,K,3) = TR2-TR3
         CH(2,K,1) = TI2+TI3
         CH(2,K,3) = TI2-TI3
         CH(1,K,2) = TR1+TR4
         CH(1,K,4) = TR1-TR4
         CH(2,K,2) = TI1+TI4
         CH(2,K,4) = TI1-TI4
  101 CONTINUE
      RETURN
  102 DO 104 K=1,L1
         DO 103 I=2,IDO,2
            TI1 = CC(I,1,K)-CC(I,3,K)
            TI2 = CC(I,1,K)+CC(I,3,K)
            TI3 = CC(I,2,K)+CC(I,4,K)
            TR4 = CC(I,2,K)-CC(I,4,K)
            TR1 = CC(I-1,1,K)-CC(I-1,3,K)
            TR2 = CC(I-1,1,K)+CC(I-1,3,K)
            TI4 = CC(I-1,4,K)-CC(I-1,2,K)
            TR3 = CC(I-1,2,K)+CC(I-1,4,K)
            CH(I-1,K,1) = TR2+TR3
            CR3 = TR2-TR3
            CH(I,K,1) = TI2+TI3
            CI3 = TI2-TI3
            CR2 = TR1+TR4
            CR4 = TR1-TR4
            CI2 = TI1+TI4
            CI4 = TI1-TI4
            CH(I-1,K,2) = WA1(I-1)*CR2+WA1(I)*CI2
            CH(I,K,2) = WA1(I-1)*CI2-WA1(I)*CR2
            CH(I-1,K,3) = WA2(I-1)*CR3+WA2(I)*CI3
            CH(I,K,3) = WA2(I-1)*CI3-WA2(I)*CR3
            CH(I-1,K,4) = WA3(I-1)*CR4+WA3(I)*CI4
            CH(I,K,4) = WA3(I-1)*CI4-WA3(I)*CR4
  103    CONTINUE
  104 CONTINUE
      RETURN
      END

C****************************************************************************
      SUBROUTINE PASSF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
      IMPLICIT NONE
      INTEGER IDO,L1,K,I
      DOUBLE PRECISION CC,CH,WA1,WA2,WA3,WA4,TR11,TI11,TR12,TI12,
     &  TI2,TI3,TI4,TI5,TR2,TR3,TR4,TR5,CI2,CI3,CI4,CI5,CR2,CR3,CR4,
     &  CR5,DI2,DI3,DI4,DI5,DR2,DR3,DR4,DR5
      DIMENSION       CC(IDO,5,L1)           ,CH(IDO,L1,5)        ,
     1         WA1(*)     ,WA2(*)     ,WA3(*)        ,WA4(*)
      DATA TR11,TI11,TR12,TI12 /.309016994374947D0,-.951056516295154D0,
     1-.809016994374947D0,-.587785252292473D0/
      IF (IDO .NE. 2) GO TO 102
      DO 101 K=1,L1
         TI5 = CC(2,2,K)-CC(2,5,K)
         TI2 = CC(2,2,K)+CC(2,5,K)
         TI4 = CC(2,3,K)-CC(2,4,K)
         TI3 = CC(2,3,K)+CC(2,4,K)
         TR5 = CC(1,2,K)-CC(1,5,K)
         TR2 = CC(1,2,K)+CC(1,5,K)
         TR4 = CC(1,3,K)-CC(1,4,K)
         TR3 = CC(1,3,K)+CC(1,4,K)
         CH(1,K,1) = CC(1,1,K)+TR2+TR3
         CH(2,K,1) = CC(2,1,K)+TI2+TI3
         CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
         CI2 = CC(2,1,K)+TR11*TI2+TR12*TI3
         CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
         CI3 = CC(2,1,K)+TR12*TI2+TR11*TI3
         CR5 = TI11*TR5+TI12*TR4
         CI5 = TI11*TI5+TI12*TI4
         CR4 = TI12*TR5-TI11*TR4
         CI4 = TI12*TI5-TI11*TI4
         CH(1,K,2) = CR2-CI5
         CH(1,K,5) = CR2+CI5
         CH(2,K,2) = CI2+CR5
         CH(2,K,3) = CI3+CR4
         CH(1,K,3) = CR3-CI4
         CH(1,K,4) = CR3+CI4
         CH(2,K,4) = CI3-CR4
         CH(2,K,5) = CI2-CR5
  101 CONTINUE
      RETURN
  102 DO 104 K=1,L1
         DO 103 I=2,IDO,2
            TI5 = CC(I,2,K)-CC(I,5,K)
            TI2 = CC(I,2,K)+CC(I,5,K)
            TI4 = CC(I,3,K)-CC(I,4,K)
            TI3 = CC(I,3,K)+CC(I,4,K)
            TR5 = CC(I-1,2,K)-CC(I-1,5,K)
            TR2 = CC(I-1,2,K)+CC(I-1,5,K)
            TR4 = CC(I-1,3,K)-CC(I-1,4,K)
            TR3 = CC(I-1,3,K)+CC(I-1,4,K)
            CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
            CH(I,K,1) = CC(I,1,K)+TI2+TI3
            CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
            CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
            CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
            CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
            CR5 = TI11*TR5+TI12*TR4
            CI5 = TI11*TI5+TI12*TI4
            CR4 = TI12*TR5-TI11*TR4
            CI4 = TI12*TI5-TI11*TI4
            DR3 = CR3-CI4
            DR4 = CR3+CI4
            DI3 = CI3+CR4
            DI4 = CI3-CR4
            DR5 = CR2+CI5
            DR2 = CR2-CI5
            DI5 = CI2-CR5
            DI2 = CI2+CR5
            CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2
            CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2
            CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3
            CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3
            CH(I-1,K,4) = WA3(I-1)*DR4+WA3(I)*DI4
            CH(I,K,4) = WA3(I-1)*DI4-WA3(I)*DR4
            CH(I-1,K,5) = WA4(I-1)*DR5+WA4(I)*DI5
            CH(I,K,5) = WA4(I-1)*DI5-WA4(I)*DR5
  103    CONTINUE
  104 CONTINUE
      RETURN
      END


text:
c
c     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine newexchg(n_spin_pola,sp,
     >                  rhoup,rhodn,
     >                  vx,enxc,vxout,excout,ro3,dz,
     >                  r_mesh,jmt,iexch)
c     ================================================================
c
      implicit   none
c
      integer    n_spin_pola
      integer    jmt
c
      integer    iexch
      integer    ir
c
      real*8     sp
      real*8     rhoup(jmt)
      real*8     rhodn(jmt)
      real*8     vx(jmt)
      real*8     enxc(jmt)
      real*8     vxout,excout,ro3,dz
      real*8     r_mesh(jmt)
c
      real*8     rhot
      real*8     alpha2
      real*8     dzr
      real*8     ro3r
      real*8     zero
      real*8     one
      real*8     three
      real*8     third
c
      parameter  (zero=0.0d0)
      parameter  (one=1.0d0)
      parameter  (three=3.0d0)
      parameter  (third=one/three)
c
c     ================================================================
c     calculate exchange-correlation potential and density...........
c     ================================================================
c
c        -------------------------------------------------------------
      vxout=0.d0
      if(dz.lt.-1.d0)write(*,*)'newexch, dz,ro3',dz,ro3
      if(ro3.lt.1.d9)vxout=alpha2(ro3,dz,sp,iexch,excout)
c        -------------------------------------------------------------
      do ir=1,jmt
         rhot=rhoup(ir)+(n_spin_pola-1)*rhodn(ir)
c        if(rhot.gt.zero) then
	 if(rhot.gt.1.d-9) then
           dzr=(rhoup(ir)-rhodn(ir))/rhot
           if(dzr.gt.+1.d0)dzr=+1.d0
           if(dzr.lt.-1.d0)dzr=-1.d0
           ro3r=(three*r_mesh(ir)**2/rhot)**third
           vx(ir)=alpha2(ro3r,dzr,sp,iexch,enxc(ir))
           if(.not.(enxc(ir).lt.0))then
             write(6,*)'newexch',ir,rhot,dzr,ro3r,enxc(ir)
             enxc(ir)=0.d0
           endif
	 else
	   vx(ir)=zero
	   enxc(ir)=zero
	 endif
      enddo
c
      return
      end

c     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine get_rho(rho_in,
     >                  rho_out, stride,
     >                  r_mesh,jmt)
c     ================================================================
c
      implicit   none
c
      integer    n_spin_pola
      integer    jmt
      integer   stride
c
      integer    ir
c
      real*8     rho_in(jmt)
      real*8     rho_out(stride*jmt)
      real*8     r_mesh(jmt)
c
      real*8     ro3r
      real*8     pi
      real*8     four
      real*8     three
      real*8     third
c
      parameter  (three=3.0d0)
      parameter  (third=1.0d0/three)
      parameter  (four = 4.0d0)
c
c     -------------------------------------------------------------
      pi = 4.d0*datan(1.d0)

      do ir=1,jmt
           ro3r=(three*r_mesh(ir)**2/rho_in(ir))**third
           rho_out(stride*(ir-1)+1) = three / (four * pi * ro3r**3)
      enddo
c
      return
      end


text:
      subroutine phase7m(r,p1,p2,p3,p4,p5,p6,p7,p8,p9,wt,*)
      implicit none
      include 'constants.f'
      include 'masses.f'
      include 'mxdim.f'
      include 'breit.f'
c******* generate phase space for 2-->4 process
c******* r(mxdim),p1(4),p2(4) are inputs reversed in sign from physical values
c---- phase space for -p1-p2 --> p3+p4+p5+p6+p7+p8+p9+p10
c---- with all 2 pi's (ie 1/(2*pi)^20)
      integer j
      double precision r(mxdim)
      double precision p1(4),p2(4),p3(4),p4(4),p5(4),p6(4),p7(4),p8(4),
     . p9(4),p12(4),
     . p345(4),p678(4),p45(4),p78(4),p345678(4),
     . smin,wt,wt0,wt129,wt345678,wt345,wt678,wt45,wt78
      parameter(wt0=1d0/twopi**5)
      wt=0d0
      do j=1,4
      p12(j)=-p1(j)-p2(j)
      enddo
      smin=100d0

      n2=1
      n3=1

      mass3=hmass
      width3=hwidth
      call phi1_2m(0d0,r(1),r(2),r(3),smin,p12,p9,p345678,wt129,*99)

      n2=0
      n3=0

      mass2=mtau
      width2=tauwidth
      mass3=mtau
      width3=tauwidth

      call phi1_2(r(4),r(5),r(6),r(7),p345678,p345,p678,wt345678,*99)

      n2=1
      n3=1

      mass3=wmass
      width3=wwidth
      call phi1_2m(0d0,r(8),r(9),r(10),smin,p345,p3,p45,wt345,*99)
      call phi1_2m(0d0,r(11),r(12),r(13),smin,p678,p6,p78,wt678,*99)

      if ((p3(4).le.0d0).or.(p6(4).le.0d0)) goto 99
      call phi3m0(r(14),r(15),p45,p4,p5,wt45,*99)
      if ((p4(4).le.0d0).or.(p5(4).le.0d0)) goto 99
      call phi3m0(r(16),r(17),p78,p7,p8,wt78,*99)
      if ((p7(4).le.0d0).or.(p8(4).le.0d0)) goto 99

      wt=wt0*wt129*wt345678*wt345*wt678*wt45*wt78

      return
 99   wt=0d0
c      write(*,*) 'wt129',wt129
c      write(*,*) 'wt345678',wt345678
c      write(*,*) 'wt345',wt345
c      write(*,*) 'wt678',wt678
c      write(*,*) 'wt45',wt45
c      write(*,*) 'wt78',wt78
      return 1
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen