text:
      subroutine qq_tchan_htqg(p,msq)
c---Matrix element squared averaged over initial colors and spins
c     u(-p1)+b(p2)->h(p3,p4)+t(p5)+d(p6)+g(p7)
      implicit none
      include 'ewcouple.f'
      include 'constants.f'
      include 'qcdcouple.f'
      include 'hdecaymode.f'
      include 'masses.f'
      include 'nwz.f'
      include 'zprods_com.f'
      integer j,nu,eta,e5
      double precision p(mxpart,4),q(mxpart,4),p5Deta,s34,hdecay,
     & msq(-nf:nf,-nf:nf),dot,msqgamgam,fac,ampsqL,ampsqH,
     & b_u,u_b,db_b,b_db,d_bb,ub_bb,bb_d,bb_ub,
     & db_g,g_b,u_g,g_db,b_g,g_u,
     & d_g,g_bb,ub_g,g_d,bb_g,g_ub

C   Deal with Higgs decay
      s34=(p(3,4)+p(4,4))**2-(p(3,1)+p(4,1))**2
     &   -(p(3,2)+p(4,2))**2-(p(3,3)+p(4,3))**2
      if (hdecaymode == 'tlta') then
          call htautaudecay(p,3,4,hdecay)
      elseif (hdecaymode == 'bqba') then
          call hbbdecay(p,3,4,hdecay)
      elseif (hdecaymode == 'gaga') then
          hdecay=msqgamgam(hmass)
      else
      write(6,*) 'Unimplemented process'
      stop
      endif
      hdecay=hdecay/((s34-hmass**2)**2+(hmass*hwidth)**2)

C   Construct demassified momentum for p5 wrt eta and store in position e5
C   Set eta=3 for now
      eta=3
      p5Deta=dot(p,5,eta)

      e5=8
      do nu=1,4
      do j=1,7
      q(j,nu)=p(j,nu)
      enddo
      q(e5,nu)=mt**2*p(eta,nu)/(2d0*p5Deta)
      q(5,nu)=p(5,nu)-q(e5,nu)
      enddo

C---zero all matrix elements
      msq(:,:)=0d0
      fac=xn**2*(2d0*CF)*gwsq**3*gsq*hdecay

      if (nwz .eq. 1) then
      call spinoru(8,q,za,zb)

      call ubhtdgsq(1,2,3,4,5,6,7,e5,ampsqL,ampsqH)
      u_b=aveqq*fac*(ampsqL+ampsqH)
      call ubhtdgsq(2,1,3,4,5,6,7,e5,ampsqL,ampsqH)
      b_u=aveqq*fac*(ampsqL+ampsqH)
      call ubhtdgsq(6,2,3,4,5,1,7,e5,ampsqL,ampsqH)
      db_b=aveqq*fac*(ampsqL+ampsqH)
      call ubhtdgsq(6,1,3,4,5,2,7,e5,ampsqL,ampsqH)
      b_db=aveqq*fac*(ampsqL+ampsqH)
c--- In g-b diagrams, remove corrections on heavy line
      call ubhtdgsq(7,2,3,4,5,6,1,e5,ampsqL,ampsqH)
      g_b=2d0*aveqg*fac*(ampsqL+zip*ampsqH)
      call ubhtdgsq(7,1,3,4,5,6,2,e5,ampsqL,ampsqH)
      b_g=2d0*aveqg*fac*(ampsqL+zip*ampsqH)
c--- In g-q diagrams, remove corrections on light line
      call ubhtdgsq(2,7,3,4,5,6,1,e5,ampsqL,ampsqH)
      g_u=aveqg*fac*(zip*ampsqL+ampsqH)
      call ubhtdgsq(1,7,3,4,5,6,2,e5,ampsqL,ampsqH)
      u_g=aveqg*fac*(zip*ampsqL+ampsqH)
      call ubhtdgsq(6,7,3,4,5,1,2,e5,ampsqL,ampsqH)
      db_g=aveqg*fac*(zip*ampsqL+ampsqH)
      call ubhtdgsq(6,7,3,4,5,2,1,e5,ampsqL,ampsqH)
      g_db=aveqg*fac*(zip*ampsqL+ampsqH)

      msq(-1,+5)=db_b
      msq(-3,+5)=db_b
      msq(+2,+5)=u_b
      msq(+4,+5)=u_b

      msq(+5,-1)=b_db
      msq(+5,-3)=b_db
      msq(+5,+2)=b_u
      msq(+5,+4)=b_u

      msq(+2,0)=u_g
      msq(+4,0)=u_g
      msq(-1,0)=db_g
      msq(-3,0)=db_g
      msq( 0,5)=g_b

      msq(0,+2)=g_u
      msq(0,+4)=g_u
      msq(0,-1)=g_db
      msq(0,-3)=g_db
      msq(5,0)=b_g

      elseif(nwz .eq. -1) then
      call spinoru(8,q,zb,za)

      call ubhtdgsq(1,2,4,3,5,6,7,e5,ampsqL,ampsqH)
      d_bb=aveqq*fac*(ampsqL+ampsqH)
      call ubhtdgsq(2,1,4,3,5,6,7,e5,ampsqL,ampsqH)
      bb_d=aveqq*fac*(ampsqL+ampsqH)
      call ubhtdgsq(6,2,4,3,5,1,7,e5,ampsqL,ampsqH)
      ub_bb=aveqq*fac*(ampsqL+ampsqH)
      call ubhtdgsq(6,1,4,3,5,2,7,e5,ampsqL,ampsqH)
      bb_ub=aveqq*fac*(ampsqL+ampsqH)

c--- In g-b diagrams, remove corrections on heavy line
      call ubhtdgsq(7,2,4,3,5,6,1,e5,ampsqL,ampsqH)
      g_bb=2d0*aveqg*fac*(ampsqL+zip*ampsqH)
      call ubhtdgsq(7,1,4,3,5,6,2,e5,ampsqL,ampsqH)
      bb_g=2d0*aveqg*fac*(ampsqL+zip*ampsqH)
c--- In g-q diagrams, remove corrections on light line
      call ubhtdgsq(2,7,4,3,5,6,1,e5,ampsqL,ampsqH)
      g_d=aveqg*fac*(zip*ampsqL+ampsqH)
      call ubhtdgsq(1,7,4,3,5,6,2,e5,ampsqL,ampsqH)
      d_g=aveqg*fac*(zip*ampsqL+ampsqH)
      call ubhtdgsq(6,7,4,3,5,1,2,e5,ampsqL,ampsqH)
      ub_g=aveqg*fac*(zip*ampsqL+ampsqH)
      call ubhtdgsq(6,7,4,3,5,2,1,e5,ampsqL,ampsqH)
      g_ub=aveqg*fac*(zip*ampsqL+ampsqH)

      msq(+1,-5)=d_bb
      msq(+3,-5)=d_bb
      msq(-2,-5)=ub_bb
      msq(-4,-5)=ub_bb

      msq(-5,+1)=bb_d
      msq(-5,+3)=bb_d
      msq(-5,-2)=bb_ub
      msq(-5,-4)=bb_ub

      msq(+1,0)=d_g
      msq(+3,0)=d_g
      msq(-2,0)=ub_g
      msq(-4,0)=ub_g
      msq(0,-5)=g_bb

      msq(0,+1)=g_d
      msq(0,+3)=g_d
      msq(0,-2)=g_ub
      msq(0,-4)=g_ub
      msq(-5,0)=bb_g

      endif
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
! Test program for the module "XMLPARSE"
! Read a set of sample XML-files and report the
! possible failure to the screen
!
! $Id: tstparse.f90,v 1.3 2006/03/26 19:05:48 arjenmarkus Exp $
!
! Per test case read:
! - The title
! - The XML file to be
! - The expected parse data:
!   - which tags and the attributes that belong to them
!   - which data lines
!
! If the XML file was read but the actually reported data did
! not match the expected ones, then the title of the test
! case is printed on screen.
!
! So, a successful test does not produce output.
!
! The test cases are defined in a file "tstpars.inp"
! which has the following format:
! - Title line
! - A line which starts with "-----"
! - The contents of the XML file
! - Another line with "-----"
! - Lines with keywords (in order) that define what
!   tags, attributes and data to expect:
!   tag: string
!   endtag: T/F
!   attribute: A=B
!   data: string
! - Another line with "-----" to separate it from the
!   next case
! - Lines with "#" in the first column are considered
!   comment and are therefore ignored
!
program tstparse

   use xmlparse

   implicit none

   character(len=60) :: title
   character(len=80) :: string
   logical           :: has_title
   integer           :: ierr

   !
   ! Start the loop
   !
   open( 10, file = 'tstparse.inp', status = 'old' )
   open( 20, file = 'tstparse.out' )

   has_title = .false.

   do while ( ierr .eq. 0 )
      read( 10, '(a)', iostat = ierr ) string
      !
      ! End of file or some error condition
      !
      if ( ierr .ne. 0 ) then
         exit
      endif
      !
      ! Skip comments
      !
      if ( string(1:1) .eq. '#' ) then
         cycle
      endif
      !
      ! Store the title and start reading the XML data
      !
      if ( .not. has_title ) then
         title = string
         has_title = .true.
      endif
      if ( has_title .and. string(1:5) .eq. '-----' ) then
         open( 30, file = 'tstparse.xml' )
         do
            read(  10, '(a)', iostat = ierr ) string
            if ( ierr .ne. 0 ) then
               exit
            endif
            if ( string(1:5) .eq. '-----' ) then
               exit
            endif
            write( 30, '(a)' ) string
         enddo

         !
         ! Now we are ready to read the XML file
         ! and analyse the report
         !
         close( 30 )
         call check_xml_file
      endif
   enddo

contains

subroutine check_xml_file
   logical           :: mustread
   type(XML_PARSE)   :: info

   character(len=80)                      :: tag
   logical                                :: endtag
   character(len=80), dimension(1:2,1:20) :: attribs
   integer                                :: no_attribs
   character(len=200), dimension(1:100)   :: data
   integer                                :: no_data
   integer                                :: i

   call xml_open( info, 'tstparse.xml', .true. )
   call xml_options( info, report_lun = 20, report_details = .true. )

   do
      call xml_get( info, tag, endtag, attribs, no_attribs, data, no_data )
      if ( xml_error(info) ) then
         write(*,*)  'Error ', title
         write(20,*) 'Error ', title
      endif

      write( 20,* ) tag, endtag
      do i = 1,no_attribs
         write( 20,'(i3,1x,3a)') &
            i, trim(attribs(1,i)), '=', trim(attribs(2,i))
      enddo
      write( 20,'(i3,1x,3a)') (i, '>',trim(data(i)), '<', i=1,no_data)
      if ( .not. xml_ok(info) ) exit
   enddo

   do
      read(10,'(a)' ) string
      if ( string(1:5) .eq. '-----' ) then
         exit
      endif
   enddo

   call xml_close( info )

end subroutine check_xml_file
end program tstparse

max_stars_repo_name: paulromano/xml-fortran

max_issues_repo_name: paulromano/xml-fortran

max_forks_repo_name: paulromano/xml-fortran


text:
Superfund is the common name for CERCLA, or Comprehensive Environmental Response, Compensation, and Liability Act,  which was enacted by the United States Congress on December 11, 1980 in response to the Love Canal disaster. This law created a tax on petroleum and chemical industries and provided broad Federal authority to respond directly to releases or threatened releases of hazardous substances that may endanger public health or the environment. Over five years, $1.6 billion was collected and the tax went to a trust fund for cleaning up abandoned or uncontrolled hazardous waste sites. (wiki:WikiPedia:CERCLA Wikipedia)

According to the EPA Website there are three Superfund sites in Davis:

Davis Transmitter Site

The Davis Transmitter Site was a outofblastradius transmission site for McClellan Air Force Base.  There was a radio tower that was the backup in case McClellan got nuked.  Just about every AFB has one.  It was decommissioned (presumably when McClellan was), and the only reason its on the Superfund is because they had a septic tank there and anything military automatically goes onto Superfund if there are any groundwater issues.  There are also minor runoffs consistant with a military installation (namely, lead).  In this case, the pollutant is primarily... airman poop.  It was discovered in 1988 and had its first and last inspection in 1994. This site is located between County Road 35 and County Road 36 East of Yolo County Grasslands Park.  It is visible by all the odd antennas sticking out of the ground. (http://cfpub1.epa.gov/supercpad/cursites/cactinfo.cfm?id0903522 EPA Database)

Frontier Fertilizer

Located on 2nd Street (aka County Road 32a) between Cantrill Drive and Mace Boulevard Mace, it is a large fencedoff area which is designated by a extremely faded sign that warns trespassers of danger and to keep out.  A http://www.atsdr.cdc.gov/HAC/PHA/frontier/frf_toc.html public health assessment was last filed in 1995 for this site. (http://cfpub1.epa.gov/supercpad/cursites/csitinfo.cfm?id0901554 EPA Database).  The EPA website contains a http://www.epa.gov/superfund/sites/npl/nar1413.htm narrative describing this site and possibly other information.  In a general way, the site has both soil and groundwater contamination due to the production of fertilizers for agriculture. The fence is to keep people away from the surface soil contamination, and there is a groundwater treatment system that is trying to clean up the groundwater.  As no one is currently being poisoned from the site, the EPA is not heavily invested in doing treatment any faster (they have other places to spend resources where people are getting sick, or at  least the potential is much higher).

The new Target is going to be placed immediately east of this site. One major source of concern surrounding this placement is that it will remove some of the wells that permit the monitoring of the rate of spread of the contamination of the groundwater. Target has said that they will replace the wells, but others say that the new wells may not provide the same quality of information as the originals. The construction of new monitoring wells is being overseen by the EPA.  Also see the Frontier Fertilizer Superfund Oversight Group.

Pilau Drum Site

This is another Superfund site, on the UC Davis Campus, but because it is south of I80 it is technically in Solano County.   It was UC Davis landfill LEHR/Old Campus Landfill which is now called the Center for Health and the Environment on Old Davis Road. This one is particularly interesting because it is the site of the infamous Radioactive Beagles that is discussed on the Urban Legends page.  It is located south of I80 at the UC Davis exit, in the area known as South Campus, near the California Raptor Center.  (http://cfpub1.epa.gov/supercpad/cursites/csitinfo.cfm?id0905944 EPA Database)

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE SOLVE1(A1,R1,RP,XI,LAM2,LAM3,LAM4,CONT)
C
C     ROUTINE TO SOLVE FOR LAMBDAS AS FNCTS. OF XI
C
C
      REAL LAM2,LAM3,LAM4
C
      IF (RP.EQ.0.0) GO TO 20
C
      SUM = A1 + XI / RP
      SINSUM = SIN(SUM)
      BB = R1 - RP * (SIN(A1) - SINSUM)
      RT = 0.0E0
      IF( SINSUM .NE. 0.0E0 ) RT = BB / SINSUM
      PSI1 = COS(SUM)
      PSI2 = -SINSUM / RP
C
C     CHECK FOR SHELL CAP CASE
      IF ( A1 .NE. 0.0 )  GO TO 40
      LAM2  = 0.0E0
      IF( BB .NE. 0.0E0 ) LAM2  = PSI1 / BB
      LAM3  =  1.0 / RP
      LAM4  = -1.0 / RP**2
      GO TO 50
C
C     ALF1 = ALF2
C
   20 SINA = SIN(A1)
      COSA = COS(A1)
      BB = R1 + XI * COSA
      RT = 0.0E0
      IF( SINA .NE. 0.0E0 ) RT = BB / SINA
      PSI1 = COSA
      PSI2 = 0.0
C
   40 LAM2 = 0.0E0
      IF( BB .NE. 0.0E0 ) LAM2 = PSI1 / BB
      LAM3 = 0.0E0
      IF( RT .NE. 0.0E0 ) LAM3 = 1.0E0 / RT
      LAM4 = 0.0E0
      IF( BB .NE. 0.0E0 ) LAM4 = PSI2 / BB
C
   50 CONTINUE
      RETURN
      END


text:
The http://www.facebook.com/pages/UCDLoveLab/250694719161 Love Lab is a mobile cart stocked with free safer sex products for students of UC Davis. The Love Lab travels around campus with the Health Education and Promotion volunteers. Every time you visit the Love Lab you can get two packets of waterbased lubricant, ten free condoms, flavored dental dams, female insertive condoms, nonlatex condoms and pamphlets for free.

For a list of condom prices at stores throughout Davis, check out Dr. Loves Condom Buying Guide for Davis.

 Where Can You Find the Love Lab?
It is located in the Student Health and Wellness Center located on La Rue, across from the Activities and Recreation Center (ARC). The Love Lab is located inside Health Education and Promotion on the third floor. It is available Monday, Tuesday, Thursday, and Friday from 8:00am to 4:30pm and Wednesday from 9:00am to 4:30pm.  The Love Lab also visits the ARC Lobby every Thursday from 6:308:30pm, unless it is raining.  Come visit to check out all our products and get a product demonstration from one of our volunteers!

 Love Lab Events:

Check back in September for Fall Quarter Events.

Products
The Love Lab carries all of the following products.

Crown Skin Less Skin  Thin and ultra sensitive. The closest thing to nothing at all. Lubricated, Latex, Noncontoured.

Fantasy Assorted Flavors  Lightly flavored with six stimulating flavors: vanilla, strawberry, mint, grape, chocolate and banana. Lubricated, Latex, Flavored.

Impulse Bare Pleasure  A ribbed and studded condom for heightened sensation. Slight rose tint. Lubricated, Latex, Ribbed and Studded.

Kimono Thin  Extra sensitive, strong, reliable and 20% thinner than competitors standard condoms. Distinguished, Japanese made, the original Kimono lubricated condom. Lubricated, Latex, Formfitting, reservoir tip.

LifeStyles Snugger Fit   Shaped for more contact and feeling. Designed to hug better so it feels better. Lubricated, Latex, Fitted Shape.

Trustex Aggie Spirit  Show your Aggie Spirit with these blue and gold condoms. Lubricated, Latex.

Trustex Assorted Colors  Bright and vibrant condoms that come in a variety of colors. Lubricated, Latex.

Trustex Extra Large  One of the most sensitive extra large condoms manufactured. Ideal for those who prefer a roomier condom. Lubricated, Latex, Larger Size.

Trustex Mint  Extra minty flavor for a fresh experience for both partners. Nonlubricated, Latex, Flared tip. Intended for oral sex.

Trustex Ribbed & Studded  Ribbed and studded texture provides the ultimate experience for both partners. Lubricated, Latex, Ribbed and Studded.

LifeStyles Skyn  A a nonlatex, polyisoprene condom. The newest latex alternative on the market. Available Upon Request. Lubricated, polyisoprine

FC2 Insertive Condom  Designed for men and women who want a pleasurable alternative contraceptive to the male condom. Available Upon Request. Lubricated, Polyurethane

Flavored Latex Dental Dams  Lightly flavored and scented. Silky thin 6X8 latex sheet specially designed for use as a barrier while performing oral sex. Available Upon Request. Nonlubricated, Latex

ID Glide Lubricant  Ultra long lasting. Nonstaining, odorless, clear formula. Water based, will not harm latex. Extensively tested for safety. Silky smooth. Doctor recommended brand. Does not contain nonoxynol9. Two Per Person

Watch It, Teach It, Win It!
Stop by the Love Lab to Watch it, Teach it, and Win it! Watch any one of the three demonstrations (male condom, insertive condom, or dental dam), teach it back to a Love Lab volunteer, and you’ll be one step closer to winning a Tshirt or other prize (while supplies last). Complete all three and visit Health Education and Promotion on the 3rd floor of the Student Health and Wellness Center to pick up your prize. On your next visit to the Love Lab talk to a volunteer for more details!


text:
      SUBROUTINE DPMINI(Lesdpm,Loudpm,Ifdpm,Pinp,Idpmvr,Idpmhk,Idpmfs)

C***********************************************************************
C                                                                      *
C     Version September 2001      by   Stefan Roesler                  *
C                                                                      *
C     Last change  10-Nov-15           Alfredo Ferrari                 *
C                                                                      *
C     This subroutine is part of the FLUKA interface to DPMJET 3,      *
C     new version (labelled 3.10 for convenience) after A.Fedynitch    *
C     work                                                             *
C     Initialization of DPMJET 3 event generation.                     *
C                                                                      *
C                                                                      *
C     argument list (FLUKA control card: codewd = dpmjet)              *
C                                                                      *
C          LESDPM  dpmjet input follows at the bottom of Fluka input   *
C                  file and a second start card is required in this    *
C                  case                                                *
C          LOUDPM  logical unit number for dpmjet output               *
C          IFDPM   this flag indicates if (and at what level?) we      *
C                  expect dpmjet output                                *
C                  for dpmjet-3 just on/off
C                                                                      *
C          PINP    initial FLUKA beam momentum to indicate maximum     *
C                  available energy in the considered system           *
C          IDPMVR  output flag indicating service is rendered by       *
C                  dpmjet-3                                            *
C          IDPMHK  carry over the dpmet(2/3) event common block size   *
C                  to enable important cross check since we have in    *
C                  principle 3 actual copies of it                     *
C          IDPMFS  string fusion flag                                  *
C                                                                      *
C***********************************************************************
C

C      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'

      INCLUDE '(BEAMCM)'
      INCLUDE '(IOIOCM)'

      DOUBLE PRECISION epn , Pinp , xdumb , xlim1 , xlim2 , xlim3
      INTEGER ibin , idp , IDPmev , Idpmfs , Idpmhk , Idpmvr , Ifdpm ,
     &        iglau , iglaub , IHEhad , IHEnuc , IHIjpr , IHMapr ,
     &        IHMata , Lesdpm , Loudpm , ninp , npchar , npmass
      INTEGER ntchar , ntmass

C emulsion treatment
      INCLUDE 'inc/dtcomp'

      INCLUDE 'inc/dtevt1'
C event flag
      INCLUDE 'inc/dtevno'
C Glauber formalism: flags and parameters for statistics
      INCLUDE 'inc/dtglgp'
C flags for input different options
      INCLUDE 'inc/dtflg1'
C
      INCLUDE 'inc/pydat1'

C histogram indices for Fluka-interface related statistics
C     CHARACTER*72 HEADER
      DIMENSION xdumb(40)
      COMMON /DTFLHX/ IHMapr , IHMata , IHIjpr , IHEnuc , IHEhad ,
     &                IDPmev

      INCLUDE 'inc/dtflka'

C
C Redirect Dpmjet and Pythia output and turn on/off output
C
      Idpmhk = NMXHKK
      Idpmvr = 310
      LPRi = Ifdpm

      LOUt = Loudpm
      MSTu(11) = Loudpm
C
C Flag for special settings needed to run the code as event
C  generator in Fluka (do not change !)
C
      ITRspt = 1
      IEMul = 0
      IFUsion = Idpmfs
      epn = Pinp                !---> should we  * 1.5
      npmass = IPRoa            !
      npchar = IPRoz            !---> last two are taken out of FLUKA common
C
C Special settings if no Dpmjet input follows the Fluka input
C
      IF ( Lesdpm.EQ.0 ) THEN
         ninp = -1
      ELSE
         ninp = 0
      END IF

C
C Initialization of Dpmjet
C
Cc
Cc ---------------------------------------------------------------------

      iglaub = 0

      CALL DT_INIT(ninp,epn,npmass,npchar,ntmass,ntchar,idp,iglau)
      CALL DT_STATIS(1)

Cc ---------------------------------------------------------------------
Cc
C     CALL DT_DTUINI(NINP,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,IEMU)
      NEVent = 0
C

      IF ( LPRi.LT.3 ) RETURN
C
C Initialization of histograms
C
      IDPmev = 0
C  mass number of projectile and target nuclei
      xlim1 = 0.5D0
      xlim2 = 200.5D0
      xlim3 = 0D0
      ibin = INT(xlim2-xlim1)
      CALL DT_NEWHGR(xlim1,xlim2,xlim3,xdumb,ibin,IHMapr)
      CALL DT_NEWHGR(xlim1,xlim2,xlim3,xdumb,ibin,IHMata)
C  index of projectile hadrons
      CALL DT_NEWHGR(xlim1,xlim2,xlim3,xdumb,ibin,IHIjpr)
C  energy of projectile nuclei and hadrons
      xlim1 = 5.0D0
      xlim2 = PBEam
C     write(0,*) ' -pbeam-',pbeam
      xlim3 = 0.D0
      ibin = -200
      CALL DT_NEWHGR(xlim1,xlim2,xlim3,xdumb,ibin,IHEnuc)
      CALL DT_NEWHGR(xlim1,xlim2,xlim3,xdumb,ibin,IHEhad)
C

C=== End of subroutine Dpmini =========================================*
      END SUBROUTINE

max_stars_repo_name: afedynitch/DPMJET

max_issues_repo_name: DPMJET/DPMJET

max_forks_repo_name: afedynitch/DPMJET


text:
module rcrlib_gnu
    use, intrinsic :: ieee_arithmetic ! requires gfortran version 5.0 or higher
    implicit none
    integer, parameter :: SP=kind(1.0), DP=selected_real_kind(9,99)
    integer, parameter :: stderr=0
    public SP, DP, is_finite, stderr, initialize

contains

subroutine initialize(infty,nan)
    real(kind=DP), intent(out) :: infty, nan
    infty = 1.0_dp/0.0_dp ! huge(1.0_dp)
    nan = 0.0_dp/0.0_dp
end subroutine initialize

elemental function is_finite(x)
    real(kind=DP), intent(in) :: x
    logical :: is_finite
    is_finite = ieee_is_finite(x) ! This call requires "ieee_arithmetic"
end function is_finite

end module rcrlib_gnu

max_stars_repo_name: bvkrauth/rcr

max_issues_repo_name: bvkrauth/rcr

max_forks_repo_name: bvkrauth/rcr


text:
c     Loads all files into memory so they don't have to be reread on
c     every routine call.
c
c     Changelog:
c       2015-04-16     Created (Joseph Nicholas)
c       2015-04-27     Bugfix in load_ap (Joseph Nicholas)

      block data common_data
          include 'common.inc'
          data FILMOD /'dgrf45.dat', 'dgrf50.dat', 'dgrf55.dat',
     &                 'dgrf60.dat', 'dgrf65.dat', 'dgrf70.dat',
     &                 'dgrf75.dat', 'dgrf80.dat', 'dgrf85.dat',
     &                 'dgrf90.dat', 'dgrf95.dat', 'dgrf00.dat',
     &                 'igrf05.dat', 'igrf05s.dat'/
          data DTEMOD / 1945., 1950., 1955., 1960., 1965., 1970.,
     &                  1975., 1980., 1985., 1990., 1995., 2000.,
     &                  2005., 2010./
      end

      subroutine load_shc(ifile, filename, ierror, errmsg)
          character filename*(*)
          include 'common.inc'

          integer   ierror
          character*256 errmsg

          parameter (iunit = 84)
          character*256 fullpath

c         Special for IRIWeb version
c1001     format('/usr/local/etc/httpd/cgi-bin/models/IRI/',A12)
c1001     format(A12)
c         For GMAT version
 1001     format('./../data/IonosphereData/',A12)

          write(fullpath, 1001) filename

          ierror = 0
          if (ifile > MAXFILES_SHC) then
c              write(*,*) 'ERROR: load_shc: ifile exceeds MAXFILES_SHC'
c              stop 1
             ierror = 1
             errmsg = 'ERROR: load_shc: ifile exceeds MAXFILES_SHC'
             return
          end if

c         Open coefficient file. Read past first header record.
          open(iunit, file=fullpath, status='old', err = 2001)

c         Read degree and order of model and Earth's radius.
          read(iunit, *)
          read(iunit, *) nmax_shc(ifile), erad_shc(ifile), yy
c          print * ,'Read ', fullpath
c          print * , nmax_shc(ifile), '  ', erad_shc(ifile), '  ', yy

c         --------------------------------------------------------------
c         Read the coefficient file, arranged as follows:
c
c                                         N     M     G     H
c                                         ----------------------
c                                     /   1     0    GH(1)  -
c                                    /    1     1    GH(2) GH(3)
c                                   /     2     0    GH(4)  -
c                                  /      2     1    GH(5) GH(6)
c             NMAX*(NMAX+3)/2     /       2     2    GH(7) GH(8)
c                records          \       3     0    GH(9)  -
c                                  \      .     .     .     .
c                                   \     .     .     .     .
c             NMAX*(NMAX+2)          \    .     .     .     .
c             elements in GH          \  NMAX  NMAX   .     .
c
c         N and M are, respectively, the degree and order of the
c         coefficient.
c         --------------------------------------------------------------
          irec = 0
          do nn = 1, nmax_shc(ifile)
              do mm = 0, nn
                  irec = irec + 1

                  if (irec > MAXRECORDS_SHC) then
c                      write(*,*)
c     &                    'ERROR: load_shc: irec exceeds MAXRECORDS_SHC'
c                      stop 1
                     ierror = 2
                errmsg = 'ERROR: load_shc: irec exceeds MAXRECORDS_SHC'
                return
                  end if

                  read(iunit, *) (shcfiles(i, irec, ifile), i=1,4)
c                  print *, shcfiles(1, irec, ifile), '  ',
c     & shcfiles(2, irec, ifile), '  ',
c     & shcfiles(3, irec, ifile), '  ',
c     & shcfiles(4, irec, ifile)
              end do
          end do

          close(iunit)
          errmsg = ''
          return

 2001     ierror = 1000 + ierror
          errmsg = 'Error: can not open "'//fullpath
          errmsg = errmsg//'"'
          return
      end

      subroutine load_igrz(filename, ierror, errmsg)
          character filename*(*)
          include 'common.inc'

          integer   ierror
          character*256 errmsg

          parameter (iunit = 84)
          character*256 fullpath

c         Special for IRIWeb version
c1001     format('/usr/local/etc/httpd/cgi-bin/models/IRI/',A)
c1001     format(A)
c         For GMAT version
 1001     format('./../data/IonosphereData/',A)

          write(fullpath, 1001) filename

          ierror = 0
          errmsg = ''
          open(iunit, file=fullpath, status='old', err = 3001)

c         Read the update date, the start date and the end date
c         (mm,yyyy), and get number of data points to read.
          read(iunit, *) iupd_igrz, iupm_igrz, iupy_igrz
          read(iunit, *) imst_igrz, iyst_igrz, imend_igrz, iyend_igrz
          inum_vals = 3 - imst_igrz + (iyend_igrz-iyst_igrz)*12
     &                + imend_igrz

          if (inum_vals > MAXRECORDLEN_IGRZ) then
c              write(*,*)
c     &           'ERROR: load_igrz: inum_vals exceeds MAXRECORDLEN_IGRZ'
c              stop 1
        ierror = 3
        errmsg = 'ERROR: load_igrz: inum_vals exceeds MAXRECORDLEN_IGRZ'
        return
          end if

c         Read data records
          read(iunit, *) (ionoindx_igrz(i), i=1,inum_vals)
          read(iunit, *) (indrz_igrz(i), i=1,inum_vals)

          close(iunit)
          return

 3001     ierror = 1000 + ierror
          errmsg = 'Error: can not open "'//fullpath
          errmsg = errmsg//'"'
          return
      end

      subroutine load_ap(filename, ierror, errmsg)
          character filename*(*)
          include 'common.inc'

          integer   ierror
          character*256 errmsg

          parameter (iunit = 84)
          character*256 fullpath

c         Special for IRIWeb version
c1001     format('/usr/local/etc/httpd/cgi-bin/models/IRI/',A)
c1001     format(A)
c         For GMAT version
 1001     format('./../data/IonosphereData/',A)

          write(fullpath, 1001) filename

          ierror = 0
          errmsg = ''
          open(iunit, file=fullpath, status='old', err = 4001)

          do irec = 1, MAXRECORDS_AP+1
            if (irec > MAXRECORDS_AP) then
c                write(*,*) 'ERROR: load_ap: irec exceeds MAXRECORDS_AP'
c                stop 1
               ierror = 4
               errmsg = 'ERROR: load_ap: irec exceeds MAXRECORDS_AP'
            end if

            read(iunit, 1002, end=2001) (ints_ap(k,irec), k=1,11),
     &                                  reals_ap(irec)
 1002       format(3I3, 8I3, F5.1)
          end do

 2001     close(iunit)
          num_records_ap = irec - 1
          return

 4001     ierror = 1000 + ierror
          errmsg = 'Error: can not open "'//fullpath
          errmsg = errmsg//'"'
          return
      end

      subroutine load_ccir(imonth, ierror, errmsg)
          include 'common.inc'

          integer   ierror
          character*256 errmsg

          parameter (iunit = 84)
          character*256 fullpath

c         Special for IRIWeb version
c1001     format('/usr/local/etc/httpd/cgi-bin/models/IRI/ccir',
c    &           I2,'.asc')
c1001     format('ccir',I2,'.asc')
 1001     format('./../data/IonosphereData/ccir',I2,'.asc')

          write(fullpath, 1001) imonth+10

          ierror = 0
          errmsg = ''
          if (imonth > MAXFILES_CCIR) then
c              write(*,*)
c     &                  'ERROR: load_ccir: imonth exceeds MAXFILES_CCIR'
c              stop 1
             ierror = 5
             errmsg = 'ERROR: load_ccir: imonth exceeds MAXFILES_CCIR'
          end if

          open(iunit, file=fullpath, status='old', err = 5001)

          read(iunit, 1002)
     &               (((f2_ccir(i,j,k,imonth), i=1,13), j=1,76), k=1,2),
     &               (((f3_ccir(i,j,k,imonth), i=1, 9), j=1,49), k=1,2)
 1002     format(1X, 4E15.8)

          close(iunit)
          return

 5001     ierror = 1000 + ierror
          errmsg = 'Error: can not open "'//fullpath
          errmsg = errmsg//'"'
          return
      end

      subroutine load_ursi(imonth, ierror, errmsg)
          include 'common.inc'

          integer   ierror
          character*256 errmsg

          parameter (iunit = 84)
          character*256 fullpath

c         Special for IRIWeb version
c1001     format('/usr/local/etc/httpd/cgi-bin/models/IRI/ursi',
c    &           I2,'.asc')
c1001     format('ursi',I2,'.asc')
 1001     format('./../data/IonosphereData/ursi',I2,'.asc')

          write(fullpath, 1001) imonth+10

          ierror = 0
          errmsg = ''
          if (imonth > MAXFILES_URSI) then
c              write(*,*)
c     &                  'ERROR: load_ursi: imonth exceeds MAXFILES_URSI'
c              stop 1
             ierror = 6
             errmsg = 'ERROR: load_ursi: imonth exceeds MAXFILES_URSI'
          end if

          open(iunit, file=fullpath, status='old', err = 6001)

          read(iunit, 1002)
     &               (((f2_ursi(i,j,k,imonth), i=1,13), j=1,76), k=1,2)
 1002     format(1X, 4E15.8)

          close(iunit)
          return

 6001     ierror = 1000 + ierror
          errmsg = 'Error: can not open "'//fullpath
          errmsg = errmsg//'"'
          return
      end

      subroutine load_all_files(ierror, errmsg)
          include 'common.inc'

          integer   ierror
          character*256 errmsg

          do i = 1, MAXFILES_SHC
              call load_shc(i, FILMOD(i), ierror, errmsg)
              if (ierror .ne. 0) then
                  return
              end if
          end do

          call load_igrz('ig_rz.dat', ierror, errmsg)
          if (ierror .ne. 0) then
             return
          end if

          call load_ap('ap.dat', ierror, errmsg)
          if (ierror .ne. 0) then
             return
          end if

          do i = 1, MAXFILES_CCIR
              call load_ccir(i, ierror, errmsg)
              if (ierror .ne. 0) then
                 return
              end if
          end do

          do i = 1, MAXFILES_URSI
              call load_ursi(i, ierror, errmsg)
              if (ierror .ne. 0) then
                 return
              end if
          end do
      end

max_stars_repo_name: supalogix/GMAT

max_issues_repo_name: supalogix/GMAT

max_forks_repo_name: supalogix/GMAT


text:
      subroutine qlI3fin(Ival0,xpi,ier)
      implicit none
      double complex Ival0
      double precision xpi(6)
      integer ier
      ier = 0
      call ffxc0(Ival0, xpi, ier)
      end


text:
CHARACTER(3) :: a, b(3)
a = b(3)(:)
a = b(3)(1:)
a = b(3)(:2)
a = b(3)(1:2)
END PROGRAM

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


!=======================================================================
SUBROUTINE FAST (INPUTFILE,iSts,iStps,PtfmMotIn,PtfmVelIn,PtfmAcceIn)

   ! This program models 2- or 3-bladed turbines of a standard configuration.

USE                             ADAMSInput
USE                             DOFs
USE                             General
USE                             InitCond
USE                             Linear
USE                             SimCont
USE                             NWTC_Library
USE                             AeroDyn
USE                             FAST_IO_Subs       ! FAST_Begin(), FAST_Input(), PrintSum(), RunTimes()
USE                             FASTsubs           ! FAST_Initialize(), TimeMarch()
USE                             FAST2ADAMSSubs     ! MakeAdm(), MakeACF(), MakeACF_Lin
USE                             FAST_Lin_Subs      ! CalcSteady(), Linearize()
USE                             HydroDyn
USE                             Noise
USE                             Seismic
! Yang: USE AQWA & Output modules for the variables
USE                             AQWAinFAST
USE                             Output
USE                             TurbConf
! Yang: End of the change

IMPLICIT                        NONE

   ! Local variables:

!REAL(ReKi)                   :: GBoxTrq                                         ! Unused gearbox torque on the LSS side in N-m.

CHARACTER (LEN=1024)         :: INPUTFILE
!INTEGER(4)                   :: L                                              ! Generic index.
INTEGER(4)                   :: iSts                                           ! status of calling FAST. 0 = first,1 = the rest, 2 = the end.
INTEGER(4)                   :: iStps                                          ! number of steps for time marching solver.= dT_AQWA/dT_FAST. In the first round, shall plus 1.

INTEGER                       :: ErrStat

! Yang: define the platform motion and velocity matrix
REAL (Reki)           :: PtfmMotIn(6),PtfmVelIn(6),PtfmAcceIn(6)                    ! Platfrom motion, velocity and acceleration in the six DoFs
! End of the change

! Set Prifile
if (iSts == 0 .and. PrtFASTRslts) then      ! adding the condition (PrtFASTRslts) in this statement is just to avoid the second initilization when stage = 2 and time = 0
    Prifile = INPUTFILE
! Get the current time.

    CALL DATE_AND_TIME ( Values=StrtTime )                                           ! Let's time the whole simulation
    CALL CPU_TIME ( UsrTime1 )                                                       ! Initial time (this zeros the start time when used as a MATLAB function)

! Set version & initialize NWTC Library (open console, set pi constants)
    ! Yang: Not output FAST notice when couping with AQWA
    CALL SetVersion
    CALL NWTC_Init()                                                                 ! sets the pi constants

       ! Tell our nice users what they're running.
    ! Yang: Not output FAST notice when couping with AQWA
    CALL DispNVD()

       ! Open and read input files, initialize global parameters.

    CALL FAST_Begin()
    CALL FAST_Input()

   ! Set up initial values for all degrees of freedom.

    CALL FAST_Initialize()

   ! Print summary information to "*.fsm"?

    IF ( SumPrint )  CALL PrintSum
endif
! Yang: Set values for Q and QD with regards to the platform DoFs
IF (iSts == 0 .or. iSts == 1)  then
   PtfmMotAQWA   = PtfmMotIn
   PtfmVelAQWA   = PtfmVelIn
   PtfmAcceAQWA   = PtfmAcceIn
  ! write(*,*) 'Running TimeMarch'
   CALL TimeMarch (iSts,iStps)
! Rotor, Nacelle, Tower mass center
   ! CoG of rotor, tower and nacelle.

   CoGTwr(1) = CoGTwr(1) * AllOuts(YawBrTDxt) + CoGTwr(3)* sin(PtfmMotAQWA(5))
   CoGTwr(2) = CoGTwr(2) * AllOuts(YawBrTDyt) + CoGTwr(3)* sin(PtfmMotAQWA(4))
   CoGTwr(3) = CoGTwr(3)
   CoGRot(1) = CoGTwr(1) + HubCM + OverHang*cos(ShftTilt/180.0*3.1415926) + FASTHH * sin(PtfmMotAQWA(5))
   CoGRot(2) = CoGTwr(2) + FASTHH * sin(PtfmMotAQWA(4))
   CoGRot(3) = FASTHH
   CoGNac(1) = NacCMxn+CoGTwr(1) + (NacCMzn+TowerHt)*sin(PtfmMotAQWA(5))
   CoGNac(2) = NacCMyn+CoGTwr(2) + (NacCMzn+TowerHt)*sin(PtfmMotAQWA(4))
   CoGNac(3) = NacCMzn+TowerHt

ENDIF

if (iSts == 2) then
   CALL RunTimes
   CALL FAST_Terminate( ErrStat )
   CALL AD_Terminate(   ErrStat )
   CALL Hydro_Terminate( )
   CALL Noise_Terminate( )
   IF ( BEEP ) CALL UsrAlarm
   CALL NormStop( )
endif

END SUBROUTINE FAST
!=======================================================================


text:
I graduated from UC Davis in June 2007 as a Political Science and Religious Studies double major.  I moved to Portland, Oregon with Users/LindseyKrabbenhoft Lindsey after graduating, where we lived for four years before moving to Vancouver, B.C. in 2011.

While I was living in Davis I served as the Chair of the County Council of the Green Party of Yolo County, Chair of the ASUCD Elections Committee and President of the Davis College Green Party. Ive also been involved in the National Campus Greens (serving as a CoChair), the Progressive Coalition, a candidate for ASUCD Senate, and was involved in numerous other political causes. Obviously Im most proud of my fictional association with the fictional ASUCD Elections Committee Slate.

Before and during my time at UCD, the Davis College Green Party had notoriously large coffers and other likeminded individuals and organizations used to come to us for money. I always thought the picture on Users/LindseyKrabbenhoft this page should have been used by all of the organizations and people that were given monetary support from the Davis College Green Party. (She wasnt my girlfriend at the time we gave them the money, I saved cronyism like that for ASUCD.)

I have a Live Journal that I used to update but havent in quite awhile. Nevertheless, you can still find incriminating things I posted on it http://www.livejournal.com/users/jtleathers/ here.

I setup an archive of past ASUCD Elections with ballot data, reports, analysis, etc. The website can be found http://www.jtleathers.com/asucd/ here.

20060530 20:15:58 nbsp I voted today at the MU polling place and since the incumbant was running for County Clerk / Recorder I wrote Jonathon Leathers in as a writein candidate. Users/RobRoy
   This is made that much funnier by the fact that you compared me to this incumbent last year when you said I was the Freddie Oakley of ASUCD.

20100311 17:52:22 nbsp Hey, just wondering if there is a story behind http://daviswiki.org/ASUCD_Elections?actiondiff&version235&version134 this edit. Typo the first time around? Or does the first number mean something else and it just got dropped there by mistake? Users/WilliamLewis
   That original number of 6052 was derived from taking the total number of votes in the Senate election (36312) and dividing it by 6 (which is how many votes each voter got). It was a very rough estimate of the number of voters in the Senate election. If you go to the CEI page and check out the image on that page you can see the turnout was much higher. I have a little longer explanation http://www.jtleathers.com/asucd/2002_Fall/fall2002overview.html here as well.

20121204 08:33:29 nbsp Hopefully youll keep that site for awhile.  Thanks for updating the links!  (Oh, and you missed the one on this page) Users/JabberWokky
   Thanks, there were a few other links I missed as well that are now fixed. I plan to keep that data up indefinitely at the current URL.

20121204 11:45:15 nbsp I didnt realize you were still running data for the elections.  Love the elections website!  Thank you! Users/MaxwellKappes
   I intended to get everything updated and moved to a new URL and hosting server, it just took me a bit longer than planned. I plan to continue updating it for future elections too and its always good to know other people look at it.


text:
*DECK DU12US
      SUBROUTINE DU12US (A, MDA, M, N, B, MDB, NB, MODE, KRANK, RNORM,
     +   H, W, IR, IC)
C***BEGIN PROLOGUE  DU12US
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DULSIA
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (U12US-S, DU12US-D)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C        Given the Householder LQ factorization of A, this
C        subroutine solves the system AX=B. If the system
C        is of reduced rank, this routine returns a solution
C        according to the selected mode.
C
C       Note - If MODE.NE.2, W is never accessed.
C
C***SEE ALSO  DULSIA
C***ROUTINES CALLED  DAXPY, DDOT, DNRM2, DSWAP
C***REVISION HISTORY  (YYMMDD)
C   810801  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C***END PROLOGUE  DU12US
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION DDOT,DNRM2
      DIMENSION A(MDA,*),B(MDB,*),RNORM(*),H(*),W(*)
      INTEGER IC(*),IR(*)
C***FIRST EXECUTABLE STATEMENT  DU12US
      K=KRANK
      KP1=K+1
C
C        RANK=0
C
      IF(K.GT.0) GO TO 410
      DO 404 JB=1,NB
      RNORM(JB)=DNRM2(M,B(1,JB),1)
  404 CONTINUE
      DO 406 JB=1,NB
      DO 406 I=1,N
      B(I,JB)=0.0D0
  406 CONTINUE
      RETURN
C
C     REORDER B TO REFLECT ROW INTERCHANGES
C
  410 CONTINUE
      I=0
  412 I=I+1
      IF(I.EQ.M) GO TO 418
      J=IR(I)
      IF(J.EQ.I) GO TO 412
      IF(J.LT.0) GO TO 412
      IR(I)=-IR(I)
      DO 413 JB=1,NB
      RNORM(JB)=B(I,JB)
  413 CONTINUE
      IJ=I
  414 DO 415 JB=1,NB
      B(IJ,JB)=B(J,JB)
  415 CONTINUE
      IJ=J
      J=IR(IJ)
      IR(IJ)=-IR(IJ)
      IF(J.NE.I) GO TO 414
      DO 416 JB=1,NB
      B(IJ,JB)=RNORM(JB)
  416 CONTINUE
      GO TO 412
  418 CONTINUE
      DO 420 I=1,M
      IR(I)=ABS(IR(I))
  420 CONTINUE
C
C     IF A IS OF REDUCED RANK AND MODE=2,
C     APPLY HOUSEHOLDER TRANSFORMATIONS TO B
C
      IF(MODE.LT.2 .OR. K.EQ.M) GO TO 440
      MMK=M-K
      DO 430 JB=1,NB
      DO 425 J=1,K
      I=KP1-J
      TT=-DDOT(MMK,A(KP1,I),1,B(KP1,JB),1)/W(I)
      TT=TT-B(I,JB)
      CALL DAXPY(MMK,TT,A(KP1,I),1,B(KP1,JB),1)
      B(I,JB)=B(I,JB)+TT*W(I)
  425 CONTINUE
  430 CONTINUE
C
C     FIND NORMS OF RESIDUAL VECTOR(S)..(BEFORE OVERWRITE B)
C
  440 DO 442 JB=1,NB
      RNORM(JB)=DNRM2((M-K),B(KP1,JB),1)
  442 CONTINUE
C
C     BACK SOLVE LOWER TRIANGULAR L
C
      DO 450 JB=1,NB
      DO 448 I=1,K
      B(I,JB)=B(I,JB)/A(I,I)
      IF(I.EQ.K) GO TO 450
      IP1=I+1
      CALL DAXPY(K-I,-B(I,JB),A(IP1,I),1,B(IP1,JB),1)
  448 CONTINUE
  450 CONTINUE
C
C
C      TRUNCATED SOLUTION
C
      IF(K.EQ.N) GO TO 462
      DO 460 JB=1,NB
      DO 460 I=KP1,N
      B(I,JB)=0.0D0
  460 CONTINUE
C
C     APPLY HOUSEHOLDER TRANSFORMATIONS TO B
C
  462 DO 470 I=1,K
      J=KP1-I
      TT=A(J,J)
      A(J,J)=H(J)
      DO 465 JB=1,NB
      BB=-DDOT(N-J+1,A(J,J),MDA,B(J,JB),1)/H(J)
      CALL DAXPY(N-J+1,BB,A(J,J),MDA,B(J,JB),1)
  465 CONTINUE
      A(J,J)=TT
  470 CONTINUE
C
C
C     REORDER B TO REFLECT COLUMN INTERCHANGES
C
      I=0
  482 I=I+1
      IF(I.EQ.N) GO TO 488
      J=IC(I)
      IF(J.EQ.I) GO TO 482
      IF(J.LT.0) GO TO 482
      IC(I)=-IC(I)
  484 CALL DSWAP(NB,B(J,1),MDB,B(I,1),MDB)
      IJ=IC(J)
      IC(J)=-IC(J)
      J=IJ
      IF(J.EQ.I) GO TO 482
      GO TO 484
  488 CONTINUE
      DO 490 I=1,N
      IC(I)=ABS(IC(I))
  490 CONTINUE
C
C        SOLUTION VECTORS ARE IN FIRST N ROWS OF B(,)
C
      RETURN
      END


text: St. Patricks Day is ostensibly a religious holidays holiday on March 17th celebrating Saint Patrick, the patron saint of Ireland, and the arrival of Christianity in Ireland.  As celebrated in the United States, its either a cultural celebration or a alcohol boozefilled day of revelry.  Many restaurants and bars in town offer greenthemed specials.  The Graduate, for instance, has a long history of opening at 6am with green eggs and ham and green beer.


text:
      subroutine project(fcm,c,nat)
      implicit double precision (a-h,o-z)

C     uses the Miller,Handy and Adams algorithm
C     to project out rotational and translational
C     components of teh force constant matrix
C     (strickly only valid at a stationary point
C      but can be useful in geometry optimisations)

C     fcm = force constant matrix
C     C = cartesian geometry
C     Nat = number of atoms

C     Caution - overwrites original force constant matrix
C     fcm with projected version

      dimension c(3,nat),fcm(nat*3,nat*3)
c     temporary arrays
c     dimension f1(nat*3,nat*3),al(nat*3,nat*3)

      real*8, allocatable :: f1(:,:),al(:,:)

      nat3=nat*3
      allocate(f1(nat3,nat3))
      allocate(al(nat3,nat3))

      call zero(al,nat3*6)
      DO  N=1,NAT
      AL(3*N-2,1)=1.0D0
      AL(3*N-1,2)=1.0D0
      AL(3*N,3)=1.0D0
      AL(3*N-2,4)=-C(2,N)
      AL(3*N-1,4)=C(1,N)
      AL(3*N-1,5)=-C(3,N)
      AL(3*N,5)=C(2,N)
      AL(3*N-2,6)=C(3,N)
      AL(3*N,6)=-C(1,N)
      ENDDO

      call  orthy(al,nat3,6)
      call dgemm('N','T',nat3,nat3,6,-1d0,al,nat3,al,nat3,0d0,f1,nat3)

      do i=1,nat3
      f1(i,i)=f1(i,i)+1d0
      enddo

C
      call dgemm('N','N',nat3,nat3,nat3,1d0,f1,nat3,fcm,nat3,0d0,al,
     1 nat3)
      call dgemm('N','N',nat3,nat3,nat3,1d0,al,nat3,f1,nat3,0d0,fcm,
     1   nat3)

      RETURN
      END
      subroutine orthy(a,m,n)
      implicit double precision (a-h,o-z)
      dimension a(m,n), work(n)

c     orthonormalises the columns of matrix A
c     caution : intended as part of subroutine projec
c     not suitable for general use
c     sets small elements of A to zero
c     small aribitrarily defined as 1d-10

      parameter (small=1d-10)

      do j=1,n
      do i=1,m
      if(dabs(a(i,j)).lt.small)a(i,j)=0d0
      enddo
      enddo

      c=dsqrt(ddot(m,a(1,1),1,a(1,1),1))

      if(c.gt.small)then
       a(1:m,1)=a(1:m,1)/c
      else
       a(1:m,1)=0d0
      endif

      do na=2,n

          do nb=1,na-1
          work(nb)=ddot(m,a(1,na),1,a(1,nb),1)
          enddo

          do nb=1,na-1
          call daxpy(m,-work(nb),a(1,nb),1,a(1,na),1)
          enddo

          do i=1,m
          if(dabs(a(i,na)).lt.small)a(i,na)=0d0
          enddo

          c=dsqrt(ddot(m,a(1,na),1,a(1,na),1))
          if(c.gt.small)then
          a(1:m,na)=a(1:m,na)/c
          else
          a(1:m,na)=0d0
          endif

      enddo
      return
      end
      subroutine zero(a,n)
      double precision a(n)
      do i=1,n
      a(i)=0d0
      enddo
      return
      end

max_stars_repo_name: mickcollins/SMFA

max_issues_repo_name: mickcollins/SMFA

max_forks_repo_name: mickcollins/SMFAPAC


text:
!******************************************************************************
!******************************************************************************
MODULE recording_simulation

  !/*	external modules	*/

    USE shared_constants

  !/*	setup	*/

    IMPLICIT NONE

    PRIVATE

    PUBLIC :: record_simulation

    !/* explicit interface   */

    INTERFACE record_simulation

        MODULE PROCEDURE record_simulation_progress, record_simulation_start, record_simulation_stop

    END INTERFACE

CONTAINS
!******************************************************************************
!******************************************************************************
SUBROUTINE record_simulation_start(num_agents_sim, seed_sim)

    !/* external objects        */

    INTEGER(our_int), INTENT(IN)    :: num_agents_sim
    INTEGER(our_int), INTENT(IN)    :: seed_sim

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

    100 FORMAT(2x,A32,1x,i8,1x,A16,1x,i8)

    OPEN(UNIT=99, FILE='sim.respy.log', ACTION='WRITE')

        WRITE(99, 100) 'Starting simulation of model for', num_agents_sim, 'agents with seed', seed_sim
        WRITE(99, *)

    CLOSE(99)

END SUBROUTINE
!******************************************************************************
!******************************************************************************
SUBROUTINE record_simulation_stop()

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

    OPEN(UNIT=99, FILE='sim.respy.log', ACCESS='APPEND', ACTION='WRITE')

        WRITE(99, *) ' ... finished'
        WRITE(99, *)

    CLOSE(99)

END SUBROUTINE
!******************************************************************************
!******************************************************************************
SUBROUTINE record_simulation_progress(i)

    !/* external objects        */

    INTEGER(our_int), INTENT(IN)    :: i

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

    IF ((i .NE. zero_int) .AND. (MOD(i, 100) == zero_dble)) THEN

        100 FORMAT(A16,i10,A7)

        OPEN(UNIT=99, FILE='sim.respy.log', ACCESS='APPEND', ACTION='WRITE')

            WRITE(99, 100) ' ... simulated ', i, ' agents'
            WRITE(99, *)

        CLOSE(99)

    ELSE

        RETURN

    END IF

END SUBROUTINE
!******************************************************************************
!******************************************************************************
END MODULE

max_stars_repo_name: boryana-ilieva/respy

max_issues_repo_name: boryana-ilieva/respy

max_forks_repo_name: boryana-ilieva/respy


text:
      SUBROUTINE EXI2
C
C     EXI2 PERFORMS EXTERNAL FORMAT SOFIN OPERATIONS
C
      EXTERNAL LSHIFT
      LOGICAL  USRMSG
      INTEGER  DRY      ,UNAME    ,UNIT     ,SYSBUF   ,A        ,
     1         Z        ,SOF      ,PRC      ,Q4       ,T3       ,
     2         SRD      ,SWRT     ,EOI      ,SP       ,BAR      ,
     3         SCR1     ,SUBR(2)  ,BUF1     ,BUF2     ,BUF3     ,
     4         BUF4     ,HDR(7)   ,RC       ,MCB(7)   ,PREC     ,
     5         DIT      ,NAME(2)  ,EOG      ,PLTS     ,OFFSET
      REAL     ZERO(6)
      DOUBLE PRECISION   DZ(1)    ,DA
      CHARACTER          UFM*23   ,UWM*25   ,UIM*29
      COMMON  /XMSSG /   UFM      ,UWM      ,UIM
      COMMON  /MACHIN/   MACH
      COMMON  /BLANK /   DRY      ,X1(3)    ,UNAME(2) ,X2(18)   ,
     1                   UNIT     ,UNIVAC   ,LBUF     ,IADD
      COMMON  /SYSTEM/   SYSBUF   ,NOUT     ,X3(6)    ,NLPP     ,
     1                   X4(2)    ,LINE
      COMMON  /ZBLPKX/   A(4)     ,IROW
      COMMON  /NAMES /   RD       ,RDREW    ,WRT      ,WRTREW   ,
     1                   REW
      COMMON  /TYPE  /   PRC(2)   ,NWORD(4)
      COMMON  /ZZZZZZ/   Z(1)
      EQUIVALENCE       (Z(1),DZ(1))        ,(A(1),DA)
      DATA     SOF      ,SRD      ,SWRT     ,EOI      ,SP       /
     1         4HSOF    ,1        ,2        ,3        ,1        /
      DATA     LEOF     ,JH       ,SCR1     ,SUBR               /
     3         4H$EOF   ,1        ,301      ,4HEXI2   ,4H       /
      DATA     DIT      ,MDI      ,EOG      ,ZERO               /
     5         4HDIT    ,4HMDI    ,2        ,6*0.0              /
      DATA     Q4       ,T3       ,BAR      ,PLTS               /
     7         2HQ4     ,2HT3     ,2HBR     ,4HPLTS             /
C
C     INITIALIZE
C
      NCORE = KORSZ(Z)
      I     = NCORE - LBUF
      IF (MACH .EQ. 12) I = I - LBUF
      NCORE = I - 1
      IRW   = IADD
      IADD  = I
      CALL EXFORT (3,UNIT,0,0,IRW,0,0)
      BUF1  = NCORE - SYSBUF + 1
      BUF2  = BUF1  - SYSBUF - 1
      BUF3  = BUF2  - SYSBUF
      BUF4  = BUF3  - SYSBUF
      NCORE = BUF4  - 1
      NOS   = 0
      IDM   = 1
      USRMSG=.TRUE.
      LCORE = NCORE
      IF (NCORE .LE. 0) GO TO 9008
      CALL SOFOPN (Z(BUF1),Z(BUF2),Z(BUF3))
      CALL PAGE
C
C     READ THE HEADER OF THE NEXT ITEM AND FETCH THE ITEM ON THE SOF
C
   10 CALL EXFORT (SRD,UNIT,JH,HDR,7,SP,0)
   20 NAME(1) = HDR(1)
      NAME(2) = HDR(2)
      ITEM    = HDR(3)
      ITEST   = HDR(7)
      IF (ITEST .EQ. EOI) ITEST = EOG
      IF (HDR(1).EQ.DIT .OR. HDR(1).EQ. MDI) GO TO 200
      IF (HDR(3).EQ. -1 .OR. HDR(1).EQ.LEOF) GO TO 300
      ITM = ITTYPE(HDR(3))
      IF (ITM .EQ. 1) GO TO 100
      RC = 3
      CALL SFETCH (HDR(1),HDR(3),SWRT,RC)
      IF (RC .EQ. 3) GO TO 60
      LINE = LINE + 2
      IF (LINE .GT. NLPP) CALL PAGE
      GO TO (30,30,60,40,50), RC
   30 WRITE (NOUT,6346) UWM,HDR(1),HDR(2),HDR(3)
      USRMSG = .FALSE.
      GO TO 60
   40 CALL EXLVL (NOS,Z(IDM),HDR,Z,LCORE)
      RC = 3
      CALL SFETCH (HDR(1),HDR(3),SWRT,RC)
      IF (RC .EQ. 3) GO TO 60
   50 CALL SMSG (RC-2,HDR(3),HDR)
      USRMSG = .FALSE.
   60 CONTINUE
C
C     TABLES
C
C
C     ELSETS TABLE CORRECTION BY G.CHAN/UNISYS   4/91
C
C     IN 91 VERSION, ELEMENT PLOT SYMBOL LINE HAS 2 WORDS, SYMBOL AND
C     NO. OF GRID POINT PER ELEMENT, NGPEL, WRITTEN OUT BY EXO2 USING
C     FORMAT 25. THE ELSETS DATA LINE COMING UP NEXT USE FORMAT 10 FOR
C     ELEMENTS WITH NO OFFSETS, FORMAT 26 FOR BAR WHICH HAS 6 OFFSET
C     VALUES, AND FORMATS 27 AND 28 FOR TRIA3 AND QUAD4 WHICH HAS 1
C     OFFSET VALUE EACH.
C     IN 90 AND EARLIER VERSIONS, ONLY ONE ELEMENT PLOT SYMBOL WORD WAS
C     WRITTEN OUT, AND ON ELSETS DATA LINE COMING UP NEXT, FORMAT 10
C     WAS USED FOR ALL ELEMENTS. NO OFFSET DATA FOR THE BAR, QUAD4 AND
C     TRIA3 ELEMENTS. NGPEL WAS THE FIRST WORD ON THE ELSETS DATA LINE.
C     ALSO, THE 90 AND EARLIER VERSIONS DID NOT COUNT PROPERTY ID, PID,
C     ON THE ELSETS DATA LINE. THUS THE TOTAL NO. OF WORDS MAY BE IN
C     ERROR AND MAY CAUSE EXTRA ZEROS AT THE END OF THE DATA LINE.
C
C     THEREFORE, IF THE 90 OR EARLIER EXTERNAL SOF FILE WAS USED, WE
C     NEED TO ADD THE OFFSETS (1 OR 6 FLOATING POINTS ZEROS) TO THE BAR,
C     QUAD4 AND TRIA3 ELEMENTS FOR THE ELSETS TABLE.
C     (AS OF 4/91, THESE CHANGES HAVE NOT BEEN TESTED)
C
      OFFSET = 0
   70 NWDS = HDR(5)
      IF (NWDS .GT. LCORE) GO TO 9008
      CALL EXFORT (SRD,UNIT,HDR(4),Z,NWDS,SP,0)
      IF (OFFSET .EQ. 0) GO TO 80
      J = 1
      CALL SUWRT (Z(1),1,J)
      NP2 = Z(1) + 2
      DO 73 K = 2,NWDS,NP2
      IF (Z(K) .EQ. 0) GO TO 75
      CALL SUWRT (Z(K),NP2,J)
      CALL SUWRT (ZERO,OFFSET,J)
   73 CONTINUE
   75 Z(1) = 0
      NWDS = 1
   80 CALL SUWRT (Z,NWDS,ITEST)
      IF (HDR(7) .EQ. EOI) GO TO 90
      CALL EXFORT (SRD,UNIT,JH,HDR,7,SP,0)
      IF (HDR(1).NE.NAME(1) .OR. HDR(2).NE.NAME(2)) GO TO 160
      IF (HDR(3) .NE. ITEM) GO TO 160
      ITEST = HDR(7)
      IF (ITEST .EQ. EOI) ITEST = EOG
      IF (ITEM.NE.PLTS .OR. HDR(5).NE.1 .OR. HDR(4).NE.10) GO TO 85
      OFFSET = 0
      IF (Z(1) .EQ. BAR) OFFSET = 6
      IF (Z(1).EQ.Q4 .OR. Z(1).EQ.T3) OFFSET = 1
   85 IF (HDR(4) .GT. 0) GO TO 70
   90 ITEST = EOI
      CALL SUWRT (0,0,ITEST)
      GO TO 140
C
C     MATRICES
C
C
C     READ TRAILER
C
  100 CALL SOFTRL (HDR(1),HDR(3),MCB(1))
      RC = MCB(1)
      IF (RC .EQ. 3) GO TO 108
      LINE = LINE + 2
      IF (LINE .GT. NLPP) CALL PAGE
      GO TO (102,102,108,104,106), RC
  102 WRITE (NOUT,6346) UWM,HDR(1),HDR(2),HDR(3)
      USRMSG = .FALSE.
      GO TO 108
  104 CALL EXLVL (NOS,Z(IDM),HDR,Z,LCORE)
      GO TO 108
  106 CALL SMSG (3,HDR(3),HDR)
      USRMSG = .FALSE.
  108 CALL EXFORT (SRD,UNIT,HDR(4),MCB(2),6,SP,0)
      NCOL   = MCB(2)
      PREC   = MCB(5)
      MCB(1) = SCR1
      MCB(2) = 0
      MCB(6) = 0
      MCB(7) = 0
      IF (USRMSG) CALL GOPEN (SCR1,Z(BUF4),WRTREW)
C
C     READ MATRIX ONE COLUMN AT A TIME AND PACK ON SCR2
C
      DO 130 J = 1,NCOL
      CALL EXFORT (SRD,UNIT,JH,HDR,7,SP,0)
      IF (HDR(1).NE.NAME(1) .OR. HDR(2).NE.NAME(2)) GO TO 160
      IF (HDR(3) .NE. ITEM) GO TO 160
      NWDS = HDR(5)
      IF (NWDS*1.4 .GT. NCORE) GO TO 9008
      CALL EXFORT (SRD,UNIT,HDR(4),Z,NWDS,PREC,DZ)
      IF (.NOT. USRMSG) GO TO 130
      CALL BLDPK (PREC,PREC,SCR1,0,0)
      IPRC = PRC(PREC)
      N    = NWORD(PREC) + IPRC
      K    = 1
  110 IF (Z(K) .LT. 0) GO TO 120
      IROW = Z(K)
      A(1) = Z(K+IPRC)
      IF (PREC .EQ. 1) GO TO 115
      A(2) = Z(K+IPRC+1)
      IF (PREC .LE. 3) GO TO 115
      A(3) = Z(K+4)
      A(4) = Z(K+5)
  115 CALL ZBLPKI
      K = K + N
      GO TO 110
  120 CALL BLDPKN (SCR1,0,MCB)
  130 CONTINUE
      IF (.NOT.USRMSG) GO TO 150
      CALL WRTTRL (MCB)
      CALL CLOSE  (SCR1,REW)
      CALL MTRXO  (SCR1,HDR,HDR(3),0,RC)
C
C     WRITE USER MESSAGE
C
  140 IF (.NOT.USRMSG) GO TO 150
      LINE = LINE + 1
      IF (LINE .GT. NLPP) CALL PAGE
      WRITE (NOUT,6357) UIM,HDR(1),HDR(2),HDR(3),UNAME,SOF
  150 USRMSG = .TRUE.
      GO TO 10
C
C     NO EOI FOR ITEM AND A NEW ITEM WAS READ
C
  160 LINE = LINE + 2
      IF (LINE .GT. NLPP) CALL PAGE
      WRITE (NOUT,6363) UWM,NAME(1),NAME(2),ITEM,UNAME
      IF (ITM .EQ. 0) CALL DELETE (NAME,ITEM,RC)
      IF (ITM .EQ. 1) CALL CLOSE (SCR1,REW)
      USRMSG = .TRUE.
      GO TO 20
C
C     READ DIT AND MDI
C
  200 NOS   = HDR(5)/2
      LCORE = NCORE - HDR(5)*4
      IDM   = LCORE + 1
      IF (6*NOS .GT. LCORE) GO TO 9008
      CALL EXFORT (SRD,UNIT,HDR(4),Z,HDR(5),SP,0)
      DO 210 I = 1,NOS
      Z(IDM+4*I-4) = Z(2*I-1)
      Z(IDM+4*I-3) = Z(2*I  )
  210 CONTINUE
      CALL EXFORT (SRD,UNIT,JH,HDR,7,SP,0)
      CALL EXFORT (SRD,UNIT,HDR(4),Z,HDR(5),SP,0)
      DO 220 I = 1,NOS
      J = IDM + 4*I - 2
      K = 6*I - 6
      Z(J  ) = LSHIFT(Z(K+1),20) + LSHIFT(Z(K+2),10) + Z(K+3)
      Z(J+1) = LSHIFT(Z(K+4),20) + LSHIFT(Z(K+5),10) + Z(K+6)
  220 CONTINUE
      GO TO 10
C
C     NORMAL MODULE COMPLETION
C
  300 CALL SOFCLS
      RETURN
C
C     ABNORMAL MODULE COMPLETION
C
 9008 CALL MESAGE (8,0,SUBR)
      DRY = -2
      CALL SOFCLS
      RETURN
C
C     MESSAGE TEXTS
C
 6346 FORMAT (A25,' 6346, SUBSTRUCTURE ',2A4,' ITEM ',A4,
     1       ' NOT COPIED.  IT ALREADY EXISTS ON THE SOF.')
 6357 FORMAT (A29,' 6357, SUBSTRUCTURE ',2A4,' ITEM ',A4,
     1       ' SUCCESSFULLY COPIED FROM ',2A4,' TO ',A4)
 6363 FORMAT (A25,' 6363, INCOMPLETE DATA FOR SUBSTRUCTURE ',2A4,
     1       ' ITEM ',A4,' ON ',2A4,'. THE ITEM WILL NOT BE COPIED.')
      END


text:
!> \brief contains the subprograms needed to build the hamiltonian and some forces
!> \author Alin M Elena
!> \date 07/11/07, 10:09:17
!
module m_Hamiltonian
  use m_Constants
  use m_Types
  use m_Useful
  use m_LinearAlgebra
  use m_SlaterKoster
  use m_Gutenberg, only: PrintMatrix
  use m_TightBinding
  use m_DensityMatrix
  private
!
  public :: BuildHamiltonian
  public :: AddBias
  public :: DiagHamiltonian
  public :: ZeroForces
  public :: RepulsiveForces
  public :: ElectronicForces
  public :: ElectronicEnergy
  public :: RepulsiveEnergy
contains
!
!> \brief builds the hamiltonian
!> \author Cristian G Sanchez, Alin M Elena
!> \date 07/11/07, 10:09:17
!> \param io type(ioType) contains all the info about I/O files
!> \param gen type(generalType) contains the info needed by the program to k_run
!> \param atomic type(atomicxType) contains all info about the atoms and basis set and some parameters
!> \param tbMod type(modelType) contains information about the tight binding model parameters
!> \param sol type(solutionType) contains information about the solution space
!
!
  subroutine BuildHamiltonian (io, gen, atomic, tbMod, sol)
    character (len=*), parameter :: myname = "BuildHamiltonian"
    type (ioType), intent (inout) :: io
    type (generalType), intent (inout) :: gen
    type (atomicxType), intent (inout) :: atomic
    type (modelType), intent (inout) :: tbMod
    type (solutionType), intent (inout) :: sol
    integer :: i, j, k, o, norbsi, norbsj
    real (kind=k_pr) :: rij, hij, l, m, n
!-------------------------------------------------!
!
    call ZeroMatrix (sol%h, io)
    if ( .not. gen%spin) then
!!!!!   !$OMP PARALLEL DO DEFAULT(shared) PRIVATE(i,j,k,o,hij,rij,l,m,n)  SCHEDULE(static)
      do i = 1, atomic%atoms%natoms - 1
        do j = i + 1, atomic%atoms%natoms
! Offsite terms
          call AtomDistance (atomic%atoms, j, i, rij, l, m, n)
          do k = 1, atomic%species%norbs(atomic%atoms%sp(i))
            do o = 1, atomic%species%norbs(atomic%atoms%sp(j))
              hij = hmn (rij, l, m, n, atomic%basis%orbitals(atomic%atoms%orbs(i, k)), atomic%basis%orbitals(atomic%atoms%orbs(j, &
             & o)), gen, tbMod, sol)
              if (Abs(hij) >= gen%hElementThreshold) then
                call SpmPut (sol%h, atomic%atoms%orbs(i, k), atomic%atoms%orbs(j, o), cmplx(hij, 0.0_k_pr, k_pr))
                call SpmPut (sol%h, atomic%atoms%orbs(j, o), atomic%atoms%orbs(i, k), cmplx(hij, 0.0_k_pr, k_pr))
              end if
            end do
          end do
        end do
!            Onsite terms
        do k = 1, atomic%species%norbs(atomic%atoms%sp(i))
          do o = 1, atomic%species%norbs(atomic%atoms%sp(i))
            hij = Onsite (atomic%basis%orbitals(atomic%atoms%orbs(i, k)), atomic%basis%orbitals(atomic%atoms%orbs(i, o)), tbMod)
            call SpmPut (sol%h, atomic%atoms%orbs(i, k), atomic%atoms%orbs(i, o), cmplx(hij, 0.0_k_pr, k_pr))
          end do
        end do
      end do
! this may look silly but the standard does not specify the value of i at the end of the loop
      i = atomic%atoms%natoms
      do k = 1, atomic%species%norbs(atomic%atoms%sp(i))
        do o = 1, atomic%species%norbs(atomic%atoms%sp(i))
          hij = Onsite (atomic%basis%orbitals(atomic%atoms%orbs(i, k)), atomic%basis%orbitals(atomic%atoms%orbs(i, o)), tbMod)
          call SpmPut (sol%h, atomic%atoms%orbs(i, k), atomic%atoms%orbs(i, o), cmplx(hij, 0.0_k_pr, k_pr))
        end do
      end do
!!!!!   !$OMP END PARALLEL DO
    else if (gen%collinear) then
      do i = 1, atomic%atoms%natoms - 1
        norbsi = atomic%species%norbs(atomic%atoms%sp(i)) / 2
        do j = i + 1, atomic%atoms%natoms
          norbsj = atomic%species%norbs(atomic%atoms%sp(j)) / 2
! Offsite terms
          call AtomDistance (atomic%atoms, j, i, rij, l, m, n)
          do k = 1, norbsi
            do o = 1, norbsj
              hij = hmn (rij, l, m, n, atomic%basis%orbitals(atomic%atoms%orbs(i, k)), atomic%basis%orbitals(atomic%atoms%orbs(j, &
             & o)), gen, tbMod, sol)
              if (Abs(hij) >= gen%hElementThreshold) then
                call SpmPut (sol%h, atomic%atoms%orbs(i, k), atomic%atoms%orbs(j, o), cmplx(hij, 0.0_k_pr, k_pr))
                call SpmPut (sol%h, atomic%atoms%orbs(j, o), atomic%atoms%orbs(i, k), cmplx(hij, 0.0_k_pr, k_pr))
              end if
              hij = hmn (rij, l, m, n, atomic%basis%orbitals(atomic%atoms%orbs(i, k+norbsi)), &
             & atomic%basis%orbitals(atomic%atoms%orbs(j, o+norbsj)), gen, tbMod, sol)
              if (Abs(hij) >= gen%hElementThreshold) then
                call SpmPut (sol%h, atomic%atoms%orbs(i, k+norbsi), atomic%atoms%orbs(j, o+norbsj), cmplx(hij, 0.0_k_pr, k_pr))
                call SpmPut (sol%h, atomic%atoms%orbs(j, o+norbsj), atomic%atoms%orbs(i, k+norbsi), cmplx(hij, 0.0_k_pr, k_pr))
              end if
            end do
          end do
        end do
!            Onsite terms
        do k = 1, norbsi
          do o = 1, norbsi
            hij = Onsite (atomic%basis%orbitals(atomic%atoms%orbs(i, k)), atomic%basis%orbitals(atomic%atoms%orbs(i, o)), tbMod)
            call SpmPut (sol%h, atomic%atoms%orbs(i, k), atomic%atoms%orbs(i, o), cmplx(hij, 0.0_k_pr, k_pr))
            hij = Onsite (atomic%basis%orbitals(atomic%atoms%orbs(i, k+norbsi)), atomic%basis%orbitals(atomic%atoms%orbs(i, &
           & o+norbsi)), tbMod)
            call SpmPut (sol%h, atomic%atoms%orbs(i, k+norbsi), atomic%atoms%orbs(i, o+norbsi), cmplx(hij, 0.0_k_pr, k_pr))
          end do
        end do
      end do
      i = atomic%atoms%natoms
      norbsi = atomic%species%norbs(atomic%atoms%sp(i)) / 2
      do k = 1, norbsi
        do o = 1, norbsi
          hij = Onsite (atomic%basis%orbitals(atomic%atoms%orbs(i, k)), atomic%basis%orbitals(atomic%atoms%orbs(i, o)), tbMod)
          call SpmPut (sol%h, atomic%atoms%orbs(i, k), atomic%atoms%orbs(i, o), cmplx(hij, 0.0_k_pr, k_pr))
          hij = Onsite (atomic%basis%orbitals(atomic%atoms%orbs(i, k+norbsi)), atomic%basis%orbitals(atomic%atoms%orbs(i, &
         & o+norbsi)), tbMod)
          call SpmPut (sol%h, atomic%atoms%orbs(i, k+norbsi), atomic%atoms%orbs(i, o+norbsi), cmplx(hij, 0.0_k_pr, k_pr))
        end do
      end do
!
!
    else
      call error ("Non-Collinear spins are not implemented yet!", myname, .true., io)
    end if
  end subroutine BuildHamiltonian
!
!
!
!> \brief adds an external field to the diagonal elements of the Hamiltonian                                 !
!> \author Cristian G. Sanchez
!> \date ~2005
!> \param atomic type(atomicType) contains all info about the atoms and basis set and some parameters
!> \param sol type(solutionType) contains information about the solution space
!> \param factor real multiplies the bias
!
  subroutine AddBias (factor, atomic, sol)
    character (len=*), parameter :: myname = 'AddBias'
    real (k_pr), intent (in) :: factor
    type (atomicxType), intent (in) :: atomic
    type (solutionType), intent (inout) :: sol
    integer :: i, k
!
    do i = 1, atomic%atoms%natoms
      do k = 1, atomic%species%norbs(atomic%atoms%sp(i))
        sol%h%a (atomic%atoms%orbs(i, k), atomic%atoms%orbs(i, k)) = sol%h%a(atomic%atoms%orbs(i, k), atomic%atoms%orbs(i, k)) + &
       & cmplx (factor*atomic%atoms%bias(i), 0.0_k_pr, k_pr)
      end do
    end do
!
  end subroutine AddBias
!
!> \brief diagonalizes the hamiltonian
!> \author Cristian G. Sanchez, Alin M Elena
!> \date 07/11/07, 13:04:42
!> \param io type(ioType) contains all the info about I/O files
!> \param gen type(generalType) contains the info needed by the program to k_run
!> \param sol type(solutionType) contains information about the solution space
!> \param atomic type(atomicxType) contains all info about the atoms and basis set and some parameters
!
  subroutine DiagHamiltonian (io, gen, atomic, sol)
    character (len=*), parameter :: myname = 'DiagHamiltonian'
    type (ioType), intent (inout) :: io
    type (generalType), intent (inout) :: gen
    type (solutionType), intent (inout) :: sol
    type (atomicxType), intent (inout) :: atomic
    integer :: n, ns, i, j
    if (gen%spin) then
      n = sol%h%dim
      ns = sol%hup%dim
      sol%eigenvals (1:n) = 0.0_k_pr
      call ZeroMatrix (sol%eigenvecs, io)
      do i = 1, ns
        do j = 1, ns
          sol%hdown%a (i, j) = sol%h%a(i, j)
          sol%buff%h%a(i, j) = sol%h%a(i, j)
        end do
      end do
      sol%buff%tmpA = 0.0_k_pr
      call ZeroMatrix (sol%buff%tmpB, io)
      call DiagonalizeMatrix (sol%buff%h, sol%buff%tmpB, sol%buff%tmpA, io)
      sol%eigenvals (1:ns) = sol%buff%tmpA(1:ns)
      do i = 1, ns
        do j = 1, ns
          sol%eigenvecs%a (i, j) = sol%buff%tmpB%a(i, j)
          sol%hup%a (i, j) = sol%h%a(i+ns, j+ns)
          sol%buff%h%a(i, j) = sol%h%a(i+ns, j+ns)
        end do
      end do
      sol%buff%tmpA = 0.0_k_pr
      call ZeroMatrix (sol%buff%tmpB, io)
      call DiagonalizeMatrix (sol%buff%h, sol%buff%tmpB, sol%buff%tmpA, io)
      sol%eigenvals (1+ns:n) = sol%buff%tmpA(1:ns)
      do i = 1, ns
        do j = 1, ns
          sol%eigenvecs%a (i+ns, j+ns) = sol%buff%tmpB%a(i, j)
        end do
      end do
!
      if (gen%lIsExcited) then
        call CreateDensityMatrixExcited (gen, atomic, sol, io)
      else
        call CreateDensityMatrixSpin (gen, atomic, sol, io)
      end if

    else
      sol%eigenvals = 0.0_k_pr
      call ZeroMatrix (sol%eigenvecs, io)
      call CopyMatrix (sol%buff%h, sol%h, io)
      call DiagonalizeMatrix (sol%h, sol%eigenvecs, sol%eigenvals, io)
      call CreateDensityMatrixNoSpin (gen, atomic, sol, io)
    end if
!
  end subroutine DiagHamiltonian
!
!> \brief resets the forces
!> \author Alin M Elena
!> \date 07/11/07, 16:43:12
!> \param atomic type(atomicxType) info about the atoms
!
  subroutine ZeroForces (atomic)
    character (len=*), parameter :: sMyName = "ZeroForces"
    type (atomicxType), intent (inout) :: atomic
!
    atomic%atoms%fx = 0.0_k_pr
    atomic%atoms%fy = 0.0_k_pr
    atomic%atoms%fz = 0.0_k_pr
!
  end subroutine ZeroForces
!
!> \brief computes the contribution to forces by the repulsive potential
!> \author Alin M Elena
!> \date 07/11/07, 17:20:31
!> \param atomic type(atomicType) info about the atoms
!> \param gen type(generalType) contains the info needed by the program to k_run
!> \param tb type(modelType) contains information about the tight binding model parameters
  subroutine RepulsiveForces (gen, atomic, tb)
!--subroutine name--------------------------------!
    character (len=*), parameter :: myname = 'RepulsiveForces'
    type (atomicType), intent (inout) :: atomic
    type (generalType), intent (inout) :: gen
    type (modelType), intent (inout) :: tb
    real (k_pr) :: rij, rbr, rib !, rir
    real (k_pr) :: rxij, ryij, rzij
    real (k_pr) :: rxbr, rybr, rzbr
    real (k_pr) :: rxib, ryib, rzib
!    real(k_pr) :: rxir, ryir, rzir
    integer :: i, j, b, r, k
    real (k_pr) :: fact, fact1
!-------------------------------------------------!
!
    if ( .not. gen%embedding) then
!
      do k = 1, atomic%nmoving
        i = atomic%moving (k)
        do j = 1, atomic%natoms
          if (i /= j) then
            call AtomDistVec (atomic, i, j, rij, rxij, ryij, rzij)
            fact = - RepP (rij, atomic%sp(i), atomic%sp(j), tb, gen) / rij
            atomic%fx (i) = atomic%fx(i) + fact * rxij
            atomic%fy (i) = atomic%fy(i) + fact * ryij
            atomic%fz (i) = atomic%fz(i) + fact * rzij
          end if
        end do
      end do
    else ! now with embedding
      do k = 1, atomic%nmoving
        b = atomic%moving (k)
        do i = 1, atomic%natoms
          if (i /= b) then
            call AtomDistVec (atomic, i, b, rib, rxib, ryib, rzib)
            fact = - RepP (rib, atomic%sp(i), atomic%sp(b), tb, gen) * EmbeddingP (argument(i), atomic%sp(b), tb) / rib
            atomic%fx (b) = atomic%fx(b) - fact * rxib
            atomic%fy (b) = atomic%fy(b) - fact * ryib
            atomic%fz (b) = atomic%fz(b) - fact * rzib
          end if
        end do
        fact1 = EmbeddingP (argument(b), atomic%sp(b), tb)
        do r = 1, atomic%natoms
          if (b /= r) then
            call AtomDistVec (atomic, b, r, rbr, rxbr, rybr, rzbr)
            fact = - fact1 * RepP (rbr, atomic%sp(b), atomic%sp(r), tb, gen) / rbr
            atomic%fx (b) = atomic%fx(b) + fact * rxbr
            atomic%fy (b) = atomic%fy(b) + fact * rybr
            atomic%fz (b) = atomic%fz(b) + fact * rzbr
          end if
        end do
      end do
    end if
!
  contains
    function argument (m)
      integer :: m
      real (k_pr) :: argument
      integer :: i
      real (k_pr) :: sum, rim
!
      sum = 0.0_k_pr
      do i = 1, atomic%natoms
        if (i /= m) then
          rim = Distance (atomic, i, m)
          sum = sum + Rep (rim, atomic%sp(i), atomic%sp(m), tb, gen)
        end if
      end do
      argument = sum
    end function argument
  end subroutine RepulsiveForces
!
!> \brief calculates force operator
!> \author Cristian G Sanchez
!> \date ~2005
!> \param alpha integer the direction 1-2-3 ->x-y-z
!> \param j integer the atom
!> \param gen type(generalType) contains the info needed by the program to k_run
!> \param atomic type(atomicType) contains all info about the atoms and basis set and some parameters
!> \param tb type(modelType) contains information about the tight binding model parameters
!> \param sol type(solutionType) contains information about the solution space
!> \param io type(ioType) contains all the info about I/O files
  subroutine ForceOperator (j, atomic, gen, tb, sol, io)
!--subroutine name--------------------------------!
    character (len=*), parameter :: myname = 'ForceOperator'
!--subroutine parameters -------------------------!
    integer, intent (inout) :: j
    type (generalType), intent (inout) :: gen
    type (atomicxType), intent (inout) :: atomic
    type (modelType), intent (inout) :: tb
    type (solutionType), intent (inout) :: sol
    type (ioType), intent (inout) :: io
!--internal variables ----------------------------!
    real (k_pr) :: rij
    integer :: i, k, o
    real (k_pr) :: fact (3), l, m, n !,fact2,ff2,
!-------------------------------------------------!
!
    call ZeroMatrix (sol%forceOpX, io)
    call ZeroMatrix (sol%forceOpY, io)
    call ZeroMatrix (sol%forceOpZ, io)
    do i = 1, atomic%atoms%natoms
      if (i /= j) then
        call AtomDistance (atomic%atoms, j, i, rij, l, m, n)
        do k = 1, atomic%species%norbs(atomic%atoms%sp(i))
          do o = 1, atomic%species%norbs(atomic%atoms%sp(j))
            fact = - DhmnXYZ (rij, l, m, n, atomic%basis%orbitals(atomic%atoms%orbs(i, k)), &
           & atomic%basis%orbitals(atomic%atoms%orbs(j, o)), gen, tb, sol)
            call SpmPut (sol%forceOpX, atomic%atoms%orbs(i, k), atomic%atoms%orbs(j, o), cmplx(fact(1), 0.0_k_pr, k_pr))
            call SpmPut (sol%forceOpX, atomic%atoms%orbs(j, o), atomic%atoms%orbs(i, k), cmplx(fact(1), 0.0_k_pr, k_pr))
            call SpmPut (sol%forceOpY, atomic%atoms%orbs(i, k), atomic%atoms%orbs(j, o), cmplx(fact(2), 0.0_k_pr, k_pr))
            call SpmPut (sol%forceOpY, atomic%atoms%orbs(j, o), atomic%atoms%orbs(i, k), cmplx(fact(2), 0.0_k_pr, k_pr))
            call SpmPut (sol%forceOpZ, atomic%atoms%orbs(i, k), atomic%atoms%orbs(j, o), cmplx(fact(3), 0.0_k_pr, k_pr))
            call SpmPut (sol%forceOpZ, atomic%atoms%orbs(j, o), atomic%atoms%orbs(i, k), cmplx(fact(3), 0.0_k_pr, k_pr))
          end do
        end do
      end if
    end do
  end subroutine ForceOperator
!
!> \brief adds electronic contribution to forces
!> \author Alin M Elena
!> \date 07/11/07, 18:30:30
!> \param io type(ioType) contains all the info about I/O files
!> \param gen type(generalType) contains the info needed by the program to k_run
!> \param atomic type(atomicxType) contains all info about the atoms and basis set and some parameters
!> \param tb type(modelType) contains information about the tight binding model parameters
!> \param sol type(solutionType) contains information about the solution space
  subroutine ElectronicForces (atomic, gen, tb, sol, io)
    character (len=*), parameter :: myname = 'ElectronicForces'
    type (generalType), intent (inout) :: gen
    type (atomicxType), intent (inout) :: atomic
    type (modelType), intent (inout) :: tb
    type (solutionType), intent (inout) :: sol
    type (ioType), intent (inout) :: io
    integer :: i, k
!
!
    do k = 1, atomic%atoms%nmoving
      i = atomic%atoms%moving(k)
      call ForceOperator (i, atomic, gen, tb, sol, io)
      atomic%atoms%fx (i) = atomic%atoms%fx(i) + ProductTrace (sol%rho, sol%forceOpX, io)
      atomic%atoms%fy (i) = atomic%atoms%fy(i) + ProductTrace (sol%rho, sol%forceOpY, io)
      atomic%atoms%fz (i) = atomic%atoms%fz(i) + ProductTrace (sol%rho, sol%forceOpZ, io)
    end do
!
  end subroutine ElectronicForces
!
!> \brief computes the electronic energy
!> \author Alin M Elena
!> \date 07/11/07, 23:15:07
!> \param io type(ioType) contains all the info about I/O files
!> \param gen type(generalType) contains the info needed by the program to k_run
!> \param sol type(solutionType) contains information about the solution space
!> \remarks \f[ E_e=Tr(\rho {\mathbf H})\f]
  real (k_pr) function ElectronicEnergy (gen, sol, io)
    character (len=*), parameter :: myname = 'ElectronicEnergy'
    type (solutionType), intent (inout) :: sol
    type (ioType), intent (inout) :: io
    type (generalType), intent (in) :: gen
!-------------------------------------------------!
!
!
    ElectronicEnergy = ProductTrace (sol%rho, sol%h, io)
!
  end function ElectronicEnergy
!
!> \brief computes repulsive energy
!> \author Alin M Elena
!> \date 07/11/07, 23:25:42
!> \param gen type(generalType) contains the info needed by the program to k_run
!> \param atomic type(atomicType) contains all info about the atoms and basis set and some parameters
!> \param tb type(modelType) contains information about the tight binding model parameters
  real (k_pr) function RepulsiveEnergy (gen, atomic, tb)
!--subroutine name--------------------------------!
    character (len=*), parameter :: myname = 'RepulsiveEnergy'
    real (k_pr) :: renergy, phi
    real (k_pr) :: rij
    type (modelType), intent (inout) :: tb
    type (generalType), intent (inout) :: gen
    type (atomicType), intent (inout) :: atomic
    integer :: i, j
!
    renergy = 0.0_k_pr
    if ( .not. gen%embedding) then
      do i = 1, atomic%natoms - 1
        do j = i + 1, atomic%natoms
          rij = Distance (atomic, i, j)
          renergy = renergy + Rep (rij, atomic%sp(i), atomic%sp(j), tb, gen)
        end do
      end do
    else
      do i = 1, atomic%natoms
        phi = 0.0_k_pr
        do j = 1, atomic%natoms
          if (i /= j) then
            rij = Distance (atomic, i, j)
            phi = phi + Rep (rij, atomic%sp(i), atomic%sp(j), tb, gen)
          end if
        end do
        renergy = renergy + embedding (phi, atomic%sp(i), tb)
      end do
    end if
    RepulsiveEnergy = renergy
!
  end function RepulsiveEnergy
!
end module m_Hamiltonian

max_stars_repo_name: naumruso/tdtbuj

max_issues_repo_name: naumruso/tdtbuj

max_forks_repo_name: naumruso/tdtbuj


text:
!$acc parallel loop reduction (max: error)
!$acc parallel loop

max_stars_repo_name: passlab/accparser

max_issues_repo_name: passlab/accparser

max_forks_repo_name: passlab/accparser


text:
!-----------------------------------------------------------------------
!     Module:        pies_fieldlines
!     Authors:       S. Lazerson (lazerson@pppl.gov)
!     Date:          11/8/2011
!     Description:   This module contains the variables utilized by the
!                    field line following routines.
!-----------------------------------------------------------------------
      MODULE pies_fieldlines
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stel_kinds
      USE EZspline_obj
      USE EZspline
!-----------------------------------------------------------------------
!     Module Variables
!          last_line    Last surface on which field lines are followed
!          hitsrf       First surface to hit the wall
!          nfollow
!          nintw        Total number of steps in phi for field line following
!          ninteqs      Size of q helper array (number of equations, two for each surface)
!          nu_spline    Number of polodial points in spline
!          nv_spline    Number of toroidal points in spline
!          hitwal       Array indicating if a field line has hit the wall
!          follow_tol   Tollerance for field line following
!          dkmin        Smallest distance to be resolved in Fourier Space
!          nfollow      Number of step in phi
!          magaxis_abs  Percision of magnetic axis solver
!          magaxis_eta  Value of |f(x)| is taken to be a zero
!          r_axis       Location of magnetic axis in phi=0 plane (in background coords)
!          z_axis       Location of magnetic axis in phi=0 plane (in background coords)
!          dphi         Change in phi when following field lines
!          psjmnc       Psi-Jacobian (1/B^v) (cos)
!          psjmns       Psi-Jacobian (1/B^v) (sin)
!          brho         B^s/B^v
!          btheta       B^u/B^v
!          rhobth       rho*B^u/B^v
!          bxsi_array   Array for field line follower
!          beta_array   Array for field line follower
!          rho_spl      Spline for rho
!          theta_spl    Spline for theta
!          beta_spl     Spline for brho
!          bxsi_spl     Spline for rhobth
!          xsiln        XSI of field line (r coordinate)
!          etaln        ETA of field line (z coordinate)
!
!----------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER                  :: last_line, hitsrf, nfollow, nintw, &
                                  ninteqs, nu_spline, nv_spline, ik_help
      INTEGER                  :: iuser(3)
      INTEGER, ALLOCATABLE     :: hitwal(:)
      REAL                     :: follow_tol
      REAL                     :: dkmin
      REAL                     :: magaxis_abs
      REAL                     :: magaxis_eta
      REAL                     :: r_axis
      REAL                     :: z_axis
      REAL                     :: dphi
      REAL                     :: phgaus
      REAL                     :: tfrac
      REAL(rprec), ALLOCATABLE :: psjmnc(:,:)
      REAL(rprec), ALLOCATABLE :: psjmns(:,:)
      REAL(rprec), ALLOCATABLE :: bthmnc(:,:)
      REAL(rprec), ALLOCATABLE :: bthmns(:,:)
      REAL(rprec), ALLOCATABLE :: brhomnc(:,:)
      REAL(rprec), ALLOCATABLE :: brhomns(:,:)
      REAL(rprec), ALLOCATABLE :: rhobth(:,:)
      REAL(rprec), ALLOCATABLE :: R_array(:)
      REAL(rprec), ALLOCATABLE :: Z_array(:)
      REAL(rprec), ALLOCATABLE :: bxsi_array(:)
      REAL(rprec), ALLOCATABLE :: beta_array(:)
      REAL(rprec), ALLOCATABLE :: Rln(:,:)
      REAL(rprec), ALLOCATABLE :: Zln(:,:)
      REAL(rprec), ALLOCATABLE :: rholn(:,:)
      REAL(rprec), ALLOCATABLE :: thetaln(:,:)
      REAL(rprec), ALLOCATABLE :: philn(:,:)
      REAL(rprec), ALLOCATABLE :: xsiln(:,:)
      REAL(rprec), ALLOCATABLE :: etaln(:,:)
      TYPE(EZspline3_r8)       :: brho_spl
      TYPE(EZspline3_r8)       :: theta_spl
      TYPE(EZspline3_r8)       :: beta_spl
      TYPE(EZspline3_r8)       :: bxsi_spl
      TYPE(EZspline3_r8)       :: R_spl
      TYPE(EZspline3_r8)       :: Z_spl
      TYPE(EZspline3_r8)       :: BS_spl
      TYPE(EZspline3_r8)       :: BU_spl
      TYPE(EZspline3_r8)       :: BV_spl

      END MODULE pies_fieldlines

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
module mod_convert_ub_to_bv
  use mod_prec
  use mod_error_id
  use mod_rotation
  use mod_orth_band_types
  use mod_band_types
  implicit none

  private

  public :: convert_ub_to_bv, d_convert_ub_to_bv, z_convert_ub_to_bv, &
       f_convert_ub_to_bv, f_d_convert_ub_to_bv, f_z_convert_ub_to_bv, &
       d_bv_of_ub, z_bv_of_ub, bv_of

  interface convert_ub_to_bv
     module procedure d_convert_ub_to_bv, z_convert_ub_to_bv
  end interface convert_ub_to_bv

  interface f_convert_ub_to_bv
     module procedure f_d_convert_ub_to_bv, f_z_convert_ub_to_bv
  end interface f_convert_ub_to_bv

  interface bv_of
     module procedure d_bv_of_ub, z_bv_of_ub
  end interface bv_of

contains

  function d_bv_of_ub(ub,error) result(bv)
    type(d_bv) :: bv
    type(d_ub), intent(in) :: ub
    type(error_info), intent(inout), optional :: error

    type(d_ub), allocatable :: ub1
    type(routine_info), parameter :: info=info_d_bv_of_ub

    if (failure(error)) return
    call push_id(info, error)

    bv=d_new_bv(get_n(ub), ub%lbw, ub%ubw)
    ub1=d_new_ub(get_n(ub),ub%lbw,  ub%ubw+1)
    call copy(ub1,ub)
    call d_convert_ub_to_bv(ub1,bv,error)

    call pop_id(error)
  end function d_bv_of_ub

  ! Errors:
  ! 0: no error
  ! 1: n<1
  ! 2: Insufficient storage in ub
  ! 3: Insufficient stroage in bv
  ! 4: ub%n /= bv%n
  subroutine d_convert_ub_to_bv(ub, bv, error)
    type(d_ub) :: ub
    type(d_bv) :: bv
    type(error_info), intent(inout), optional :: error
    type(routine_info), parameter :: info=info_d_convert_ub_to_bv

    if (failure(error)) return
    call push_id(info, error)

    if (get_n(ub) < 1) then
       call set_error(1, info, error); return
    end if
    if (get_ubwmax(ub) < ub%ubw+1 .and. ub%ubw < get_n(ub)-1) then
       print *, ub%ubw
       print *, get_ubwmax(ub)
       call set_error(2, info, error); return
    end if
    if (get_lbwmax(bv) < ub%lbw .or. get_ubwmax(bv) < ub%ubw) then
       call set_error(3, info, error); return
    end if
    if (get_n(ub) /= get_n(bv)) then
       call set_error(4, info, error); return
    end if
    call f_d_convert_ub_to_bv(ub%bc, get_n(ub), ub%lbw, ub%ubw, get_lbwmax(ub), &
         get_ubwmax(ub), ub%numrotsu, ub%jsu, ub%csu, ub%ssu, bv%br, bv%lbw, &
         bv%ubw, get_lbwmax(bv), &
         get_ubwmax(bv), bv%numrotsv, bv%ksv, bv%csv, bv%ssv)
    call pop_id(error)
  end subroutine d_convert_ub_to_bv

  subroutine f_d_convert_ub_to_bv(b_ub, n, lbw, ubw, lbwmax_ub, ubwmax_ub, numrotsu, &
       jsu, csu, ssu, b_bv, lbw_bv, ubw_bv, lbwmax_bv, ubwmax_bv, numrotsv, ksv, &
       csv, ssv)
    real(kind=dp), dimension(lbwmax_ub+ubwmax_ub+1,n), intent(inout) :: b_ub
    integer(kind=int32), intent(in) :: n, lbw, ubw, lbwmax_ub, ubwmax_ub, lbwmax_bv, ubwmax_bv
    integer(kind=int32), dimension(n), intent(in) :: numrotsu
    integer(kind=int32), dimension(ubwmax_ub,n), intent(in) :: jsu
    real(kind=dp), dimension(ubwmax_ub,n), intent(in) :: csu, ssu

    real(kind=dp), dimension(n,lbwmax_bv+ubwmax_bv+1), intent(out) :: b_bv
    integer(kind=int32), dimension(n), intent(out) :: numrotsv
    integer(kind=int32), dimension(n,ubwmax_bv), intent(out) :: ksv
    real(kind=dp), dimension(n,ubwmax_bv), intent(out) :: csv, ssv
    integer(kind=int32), intent(out) :: lbw_bv, ubw_bv

    integer(kind=int32) :: j, k, k0, k1, ubw1, lbw1
    type(d_rotation) :: rot
    logical :: full_ubw

    b_bv(:,1:lbw+ubw+1)=0.0_dp; numrotsv=0
    ssv(:,1:ubw)=0.0_dp; csv(:,1:ubw)=0.0_dp
    ksv(:,1:ubw)=0
    lbw_bv=lbw; ubw_bv=ubw
    if (n == 1) then
       b_bv(1,1)=b_ub(1,1);
       lbw_bv=0; ubw_bv=0; numrotsv=0; return
    end if
    ! must allow for temporary fill-in
    lbw1=lbw
    if (ubw < n-1) then
       ubw1=ubw+1
       b_ub=eoshift(b_ub,-1,dim=1_int64)
       full_ubw=.false.
    else
       ubw1=ubw
       full_ubw=.true.
    end if
    do k=1,n-2 ! size of trailing principal submatrix
       do j=1,numrotsu(n-k)
          call f_rotation_times_tbc(csu(j,n-k),ssu(j,n-k),b_ub,n,lbw1,ubw1,&
               n-k,0,jsu(j,n-k))
       end do
       k0=max(n-k+1,ubw+2)
       k1=min(n-k+ubw,n)
       numrotsv(n-(k+1))=max(k1-k0+1,0)
       do j=k0,k1
          rot=rgivens(get_el_bc(b_ub,ubw1,j-ubw1,j-1), get_el_bc(b_ub,ubw1,j-ubw1,j))
          ksv(n-(k+1),k1-j+1)=j-1
          csv(n-(k+1),k1-j+1)=rot%cosine; ssv(n-(k+1),k1-j+1)=rot%sine
          call tbc_times_rotation(b_ub,n,lbw1,ubw1,0,k+1,rot,j-1)
       end do
    end do
    ! Store the results in b_bv
    if (.not. full_ubw) then
       b_ub=eoshift(b_ub,1,dim=1_int64)
    end if
    call bc_to_br(b_ub, b_bv, lbw, ubw)
  end subroutine f_d_convert_ub_to_bv

  function z_bv_of_ub(ub,error) result(bv)
    type(z_bv) :: bv
    type(z_ub), intent(in) :: ub
    type(error_info), intent(inout), optional :: error

    type(z_ub), allocatable :: ub1
    type(routine_info), parameter :: info=info_z_bv_of_ub

    if (failure(error)) return
    call push_id(info, error)

    bv=z_new_bv(get_n(ub), ub%lbw, ub%ubw)
    ub1=z_new_ub(get_n(ub),ub%lbw,  ub%ubw+1)
    call copy(ub1,ub)
    call z_convert_ub_to_bv(ub1,bv,error)
    deallocate(ub1)
    call pop_id(error)
  end function z_bv_of_ub

  subroutine z_convert_ub_to_bv(ub, bv, error)
    type(z_ub) :: ub
    type(z_bv) :: bv
    type(error_info), intent(inout), optional :: error
    type(routine_info), parameter :: info=info_z_convert_ub_to_bv

    if (failure(error)) return
    call push_id(info, error)

    if (get_n(ub) < 1) then
       call set_error(1, info, error); return
    end if
    if (get_ubwmax(ub) < ub%ubw+1 .and. ub%ubw < get_n(ub)-1) then
       call set_error(2, info, error); return
    end if
    if (get_lbwmax(bv) < ub%lbw .or. get_ubwmax(bv) < ub%ubw) then
       call set_error(3, info, error); return
    end if
    if (get_n(ub) /= get_n(bv)) then
       call set_error(4, info, error); return
    end if
    call f_z_convert_ub_to_bv(ub%bc, get_n(ub), ub%lbw, ub%ubw, get_lbwmax(ub), &
         get_ubwmax(ub), ub%numrotsu, ub%jsu, ub%csu, ub%ssu, bv%br, bv%lbw, &
         bv%ubw, get_lbwmax(bv), &
         get_ubwmax(bv), bv%numrotsv, bv%ksv, bv%csv, bv%ssv)
    call pop_id(error)
  end subroutine z_convert_ub_to_bv

  subroutine f_z_convert_ub_to_bv(b_ub, n, lbw, ubw, lbwmax_ub, ubwmax_ub, numrotsu, &
       jsu, csu, ssu, b_bv, lbw_bv, ubw_bv, lbwmax_bv, ubwmax_bv, numrotsv, ksv, &
       csv, ssv)
    complex(kind=dp), dimension(lbwmax_ub+ubwmax_ub+1,n), intent(inout) :: b_ub
    integer(kind=int32), intent(in) :: n, lbw, ubw, lbwmax_ub, ubwmax_ub, &
         lbwmax_bv, ubwmax_bv
    integer(kind=int32), dimension(n), intent(in) :: numrotsu
    integer(kind=int32), dimension(ubwmax_ub,n), intent(in) :: jsu
    real(kind=dp), dimension(ubwmax_ub,n), intent(in) :: csu
    complex(kind=dp), dimension(ubwmax_ub,n), intent(in) :: ssu

    complex(kind=dp), dimension(n,lbwmax_bv+ubwmax_bv+1), intent(out) :: b_bv
    integer(kind=int32), dimension(n), intent(out) :: numrotsv
    integer(kind=int32), dimension(n,ubwmax_bv), intent(out) :: ksv
    real(kind=dp), dimension(n,ubwmax_bv), intent(out) :: csv
    complex(kind=dp), dimension(n,ubwmax_bv), intent(out) :: ssv
    integer(kind=int32), intent(out) :: lbw_bv, ubw_bv

    integer(kind=int32) :: j, k, k0, k1, ubw1, lbw1
    type(z_rotation) :: rot
    logical :: full_ubw

    b_bv(:,1:lbw+ubw+1)=(0.0_dp,0.0_dp); numrotsv=0
    ssv(:,1:ubw)=(0.0_dp,0.0_dp); csv(:,1:ubw)=0.0_dp
    ksv(:,1:ubw)=0
    lbw_bv=lbw; ubw_bv=ubw
    if (n == 1) then
       b_bv(1,1)=b_ub(1,1);
       lbw_bv=0; ubw_bv=0; numrotsv=0; return
    end if
    ! must allow for temporary fill-in
    lbw1=lbw
    if (ubw < n-1) then
       ubw1=ubw+1
       b_ub=eoshift(b_ub,-1,dim=1_int64)
       full_ubw=.false.
    else
       ubw1=ubw
       full_ubw=.true.
    end if
    do k=1,n-2 ! size of trailing principal submatrix
       do j=1,numrotsu(n-k)
          call f_rotation_times_tbc(csu(j,n-k),ssu(j,n-k),b_ub,n,lbw1,ubw1,&
               n-k,0,jsu(j,n-k))
       end do
       k0=max(n-k+1,ubw+2)
       k1=min(n-k+ubw,n)
       numrotsv(n-(k+1))=max(k1-k0+1,0)
       do j=k0,k1
          rot=rgivens(get_el_bc(b_ub,ubw1,j-ubw1,j-1), get_el_bc(b_ub,ubw1,j-ubw1,j))
          ksv(n-(k+1),k1-j+1)=j-1
          csv(n-(k+1),k1-j+1)=rot%cosine; ssv(n-(k+1),k1-j+1)=rot%sine
          call tbc_times_rotation(b_ub,n,lbw1,ubw1,0,k+1,rot,j-1)
       end do
    end do
    ! Store the results in b_bv
    if (.not. full_ubw) then
       b_ub=eoshift(b_ub,1,dim=1_int64)
    end if
    call bc_to_br(b_ub, b_bv, lbw, ubw)
  end subroutine f_z_convert_ub_to_bv

end module mod_convert_ub_to_bv

max_stars_repo_name: m-a-stewart/BandGivensWeight

max_issues_repo_name: m-a-stewart/BandGivensWeight

max_forks_repo_name: m-a-stewart/BandGivensWeight


text:
C-----------------------------------------------------------------------
      SUBROUTINE SPWGET(IROMB,MAXWV,EPS,EPSTOP,ENN1,ELONN1,EON,EONTOP)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:  SPWGET     GET WAVE-SPACE CONSTANTS
C   PRGMMR: IREDELL       ORG: W/NMC23       DATE: 96-02-29
C
C ABSTRACT: THIS SUBPROGRAM GETS WAVE-SPACE CONSTANTS.
C
C PROGRAM HISTORY LOG:
C   96-02-29  IREDELL
C
C USAGE:    CALL SPWGET(IROMB,MAXWV,EPS,EPSTOP,ENN1,ELONN1,EON,EONTOP)
C   INPUT ARGUMENTS:
C     IROMB    - INTEGER SPECTRAL DOMAIN SHAPE
C                (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL)
C     MAXWV    - INTEGER SPECTRAL TRUNCATION
C   OUTPUT ARGUMENTS:
C     EPS      - REAL ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
C     EPSTOP   - REAL (MAXWV+1)
C     ENN1     - REAL ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
C     ELONN1   - REAL ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
C     EON      - REAL ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
C     EONTOP   - REAL (MAXWV+1)
C
C SUBPROGRAMS CALLED:
C   SPEPS        COMPUTE UTILITY SPECTRAL FIELDS
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C
C$$$
      REAL EPS((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EPSTOP(MAXWV+1)
      REAL ENN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
      REAL ELONN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
      REAL EON((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EONTOP(MAXWV+1)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
      MXTOP=MAXWV+1
      CALL SPEPS(IROMB,MAXWV,EPS,EPSTOP,ENN1,ELONN1,EON,EONTOP)
      END

max_stars_repo_name: GMAO-SI-Team/jedi-stack

max_issues_repo_name: GMAO-SI-Team/jedi-stack