text:
      double complex function atrLRL(j1,j2,j3,j4,j5,j6,za,zb)
      implicit none
c---atrLRL is the amplitude for
c---q+(-p4)+Q-(-p2)+l+(-p5) ---> q-(p1)+Q+(p3)+l-(p6)
c---All outgoing particles are except p3 left-handed
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      include 'masses.f'
      integer j1,j2,j3,j4,j5,j6
      double complex atree
      double precision prop
      prop=s(5,6)/dsqrt((s(5,6)-wmass**2)**2+(wmass*wwidth)**2)
c---Note interchange of za,zb to effect complex conjugation
      atrLRL=atree('pp',j1,j2,j3,j4,j5,j6,zb,za)*prop
c      write(6,*) 'atrLRL',atrLRL
      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
      subroutine set_fraction( lurep    , notot   , syname, nomult, imultp,
     +                         procesdef, allitems, no_act, actlst, nbpr  )

!     Deltares Software Centre

!>/File
!>      set the substance fractions intothe processes structure

!     Created   : Aug   2012 by Jan van Beek

      use timers         !< performance timers
      use processet      !< use processet definitions
      implicit none

      ! arguments

      integer             , intent(in   ) :: lurep                  !< unit number report file
      integer             , intent(in   ) :: notot                  !< number of substances
      character(len=*)    , intent(in   ) :: syname(*)              !< substance names
      integer             , intent(in   ) :: nomult                 !< number of multiple substances
      integer             , intent(in   ) :: imultp(2,nomult)       !< multiple substance administration
      type(procespropcoll), intent(inout) :: procesdef              !< the proces definition
      type(itempropcoll)  , intent(inout) :: allitems               !< all items of the proces system
      integer             , intent(inout) :: no_act                 !< number of activated processes
      character(len=*)    , intent(inout) :: actlst(*)              !< list of activated processes
      integer             , intent(inout) :: nbpr                   !< number of processes

      ! local

      type(sfracsprop)                    :: sfracs          ! substance fraction properties
      integer(4)                          :: ithndl = 0      ! handle for performance timer
      if (timon) call timstrt( "set_fractions", ithndl )

      ! set the fractions structure

      call get_sfrac ( lurep , notot , syname, nomult, imultp,
     +                 sfracs)

      ! add atributes to processes ( and sfracs ?)

      call add_atrfrc( lurep , procesdef, sfracs)

      ! add calculation of the sum of the fractions

      call add_sumfrc( lurep , procesdef, allitems, sfracs, no_act,
     +                 actlst, nbpr     )

      ! add the fluxes to the fractions by expanding processes

      call expand_frc( lurep , procesdef, allitems, sfracs)

      ! add the fluxes to the fractions by copying processes

      call add_prcfrc( lurep , procesdef, allitems, sfracs, no_act,
     +                 actlst, nbpr     )

      ! add the fluxes to the fractions by adding a distribution process

      call add_flxfrc( lurep , procesdef, allitems, sfracs, no_act,
     +                 actlst, nbpr     )

      ! add the dispersion stochi

      call add_dspfrc( lurep , procesdef, sfracs)

      if (timon) call timstop( ithndl )
      return
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


!=====================================================|
module geometry_def
!=====================================================|
  integer, parameter :: nx = 400
  integer, parameter :: margin = 2
  integer, parameter :: nxg = nx + 2*margin
  integer, parameter :: mtype = 8

  ! geometry
  real(KIND(0.d0)), save :: xmax, xmin, dx
  real(KIND(0.d0)), dimension(nxg), save :: x  ! geometry
  real(KIND(0.d0)), dimension(nxg,mtype), save :: qq, qqm, qqp ! variable
  real(KIND(0.d0)), dimension(nxg), save :: cc ! characteristic velocity

  ! time variable
  integer, save :: ns, nd
  real(KIND(0.d0)), save :: time,dt

  ! data directory
  character (len = 10),save :: data_dir
end module geometry_def


text:
module test2007_176_mod1
     integer, protected :: val
   ! integer :: val
   ! integer, protected :: age
     integer :: age

   ! The bug here is that "protected" as a statement is not equivalent to its specification in the variable declaration.
     protected :: age

  contains
     subroutine set_val(arg)
          integer arg
          val = arg
     end subroutine
     subroutine set_age(arg)
          integer arg
          age  = arg
     end subroutine
end module

program dt_init01
     use test2007_176_mod1
     implicit none
     integer :: value, his_age
     call set_val(88)
     call set_age(38)

   ! This should fail (and does in ROSE) if the declaration used the "protected" attribute (but we don't have the semantic analysis in place).
     value = val

   ! This should fail if the "protected" statement is modified but not if the variable is read only.
     his_age = age

     print *, value, his_age
end program

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
SUBROUTINE solrem(refset)
!-----------------------------------------------------------------------
! remeshing routine for solid elements
!-----------------------------------------------------------------------
USE damp_db, ONLY:  updlon_damp
USE c_input, ONLY: openfi, ludat
USE npo_db,ONLY: coora, velnp, oldlb
USE meshmo_db,ONLY: itransf, numpn, coorn, v_new, strnd, p_new,   &
                    sg_old, sg_new, l_new, ne_new, elemp, nodset, headp,  &
                    tailp, delete_pline, &
                    nelbnd, cdbnd, lnbnd, nlbnd, fside, nodlb   !remeshing database
USE ctrl_db, ONLY: ntype,ndime,ndofn,nload,echo_chnode
USE kinc_db, ONLY : nvelr,naris
USE surf_db,ONLY: surfa, new_srf, dalloc_srf
USE param_db,ONLY: mlen
USE name_db,ONLY: output
USE SPR_Constants, ONLY: nvarg
USE cms_db, ONLY : nconm,updlon_cm
!USE outp_db, ONLY : updlon_outp
USE loa_db, ONLY : updlon_lo
!#if UNIX
!     USE ifport
!#else
   USE dflib
   USE dfport
!#endif
IMPLICIT NONE

  !--- Dummy variables
  CHARACTER(len=*),INTENT(IN):: refset   !label of the element set to modify
  !--- Local variables
  LOGICAL          :: found, axisym
  CHARACTER (len=mlen) :: command
  CHARACTER(len=mlen), PARAMETER :: dirgid ='C:\Utils\GiD\Compack1.0.5\'
  CHARACTER(len=mlen):: giddir
  CHARACTER(len=mlen):: gidv = 'SIMPACT'   !len=32
  INTEGER (kind=4) :: lengn = 27
  INTEGER (kind=4) :: istop=0, ng, nelm, nnod, ludat_keep
  INTEGER (kind=4), POINTER :: ifx_b(:,:)   !  boundary fixities
  INTEGER(kind=4),POINTER:: l_old(:,:)      ! (nnode,nelem) old elements
  REAL    (kind=8) :: fact

  REAL(kind=8):: scalef

  INTERFACE
    INCLUDE 'contac.h'
    INCLUDE 'elemnt.h'
    INCLUDE 'get_ifxb.h'
    INCLUDE 'plines.h'
    INCLUDE 'pregid.h'
    INCLUDE 'qloop.h'
    INCLUDE 'rmctrln.h'
    INCLUDE 'rnmesh.h'
    INCLUDE 'transf.h'
    INCLUDE 'transf1.h'
    INCLUDE 'transf_g.h'
  END INTERFACE

  found = .FALSE.
  axisym = ntype == 3

  !Extract a set of nodes used in the definition of refset
  CALL elemnt('NODSET',name=refset,istop=istop, flag2=found)  !==>nodset(numpo)
  !  Get conectivities of the set
  CALL elemnt('SLNODS',name=refset,lnod=l_old,flag2=found) ! ==>  l_old(nnod,nelm)
  nelm = SIZE(l_old,2)   !Get number of elements in the set
  nnod = SIZE(l_old,1)   !Get number of nodes by element
  ! get external surface/line definition
  CALL new_srf(surfa)         !Get memory for surface
  CALL elemnt('BOUNDA', name=refset, flag2=found)  !==>surfa
  !--- Generation of a list of nodes in the boundary
  CALL rmctrln(surfa%nelem,surfa%head,nelbnd,cdbnd,lnbnd,nlbnd,fside)
                                    !==> lnbnd(nelbnd)       Conectivities of segments
                                    !==> nlbnd(nline)        Number of lines in boundary
                                    !==> cdbnd(ndime,nelbnd) Coordinates of nodes
  CALL dalloc_srf(surfa)   !Release 'surfa' from memory
  ! Eliminates loops
  CALL qloop(nelbnd,cdbnd,nlbnd)
  !get boundary conditions on the boundary line ==> ifx_d
  CALL get_ifxb(nelbnd,lnbnd,ifx_b)
  ! analyze contour fixity codes and vertices to find polylines
  CALL plines(cdbnd,nlbnd,fside,ifx_b)
  ! prepare GiD batch file
  fact = scalef     !keep scale factor
  scalef = 1d4      !modify factor previous to writing background mesh
  CALL pregid(scalef,nnod,nelbnd,cdbnd,nlbnd,ifx_b)  !write background mesh

  CALL delete_pline(headp,tailp)  !Release memory of the boundary arrays
  ! generation of a new mesh
  lengn = GETENVQQ(gidv, giddir)     !get enviromental variable ==> giddir
  IF( lengn <= 0 )THEN
    giddir = dirgid
    lengn = LEN_TRIM (giddir)
  END IF
!#if UNIX
!  command='".'//giddir(1:lengn)//'gid.exe" -n -b '//TRIM(output)//'.bgi'
!  istop = SYSTEM(TRIM(command))
!#else
  command='"'//giddir(1:lengn)//'gid.exe" -n2 -b '//TRIM(output)//'.bgi' ! -n... not work in GiD 9 and later for Windows
  istop = SYSTEM(TRIM(command))
!#endif

!#if UNIX
!  command = 'rm -rf .//'//TRIM(output)//'.gid'
!  istop = SYSTEM(TRIM(command))   ! delete GiD files (& directory)
!#else
  command = 'rmdir /S /Q .\'//TRIM(output)//'.gid'
  istop = SYSTEM(TRIM(command))   ! delete GiD files (& directory)
!#endif

  CALL openfi(54)         !open mesh file to read
  ludat_keep = ludat      !remember standard value of LUDAT
  ludat = 54              !changed temporarily to allow the use of listen with this unit

  !read new coordinates and connectivities
  CALL rnmesh(scalef,f_axsym=axisym)
  scalef = fact !restore value once coordinates have been read

  ! transfer all the nodal variables from the old mesh to the new one
  IF (nnod ==3)ng=1          !number of Gauss points (1 for triangles)
  IF (nnod ==4)ng=4          !number of Gauss points (4 for quads)
  IF (itransf == 2) THEN      ! transfer from nodal smoothed vars
    CALL transf(numpn, coorn, coora, l_old, velnp, v_new ) ! transfer nodal vars
    !Smoothing of Gaussian variables in the remeshed set of elements
    CALL elemnt('SMOOTH', name=refset, istop=istop, flag2=found) !==>strnd(nstre,npoin)
    nvarg = SIZE(strnd,1)     !number of Gauss variables to transfer
    ALLOCATE(sg_new(nvarg,ng*ne_new))  !get memory for new Gauss points
    CALL transf1(ne_new,nnod,ng,coorn,coora,l_old,l_new,sg_new,strnd) ! transfer gauss vars
  ELSE IF (itransf == 1) THEN  !transfer all GP variables using
                               !superconvergent patch recovery (SCPR) algorithm
    CALL transf(numpn, coorn, coora, l_old, velnp, v_new)     ! transfer nodal vars
    CALL elemnt('INTVAR', name=refset, istop=istop, flag2=found) !==>sg_old(nvarg,nelem)
    nvarg = SIZE(sg_old,1)     !number of Gauss variables to transfer
    ALLOCATE(sg_new(nvarg,ng*ne_new))  !get memory for new Gauss points
    CALL transf_g(numpn, coorn,  coora, l_old, sg_old, l_new, sg_new)  ! transfer GP vars
  END IF
  ! exchange nodes - delete nodes of the old mesh and add new nodes ==> modify nodal database
  CALL exnods(refset,.TRUE.)   !Updated lagrangian
  ! reallocates the element set after remeshing, compute Gaussian variables
  ! exchange mesh - deallocate old mesh and put new mesh,
  CALL elemnt('REALLO', name=refset, flag2=found)
  ! update boundary conditions, delete those associated with old mesh & add new
  CALL newfix( )
  CLOSE(54,STATUS='DELETE')      !close and delete new mesh data file
  ludat = ludat_keep    !restore read file to original

  !--- Release memory

  DEALLOCATE(coorn, l_old, ifx_b)
  DEALLOCATE(lnbnd,cdbnd,nlbnd,fside)
  IF (ASSOCIATED(nodlb))  DEALLOCATE(nodlb)
  IF (ASSOCIATED(l_new))  DEALLOCATE(l_new)
  IF (ASSOCIATED(nodset)) DEALLOCATE(nodset)
  IF (ASSOCIATED(elemp))  DEALLOCATE(elemp)
  IF (ASSOCIATED(strnd))  DEALLOCATE(strnd)
  IF (ASSOCIATED(sg_new)) DEALLOCATE(sg_new)
  IF (ASSOCIATED(sg_old)) DEALLOCATE(sg_old)
  IF (ASSOCIATED(v_new))  DEALLOCATE(v_new)
  IF (ASSOCIATED(p_new))  DEALLOCATE(p_new)

  ! updates data bases
  echo_chnode = .FALSE.
  CALL elemnt ('UPDLON')   ! elements/segments connectivities
  CALL updlon_kc (naris,nvelr,ndime)   ! KINEMATIC conditions
  CALL updlon_cm (ndofn)               ! CONCENTRATED MASSES
  !CALL updlon_outp (oldlb)            ! output nodes (unnecesary)
  CALL updlon_damp ()                  ! DAMPING nodes and sets
  CALL updlon_lo (nload,oldlb)         ! LOADS
  CALL contac (task='UPDLON',iwrit=0)             ! CONTACT nodes (CONT 2 & 3 & 4 only)
  echo_chnode = .TRUE.

  RETURN
END SUBROUTINE solrem

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


text:
      SUBROUTINE DT_AAEVT(Nevts,Epn,Npmass,Npchar,Ntmass,Ntchar,Idp,
     &                    Iglau)

C***********************************************************************
C This version dated 22.03.96 is written by S. Roesler.                *
C***********************************************************************

      IMPLICIT NONE
      DOUBLE PRECISION dum , Epn
      INTEGER i , idmnyr , Idp , ievt , Iglau , irej , kkmat , Nevts ,
     &        nmsg , Npchar , Npmass , Ntchar , Ntmass
      SAVE

      INCLUDE 'inc/dtflka'
C emulsion treatment
      INCLUDE 'inc/dtcomp'
C event flag
      INCLUDE 'inc/dtevno'

      INCLUDE 'inc/dtflg1'

C      COMMON /PRTANU/ IAPROJ,IATARG

      CHARACTER*8 date , hhmmss
      DIMENSION idmnyr(3)

      kkmat = 1
      nmsg = MAX(Nevts/100,1)

C initialization of run-statistics and histograms
      CALL DT_STATIS(1)

      CALL PHO_PHIST(1000,dum)

C initialization of Glauber-formalism
      IF ( NCOmpo.LE.0 ) THEN
         CALL DT_SHMAKI(Npmass,Npchar,Ntmass,Ntchar,Idp,Epn,Iglau)
      ELSE
         DO i = 1 , NCOmpo
            CALL DT_SHMAKI(Npmass,Npchar,IEMuma(i),IEMuch(i),Idp,Epn,0)
         END DO
      END IF
      CALL DT_SIGEMU

      CALL IDATE(idmnyr)
      WRITE (date,'(I2,''/'',I2,''/'',I2)') idmnyr(1) , idmnyr(2) ,
     &       MOD(idmnyr(3),100)
      CALL ITIME(idmnyr)
      WRITE (hhmmss,'(I2,'':'',I2,'':'',I2)') idmnyr(1) , idmnyr(2) ,
     &       idmnyr(3)

      IF ( LPRi.GT.4 ) WRITE (LOUt,99010) date , hhmmss
99010 FORMAT (/,' DT_AAEVT: Initialisation finished. ( Date: ',A8,
     &        '   Time: ',A8,' )')

C generate NEVTS events
      DO ievt = 1 , Nevts

C  print run-status message
         IF ( MOD(ievt,nmsg).EQ.0 ) THEN
            CALL IDATE(idmnyr)
            WRITE (date,'(I2,''/'',I2,''/'',I2)') idmnyr(1) , idmnyr(2)
     &             , MOD(idmnyr(3),100)
            CALL ITIME(idmnyr)
            WRITE (hhmmss,'(I2,'':'',I2,'':'',I2)') idmnyr(1) ,
     &             idmnyr(2) , idmnyr(3)

            IF ( LPRi.GT.4 ) WRITE (LOUt,99020) ievt - 1 , Nevts ,
     &           date , hhmmss
99020       FORMAT (/,1X,I8,' out of ',I8,' events sampled ( Date: ',A,
     &              '   Time: ',A,' )',/)
C           WRITE(LOUT,1000) IEVT-1
C1000       FORMAT(1X,I8,' events sampled')
         END IF
         NEVent = ievt
C  treat nuclear emulsions
C  composite targets only
         IF ( IEMul.GT.0 ) CALL DT_GETEMU(Ntmass,Ntchar,kkmat,0)
         kkmat = -kkmat
C  sample this event
         CALL DT_KKINC(Npmass,Npchar,Ntmass,Ntchar,Idp,Epn,kkmat,irej)

C       chain fusion
C       WRITE(6,*)' IFUSION,IAPROJ,IATARG ',IFUSION,IAPROJ,IATARG
         IF ( (IFUsion.EQ.1) .AND. (Npmass.GT.12) .AND. (Ntmass.GT.12) )
     &        CALL DT_DENSITY
C
C
         CALL PHO_PHIST(2000,dum)

      END DO

C print run-statistics and histograms to output-unit 6

      CALL PHO_PHIST(3000,dum)

      CALL DT_STATIS(2)

      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
! RUN: %python %S/test_errors.py %s %flang_fc1

! Test use of implicitly declared variable in specification expression

subroutine s1()
  m = 1
contains
  subroutine s1a()
    implicit none
    !ERROR: No explicit type declared for 'n'
    real :: a(m, n)
  end
  subroutine s1b()
    !ERROR: Implicitly typed local entity 'n' not allowed in specification expression
    real :: a(m, n)
  end
end

subroutine s2()
  type :: t(m, n)
    integer, len :: m
    integer, len :: n
  end type
  n = 1
contains
  subroutine s2a()
    !ERROR: Implicitly typed local entity 'm' not allowed in specification expression
    type(t(m, n)) :: a
  end
  subroutine s2b()
    implicit none
    !ERROR: No explicit type declared for 'm'
    character(m) :: a
  end
end

subroutine s3()
  m = 1
contains
  subroutine s3a()
    implicit none
    real :: a(m, n)
    !ERROR: No explicit type declared for 'n'
    common n
  end
  subroutine s3b()
    ! n is okay here because it is in a common block
    real :: a(m, n)
    common n
  end
end

subroutine s4()
  implicit none
contains
  subroutine s4a()
    !ERROR: No explicit type declared for 'n'
    real :: a(n)
  end
end

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm


text:
module singleton_module

    implicit none
    private

    public :: single, get_instance, dispose_instance

    logical :: lock = .false.

    type single_t
        private
        integer :: value
    end type single_t

    type(single_t) :: single

contains

    function get_instance(value) result(single)
        integer, intent(in) :: value
        type(single_t) :: single
        if (lock) then
            print *, "Single instance already created."
            return
        else
            print *, "Creating single instance now."
            single%value = value
            lock = .true.
        end if
    end function get_instance

    subroutine dispose_instance(single)
        type(single_t), intent(inout) :: single
        print *, "Disposing single instance now."
        single%value = 0
        lock = .false.
    end subroutine dispose_instance

end module singleton_module

max_stars_repo_name: zoziha/Fortran-Design-Patterns

max_issues_repo_name: zoziha/Fortran-Design-Patterns

max_forks_repo_name: zoziha/Fortran-Design-Patterns


text:
program catalog_no21

  integer,parameter :: n=139, KB=8, IB=8

  real*8 a(n, n), b(n, n), c(n, n)

  do kk=1,n,KB
  do ii=1,n,IB
  do j=1,n,4
     do k=kk,min(n,kk+KB-1),4
        do i=ii,min(n,ii+IB-1)
           a(i,j) = a(i,j) + b(i,k)*c(k,j) &
                          + b(i,k+1)*c(k+1,j) &
                          + b(i,k+2)*c(k+2,j) &
                          + b(i,k+3)*c(k+3,j)
           a(i,j+1) = a(i,j+1) + b(i,k)*c(k,j+1) &
                          + b(i,k+1)*c(k+1,j+1) &
                          + b(i,k+2)*c(k+2,j+1) &
                          + b(i,k+3)*c(k+3,j+1)
           a(i,j+2) = a(i,j+2) + b(i,k)*c(k,j+2) &
                          + b(i,k+1)*c(k+1,j+2) &
                          + b(i,k+2)*c(k+2,j+2) &
                          + b(i,k+3)*c(k+3,j+2)
           a(i,j+3) = a(i,j+3) + b(i,k)*c(k,j+3) &
                          + b(i,k+1)*c(k+1,j+3) &
                          + b(i,k+2)*c(k+2,j+3) &
                          + b(i,k+3)*c(k+3,j+3)
         end do
      end do
   end do
   end do
   end do

end program catalog_no21

max_stars_repo_name: xevolver/xevxml

max_issues_repo_name: xevolver/xevxml

max_forks_repo_name: xevolver/xevxml


text:
Many restaurants have lunch specials to keep business flowing during lunch hours.  They are often helpful to students who do not have much money and limited time to cook or pack a lunch.  This is a page designed for people to list the lunch specials around town. (See also Cheap Food and lunch buffets)

 California Eclectic

  Farmers Kitchen Cafe Bioregional gourmet lunches deeply discounted to order online for prepaid subscribers. $1 off on any lunch special with UDash ticket. Month of September, UDash riders get subscriber price for lunch..go to www.http//farmerskitchencafe.com, select subscription meals and place your online order. (Beautiful gourmet lunches for as little as $4.95$6.95, depending on daily menu). In place of subscriber #, type udash. Also, see regular cafe menu online at the same site.
  Seasons  $ 6.95  11.95.  This high class restaurant offers great lunch prices for an upscale restaurant. Went there May 1st 2012 and no reduced lunch specials, just the special of the day.

 Chinese

  Davis Noodle city  $4.957.50 (lunch specials served until 5pm!) Many of their dishes come in cheaper lunch portions.  Their eggplant is very good.
  Ding How  $ 3.506.50 (3.50 take out special served with rice only).  the regular lunch specials 4.956.50 are served with fried rice, chow mein, soup, egg roll, and fruit.  They are known for their fake meats, ask for the vegetarian lunch special menu.
  House of Chang  $5.50$6.25 Soup, eggroll, fried wonton, steamed or fried rice, orange, meat entree.
  Hunan  $4.256.25,  Comes with either soup or salad (your choice) and a fried won ton.  However many people seem to have differing opinions on taste, quality, and shady business practices.
  Jade Garden  $ 4.505.95 Soup, eggroll, and your choice of chow mein, fried rice or white rice.
  Silver Dragon  $5.50  $6.75, served with hot and sour soup, egg roll, chow mein, and pork fried rice
  Wok of Flame  $5.25  5.95, served with pork fried rice and pork chow mein.
  Davis Lunchbox  $8.00, traditional Szechuan style, served with a main and a side dish, and a bottle of water.

 Italian

  Osteria Fasulo $ 15.00  3 courses Pasta Special (cup of soup or salad,choices of pasta.dessert) Monday though Friday 11:30am to 1:30pm
    ( reservations suggested 530 758 1324)

 Japanese

  Nobu Hiro  $5.507.99 various japanese lunch specials, including a bento box that has the option of double sashimi (any combination of tuna and salmon).  They also have a sushi buffet for $11.99+ tax

 Mediterranean

  Sams Mediterranean  Daily lunch special in which the item of the day comes with a free drink.

 Persian

  Downloaddish  Up to $8 full lunch portion with a free drink and $1 delivery service.

 Mexican

  El Mariachi  Daily lunch special, includes drink

 Thai

  Sophias  $7.95, lunch served with salad, spring roll and rice
  Thai Nakorn  $7.45  8.95, veggie and meat lunch specials, most around $7.95  (comes with rice, soup, wontons, and orange)
  Thai Recipes  $7.00  8.50 served with rice and soup.  Their lunch specials do not discount much from their normal prices (not much more than a dollar off per entree) but the food is quite good.
  KetMoRee  $11.00 non noodle dishes are served with rice. Same portions as dinner but with discounted price.
  Taste of Thai

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
SUBROUTINE pbj(bit,tsolve)

!-------------------------------------------------------------
!
!  Parallel Block Jacobi
!
!  Use the scalar flux solution to check for convergence
!  If not converged, move to next iteration by passing angular
!   flux moments to neighbors. Neighbors store in psi* vectors
!  If not converged and out of iterations, store everything
!   and exit this group
!  If converged, store everything and move to next group
!
!-------------------------------------------------------------

USE invar
USE solvar
IMPLICIT NONE
INTEGER, INTENT(IN) :: bit
INTEGER :: myztag, myytag, myxtag, znxt, zprv, ynxt, yprv, xnxt, xprv
INTEGER :: znxtag, zprtag, ynxtag, yprtag, xnxtag, xprtag, neqz, neqy, neqx, ierr
INTEGER, DIMENSION(3) :: istat
REAL*8, INTENT(OUT) :: tsolve
REAL*8 :: df, dfmx

INCLUDE 'mpif.h'

! Set up some tags
myztag = 100 + zrank
myytag = 100 + yrank
myxtag = 100 + xrank
znxt = zrank + 1
znxtag = 100 + znxt
zprv = zrank - 1
zprtag = 100 + zprv
ynxt = yrank + 1
ynxtag = 100 + ynxt
yprv = yrank - 1
yprtag = 100 + yprv
xnxt = xrank + 1
xnxtag = 100 + xnxt
xprv = xrank - 1
xprtag = 100 + xprv

! Each process checks its own fluxes, then perform MPI Reduction
IF (MINVAL(ABS(phiold)) >= tolr) THEN
   df = MAXVAL(ABS((phi - phiold)/phiold))
ELSE
   df = MAXVAL(ABS(phi - phiold))
END IF

IF (isize == 1) THEN
   IF (bc(1)/=1 .AND. bc(2)/=1 .AND. bc(3)/=1 .AND. bc(4)/=1 .AND. bc(5)/=1 .AND. bc(6)/=1) THEN
      df = 0.0
   END IF
END IF

! Get the overall max difference
CALL MPI_ALLREDUCE(df,dfmx,1,MPI_DOUBLE_PRECISION,MPI_MAX,allcomm,ierr)

IF (dfmx > err .AND. bit < bitmx) THEN
   ! Only print iteration information if requested
   IF (irank == root) THEN
      IF (itp == 1) WRITE(8,111) bit, dfmx
   END IF
   ! Set the previous iterate of the flux equal to the current
   phiold = phi
   ! Send psio to neighbors. Inner blocks have 24 communications (8 octants to
   ! 3 neigbors each). Then set the incoming data into psii

   ! Perform the sends 'up z'
   neqz = apo*xys
   ! Octant 1
   IF (zrank /= npz-1) CALL MPI_SEND(psioz(1,1,1),neqz,MPI_DOUBLE_PRECISION,znxt,znxtag,zcomm,ierr)
   IF (zrank /= 0) CALL MPI_RECV(psiiz(1,1,1),neqz,MPI_DOUBLE_PRECISION,zprv,myztag,zcomm,istat,ierr)
   ! Octant 2
   IF (zrank /= npz-1) CALL MPI_SEND(psioz(1,1,2),neqz,MPI_DOUBLE_PRECISION,znxt,znxtag,zcomm,ierr)
   IF (zrank /= 0) CALL MPI_RECV(psiiz(1,1,2),neqz,MPI_DOUBLE_PRECISION,zprv,myztag,zcomm,istat,ierr)
   ! Octant 3
   IF (zrank /= npz-1) CALL MPI_SEND(psioz(1,1,3),neqz,MPI_DOUBLE_PRECISION,znxt,znxtag,zcomm,ierr)
   IF (zrank /= 0) CALL MPI_RECV(psiiz(1,1,3),neqz,MPI_DOUBLE_PRECISION,zprv,myztag,zcomm,istat,ierr)
   ! Octant 4
   IF (zrank /= npz-1) CALL MPI_SEND(psioz(1,1,4),neqz,MPI_DOUBLE_PRECISION,znxt,znxtag,zcomm,ierr)
   IF (zrank /= 0) CALL MPI_RECV(psiiz(1,1,4),neqz,MPI_DOUBLE_PRECISION,zprv,myztag,zcomm,istat,ierr)

   ! Perform the sends 'down z'
   ! Octant 5
   IF (zrank /= 0) CALL MPI_SEND(psioz(1,1,5),neqz,MPI_DOUBLE_PRECISION,zprv,zprtag,zcomm,ierr)
   IF (zrank /= npz-1) CALL MPI_RECV(psiiz(1,1,5),neqz,MPI_DOUBLE_PRECISION,znxt,myztag,zcomm,istat,ierr)
   ! Octant 6
   IF (zrank /= 0) CALL MPI_SEND(psioz(1,1,6),neqz,MPI_DOUBLE_PRECISION,zprv,zprtag,zcomm,ierr)
   IF (zrank /= npz-1) CALL MPI_RECV(psiiz(1,1,6),neqz,MPI_DOUBLE_PRECISION,znxt,myztag,zcomm,istat,ierr)
   ! Octant 7
   IF (zrank /= 0) CALL MPI_SEND(psioz(1,1,7),neqz,MPI_DOUBLE_PRECISION,zprv,zprtag,zcomm,ierr)
   IF (zrank /= npz-1) CALL MPI_RECV(psiiz(1,1,7),neqz,MPI_DOUBLE_PRECISION,znxt,myztag,zcomm,istat,ierr)
   ! Octant 8
   IF (zrank /= 0) CALL MPI_SEND(psioz(1,1,8),neqz,MPI_DOUBLE_PRECISION,zprv,zprtag,zcomm,ierr)
   IF (zrank /= npz-1) CALL MPI_RECV(psiiz(1,1,8),neqz,MPI_DOUBLE_PRECISION,znxt,myztag,zcomm,istat,ierr)

   ! Update for reflective boundary conditions at zrank = 0 and npz-1. Only
   ! apply new boundary conditions to those on the z-face for each incoming octant
   IF (zrank == 0 .AND. bc(1) == 1) THEN         ! Reflective at low-z
      psiiz(:,:,1) = psioz(:,:,5)
      psiiz(:,:,2) = psioz(:,:,6)
      psiiz(:,:,3) = psioz(:,:,7)
      psiiz(:,:,4) = psioz(:,:,8)
   END IF
   IF (zrank == npz-1 .AND. bc(2) == 1) THEN     ! Reflective at high-z
      psiiz(:,:,5) = psioz(:,:,1)
      psiiz(:,:,6) = psioz(:,:,2)
      psiiz(:,:,7) = psioz(:,:,3)
      psiiz(:,:,8) = psioz(:,:,4)
   END IF

   ! Perform the sends 'up y'
   neqy = apo*xzs
   ! Quadrant 1
   IF (yrank /= npy-1) CALL MPI_SEND(psioy(1,1,1),neqy,MPI_DOUBLE_PRECISION,ynxt,ynxtag,ycomm,ierr)
   IF (yrank /= 0) CALL MPI_RECV(psiiy(1,1,1),neqy,MPI_DOUBLE_PRECISION,yprv,myytag,ycomm,istat,ierr)
   ! Quadrant 2
   IF (yrank /= npy-1) CALL MPI_SEND(psioy(1,1,2),neqy,MPI_DOUBLE_PRECISION,ynxt,ynxtag,ycomm,ierr)
   IF (yrank /= 0) CALL MPI_RECV(psiiy(1,1,2),neqy,MPI_DOUBLE_PRECISION,yprv,myytag,ycomm,istat,ierr)
   ! Quadrant 5
   IF (yrank /= npy-1) CALL MPI_SEND(psioy(1,1,5),neqy,MPI_DOUBLE_PRECISION,ynxt,ynxtag,ycomm,ierr)
   IF (yrank /= 0) CALL MPI_RECV(psiiy(1,1,5),neqy,MPI_DOUBLE_PRECISION,yprv,myytag,ycomm,istat,ierr)
   ! Quadrant 6
   IF (yrank /= npy-1) CALL MPI_SEND(psioy(1,1,6),neqy,MPI_DOUBLE_PRECISION,ynxt,ynxtag,ycomm,ierr)
   IF (yrank /= 0) CALL MPI_RECV(psiiy(1,1,6),neqy,MPI_DOUBLE_PRECISION,yprv,myytag,ycomm,istat,ierr)

   ! Perform the sends 'down y'
   ! Quadrant 3
   IF (yrank /= 0) CALL MPI_SEND(psioy(1,1,3),neqy,MPI_DOUBLE_PRECISION,yprv,yprtag,ycomm,ierr)
   IF (yrank /= npy-1) CALL MPI_RECV(psiiy(1,1,3),neqy,MPI_DOUBLE_PRECISION,ynxt,myytag,ycomm,istat,ierr)
   ! Quadrant 4
   IF (yrank /= 0) CALL MPI_SEND(psioy(1,1,4),neqy,MPI_DOUBLE_PRECISION,yprv,yprtag,ycomm,ierr)
   IF (yrank /= npy-1) CALL MPI_RECV(psiiy(1,1,4),neqy,MPI_DOUBLE_PRECISION,ynxt,myytag,ycomm,istat,ierr)
   ! Quadrant 7
   IF (yrank /= 0) CALL MPI_SEND(psioy(1,1,7),neqy,MPI_DOUBLE_PRECISION,yprv,yprtag,ycomm,ierr)
   IF (yrank /= npy-1) CALL MPI_RECV(psiiy(1,1,7),neqy,MPI_DOUBLE_PRECISION,ynxt,myytag,ycomm,istat,ierr)
   ! Quadrant 8
   IF (yrank /= 0) CALL MPI_SEND(psioy(1,1,8),neqy,MPI_DOUBLE_PRECISION,yprv,yprtag,ycomm,ierr)
   IF (yrank /= npy-1) CALL MPI_RECV(psiiy(1,1,8),neqy,MPI_DOUBLE_PRECISION,ynxt,myytag,ycomm,istat,ierr)

   ! Update for reflective boundary conditions at yrank = 0 and npy-1. Only
   ! apply new boundary conditions to those on the y-face for each incoming octant
   IF (yrank == 0 .AND. bc(3) == 1) THEN         ! Reflective at low-y
      psiiy(:,:,1) = psioy(:,:,4)
      psiiy(:,:,2) = psioy(:,:,3)
      psiiy(:,:,5) = psioy(:,:,8)
      psiiy(:,:,6) = psioy(:,:,7)
   END IF
   IF (yrank == npy-1 .AND. bc(4) == 1) THEN     ! Reflective at high-y
      psiiy(:,:,4) = psioy(:,:,1)
      psiiy(:,:,3) = psioy(:,:,2)
      psiiy(:,:,8) = psioy(:,:,5)
      psiiy(:,:,7) = psioy(:,:,6)
   END IF

   ! Perform the sends 'up x'
   neqx = apo*yzs
   ! Quadrant 1
   IF (xrank /= npx-1) CALL MPI_SEND(psiox(1,1,1),neqx,MPI_DOUBLE_PRECISION,xnxt,xnxtag,xcomm,ierr)
   IF (xrank /= 0) CALL MPI_RECV(psiix(1,1,1),neqx,MPI_DOUBLE_PRECISION,xprv,myxtag,xcomm,istat,ierr)
   ! Quadrant 4
   IF (xrank /= npx-1) CALL MPI_SEND(psiox(1,1,4),neqx,MPI_DOUBLE_PRECISION,xnxt,xnxtag,xcomm,ierr)
   IF (xrank /= 0) CALL MPI_RECV(psiix(1,1,4),neqx,MPI_DOUBLE_PRECISION,xprv,myxtag,xcomm,istat,ierr)
   ! Quadrant 5
   IF (xrank /= npx-1) CALL MPI_SEND(psiox(1,1,5),neqx,MPI_DOUBLE_PRECISION,xnxt,xnxtag,xcomm,ierr)
   IF (xrank /= 0) CALL MPI_RECV(psiix(1,1,5),neqx,MPI_DOUBLE_PRECISION,xprv,myxtag,xcomm,istat,ierr)
   ! Quadrant 8
   IF (xrank /= npx-1) CALL MPI_SEND(psiox(1,1,8),neqx,MPI_DOUBLE_PRECISION,xnxt,xnxtag,xcomm,ierr)
   IF (xrank /= 0) CALL MPI_RECV(psiix(1,1,8),neqx,MPI_DOUBLE_PRECISION,xprv,myxtag,xcomm,istat,ierr)

   ! Perform the sends 'down x'
   ! Quadrant 2
   IF (xrank /= 0) CALL MPI_SEND(psiox(1,1,2),neqx,MPI_DOUBLE_PRECISION,xprv,xprtag,xcomm,ierr)
   IF (xrank /= npx-1) CALL MPI_RECV(psiix(1,1,2),neqx,MPI_DOUBLE_PRECISION,xnxt,myxtag,xcomm,istat,ierr)
   ! Quadrant 3
   IF (xrank /= 0) CALL MPI_SEND(psiox(1,1,3),neqx,MPI_DOUBLE_PRECISION,xprv,xprtag,xcomm,ierr)
   IF (xrank /= npx-1) CALL MPI_RECV(psiix(1,1,3),neqx,MPI_DOUBLE_PRECISION,xnxt,myxtag,xcomm,istat,ierr)
   ! Quadrant 6
   IF (xrank /= 0) CALL MPI_SEND(psiox(1,1,6),neqx,MPI_DOUBLE_PRECISION,xprv,xprtag,xcomm,ierr)
   IF (xrank /= npx-1) CALL MPI_RECV(psiix(1,1,6),neqx,MPI_DOUBLE_PRECISION,xnxt,myxtag,xcomm,istat,ierr)
   ! Quadrant 7
   IF (xrank /= 0) CALL MPI_SEND(psiox(1,1,7),neqx,MPI_DOUBLE_PRECISION,xprv,xprtag,xcomm,ierr)
   IF (xrank /= npx-1) CALL MPI_RECV(psiix(1,1,7),neqx,MPI_DOUBLE_PRECISION,xnxt,myxtag,xcomm,istat,ierr)

   ! Update for reflective boundary conditions at xrank = 0 and npx-1. Only
   ! apply new boundary conditions to those on the x-face for each incoming octant
   IF (xrank == 0 .AND. bc(5) == 1) THEN         ! Reflective at low-x
      psiix(:,:,1) = psiox(:,:,2)
      psiix(:,:,4) = psiox(:,:,3)
      psiix(:,:,5) = psiox(:,:,6)
      psiix(:,:,8) = psiox(:,:,7)
   END IF
   IF (xrank == npx-1 .AND. bc(6) == 1) THEN     ! Reflective at high-x
      psiix(:,:,2) = psiox(:,:,1)
      psiix(:,:,3) = psiox(:,:,4)
      psiix(:,:,6) = psiox(:,:,5)
      psiix(:,:,7) = psiox(:,:,8)
   END IF

ELSE IF (dfmx < err) THEN
   bcnvf = 1
   tsolve = MPI_WTIME()
   IF (irank == root) THEN
      WRITE(8,'(2X,A,I5,A,ES11.3,/)') "System converged in ", bit, " iterations. Max error = ", dfmx
  END IF
ELSE IF (bit == bitmx) THEN
   bcnvf = 0
   tsolve = MPI_WTIME()
   IF (irank == root) THEN
      WRITE (8,'(/,2X,A,ES11.3)') "ERROR: system did not converge. Max error= ", dfmx
      WRITE (8,*) "  Will adversely affect other results. Aborting program."
   END IF
END IF

111 FORMAT(2X,'It ', I5,' Dfmx ', ES11.3)

RETURN
END SUBROUTINE pbj

max_stars_repo_name: lanl/pidots

max_issues_repo_name: lanl/pidots

max_forks_repo_name: lanl/pidots


text:
C    @(#)bldmdc.f	20.4 2/13/96
        subroutine bldmdc (kxx, kerr)
c
c       This subroutine generates entities into arrays "dcmtbs" and
c       "dcmtln".
c
        include 'ipfinc/parametr.inc'

        include 'ipfinc/blank.inc'
c	Global variables used:
c		mtdcln, mtdcbs
        include 'ipfinc/branch.inc'
c	Global variables used:
c		kx, ky, brnch
        include 'ipfinc/bus.inc'
c	Global variables used:
c		busdta, base
        include 'ipfinc/dcmt.inc'
c	Global variables used:
c		dcmtln(r*8), dcmtbs(r*8)
        include 'ipfinc/lfiles.inc'
c	Global variables used:
c		None
        include 'ipfinc/prt.inc'
c	Global variables used:
c		errbuf
c
c	Local variables
c

        common /is_batch / is_batch
c
        character dccode*1, xbuf*120
        integer p, char
        save
c
c**********************************************************************
c**********************************************************************
c       bus entry
c
        entry dcntbs (nb,kerr)
        kerr = 0
        mtdcbs = mtdcbs + 1
        if (mtdcbs .gt. MAXMDC) then
           write (errbuf(1),100) MAXMDC
  100      format (' More than ',i3,' type "BM" buses.')
           if (is_batch .eq. 0) then
             call prterx ('E',1)
           else
             call prterx ('F',1)
           endif
           kerr = 1
           mtdcbs = 1
        endif
        do 120 j = 1,36
  120   dcmtbs(j,mtdcbs) = 0.0d0
        call getchr(1,dccode,kbsdta(12,nb))
        call putchr_8(1,dccode,dcmtbs(4,mtdcbs))
        dcmtbs(1,mtdcbs) = nb
        dcmtbs(2,mtdcbs) = dble(base(nb))
        dcmtbs(3,mtdcbs) = kbsdta(9,nb)
        dcmtbs(5,mtdcbs) = dble(busdta(13,nb))
        dcmtbs(6,mtdcbs) = dble(busdta(14,nb))
        dcmtbs(7,mtdcbs) = dble(busdta(3,nb))
        dcmtbs(8,mtdcbs) = dble(busdta(7,nb))
        dcmtbs(9,mtdcbs) = dble(busdta(8,nb))
        dcmtbs(10,mtdcbs) = 0.01745329252d0*dble(busdta(5,nb))
        dcmtbs(11,mtdcbs) = 0.01745329252d0*dble(busdta(6,nb))
        dcmtbs(12,mtdcbs) = 0.01745329252d0*dble(busdta(10,nb))
        dcmtbs(29,mtdcbs) = 0.01745329252d0*dble(busdta(11,nb))
        dcmtbs(14,mtdcbs) = dble(busdta(4,nb))
c
c       check passive dc bus
c
        if (dccode .eq. ' ') then
           if (dcmtbs(3,mtdcbs) .ne. 0) then
              write (errbuf(1),150) dccode
  150         format (' Converter type "',a1,'" is reserved for ',
     1               'passive "BM" with blank commutator bus name.')
              errbuf(2)=' '
              call bcdbus (nb,xbuf)
              write (errbuf(3),170) xbuf(1:80)
  170         format (2x,'(',a80,')')
              call prterx ('W',3)
           endif
c
c       check non-passive dc bus
c
        else
c
c          Insert default angles if omitted
c
           if (dcmtbs(10,mtdcbs) .eq. 0.0d0) then
              dcmtbs(10,mtdcbs) = 0.01745329252d0 * 5.0d0
           endif
           if (dcmtbs(11,mtdcbs) .eq. 0.0d0) then
              dcmtbs(11,mtdcbs) = 0.01745329252d0 * 105.0d0
           endif
           if (dccode .eq. 'R') then
              if (dcmtbs(12,mtdcbs) .eq. 0.0d0) then
                 dcmtbs(12,mtdcbs) = 0.01745329252d0 * 15.0d0
              endif
           else if (dccode .eq. 'I') then
              if (dcmtbs(12,mtdcbs) .eq. 0.0d0) then
                 dcmtbs(12,mtdcbs) = 0.01745329252d0 * 21.0d0
              endif
           else if (dccode .eq. 'M') then
              if (dcmtbs(12,mtdcbs) .eq. 0.0d0) then
                 dcmtbs(12,mtdcbs) = 0.01745329252d0 * 16.0d0
              endif
           endif
           if (dcmtbs(29,mtdcbs) .eq. 0.0d0) then
              dcmtbs(29,mtdcbs) = 0.01745329252d0 * 16.0d0
           endif
           if (dcmtbs(10,mtdcbs) .ge. dcmtbs(11,mtdcbs)) then
              write (errbuf(1),190) dccode,busdta(5,nb),busdta(6,nb)
  190         format (' Type "BM" bus code "',a1,'" has amin (',f5.1,
     1                ') .ge. astop (',f5.1,')')
              errbuf(2)=' '
              call bcdbus (nb,xbuf)
              write (errbuf(3),170) xbuf(1:80)
              call prterx ('W',3)
           endif
           if (dccode .eq. 'I') then
              if (dcmtbs(29,mtdcbs) .ge. dcmtbs(12,mtdcbs)) then
                 write (errbuf(1),210) dccode,busdta(11,nb),
     &              busdta(10,nb)
  210            format (' Type "BM" bus code "',a1,'" has gmin (',
     1                   f5.1,') .ge. gnom (',f5.1,')')
                 errbuf(2)=' '
                 call bcdbus (nb,xbuf)
                 write (errbuf(3),170) xbuf(1:80)
                 call prterx ('W',3)
              endif
           else if (dccode .eq. 'M') then
              if (dcmtbs(29,mtdcbs) .gt. dcmtbs(12,mtdcbs)) then
                 write (errbuf(1),220) dccode,busdta(11,nb),
     &              busdta(10,nb)
  220            format ('0 type "BM" bus code "',a1,'" has gmin (',
     1                   f5.1,') .gt. gnom (',f5.1,')')
                 errbuf(2)=' '
                 call bcdbus (nb,xbuf)
                 write (errbuf(3),170) xbuf(1:80)
                 call prterx ('W',3)
              endif

           else if (dccode .eq. 'R') then

           else
              write (errbuf(1),230) dccode
  230         format (' type "BM" bus has illegal code "',a1,'".')
              errbuf(2)=' '
              call bcdbus (nb,xbuf)
              write (errbuf(3),170) xbuf(1:80)
              if (is_batch .eq. 0) then
                call prterx ('E',3)
              else
                call prterx ('F',3)
              endif
           endif

           if (dcmtbs(7,mtdcbs) .lt. 1.0d0) then
              write (errbuf(1),270) dcmtbs(7,mtdcbs)
  270         format (' Type "BM" bus has an illegal no. of bridges :',
     1                f4.1,'. Default is 1.0 bridge.')
              errbuf(2)=' '
              call bcdbus (nb,xbuf)
              write (errbuf(3),170) xbuf(1:80)
              call prterx ('W',3)
              dcmtbs(7,mtdcbs) = 1.0d0
           endif

        endif

        dcmtbs(8,mtdcbs) = 0.001d0*dcmtbs(7,mtdcbs)*dcmtbs(8,mtdcbs)

        return
c
c**********************************************************************
c**********************************************************************
c       line entry
c
        entry dcntln (p, kerr)

        k1  = kx(p)
        k2  = ky(p)
        nbr = brnch_ptr(p)
        jbr = iabs(nbr)
        kerr = 0
        do i = 1,mtdcln
           if (dcmtln(1,i) .eq. k1 .and. dcmtln(2,i) .eq. k2) goto 900
           if (dcmtln(1,i) .eq. k2 .and. dcmtln(2,i) .eq. k1) goto 900
        enddo

        mtdcln = mtdcln + 1
        if (mtdcln .gt. MAXDCL) then
           write (errbuf(1),320) MAXDCL
  320      format (' More than ',i3,' type "LM" branches.')
           mtdcln = 1
           kerr = 1
           if (is_batch .eq. 0) then
             call prterx ('E',1)
           else
             call prterx ('F',1)
           endif
        endif
        do 340 j = 7,10
  340   dcmtln(j,mtdcln) = 0
        dcmtln(1,mtdcln) = kx(p)
        dcmtln(2,mtdcln) = ky(p)
        do 350 j = 4,6
  350   dcmtln(j,mtdcln) = dble(brnch(j+1,jbr))
        if (dcmtln(4,mtdcln) .le. 0.0d0) then
           write (errbuf(1),360) dcmtln(4,mtdcln)
  360      format (' Type "LM" line has illegal resistance: ',f7.1,
     &             '.  Chaged to "0.001".')
           errbuf(2)=' '
           call bcdbrn (p,xbuf)
           write (errbuf(3),170) xbuf(1:80)
           call prterx ('W',3)
           dcmtln(4,mtdcln) = 0.001d0
        endif
  900   return
        end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


!=====================================
! projectile.f90
! Author: Devansh Shukla
!-------------------------------------

program projectile
    ! Program to compute motion of a particle moving in projectile motion.

    implicit none
    real*8, parameter :: PI=3.141592, g=9.81
    real*8 :: x, y, vx, vy
    real*8 :: theta, v0x, v0y, v0
    real*8 :: t0, tf, dt, t
    character(len=*), parameter :: fmt1 = "(F12.4, F12.4, F12.4, F12.4, F12.4)"

    open(unit=8, file="Projectile.dat")

    print *, "----------------------------------------"
    print *, "Enter v0, theta(deg)"
    read *, v0, theta
    print *, "Enter t0, tf, dt"
    read *, t0, tf, dt
    print *, "----------------------------------------"

    print *, "v0, theta =", v0, theta
    print *, "t0, tf, dt=", t0, tf, dt

    if (v0 .le. 0.0) stop "Illegal value of v0<=0"
    if (theta .lt. 0.0 .or. theta .gt. 90.0 ) &
    stop "Illegal value of theta"
    print *, "----------------------------------------"

    theta = (PI/180.0) * theta      ! theta radians
    v0x = v0*cos(theta)
    v0y = v0*sin(theta)

    print "(A12, A12, A12, A12, A12)", "time", "x(t)", "y(t)", "vx(t)", "vy(t)"
    t = t0
    do while(t <= tf)
        x = v0x * t
        y = v0y * t - 0.5*g*t*t
        vx = v0x
        vy = v0y - g*t
        write (*, fmt1) t, x, y, vx, vy
        write (8, fmt1) t, x, y, vx, vy
        if (y < 0.0) stop "y-ve"
        t = t + dt
    enddo
    print *, "----------------------------------------"
    close(8)
end program projectile


text:
submodule (parent_hier:grandchild_hier) greatgrandchild_hier

implicit none

contains

module procedure hi3

print *,"I'm a greatgrandchild, tau=",2*pi

end procedure hi3

end submodule greatgrandchild_hier

max_stars_repo_name: scivision/fortran-submodule

max_issues_repo_name: scivision/fortran-submodule

max_forks_repo_name: scivision/fortran-submodule


!!!-------------------------------------------!!!
!!!	Troy Joseph Raen
!!!	June 2017
!!!
!!!	add-on module for MESA
!!!	calculates extra heat transported by DM
!!!-------------------------------------------!!!
!!! controls:
!!!
!!! mxGeV = s% X_CTRL(6)	! DM mass in GeV
!!! sigmaxp = s% X_CTRL(7)	! DM-proton cross section, cm^2
!!!	cboost = s% X_CTRL(1)
!!!	spindep = s% X_LOGICAL_CTRL(1)
			! .true. = spin dependent; .false. = spin independent
!!!			! Spin independent option is not completely set up.
!!!			! Partial derivatives need to be calculated for
!!!			! use_other_energy_implicit.
!!!
!!!	extra history columns values = s% X_CTRL(2:5)
!!!
!!!
!!! Nx = s% xtra(1)
!!!			! in extras_finish_step (run_star_extras) s% xtra(1) = s% xtra(2)
!!! 		! so DM is not collected when step is not accepted
!!!
!!!-------------------------------------------!!!

	MODULE DM_module

	use star_def
	use const_def
	use chem_def
	use DM_num
	IMPLICIT NONE
	INCLUDE 'DM_vars.h'

	CONTAINS

!!----------------------------
!!	main routine, called by MESA run_star_extras.f
!!----------------------------
	SUBROUTINE DM_energy_transport(id,ierr)
	IMPLICIT NONE
	INTEGER, INTENT(IN) :: id
	INTEGER, INTENT(OUT) :: ierr
	INTEGER :: itr
	TYPE (star_info), pointer :: s ! pointer to star type
	ierr=0
	CALL GET_STAR_PTR(id, s, ierr)
	IF ( ierr /= 0 ) RETURN

	cboost = s% X_CTRL(1)  ! boost in capture rate of DM compared to the local capture rate near the Sun, \propto density/sigma_v
	spindep = s% X_LOGICAL_CTRL(1)

	CALL get_star_variables(id,ierr)
	CALL set_DM_variables(id,ierr)
	CALL calc_xheat(Tx)

	DO itr = 1,kmax
		s% extra_heat(itr) = xheat(itr)
		s% d_extra_heat_dlndm1(itr) = 0.D0
        s% d_extra_heat_dlnd00(itr) = d_xheat_dlnd00(itr)
        s% d_extra_heat_dlndp1(itr) = 0.D0
        s% d_extra_heat_dlnTm1(itr) = 0.D0
        s% d_extra_heat_dlnT00(itr) = d_xheat_dlnT00(itr)
        s% d_extra_heat_dlnTp1(itr) = 0.D0
        s% d_extra_heat_dlnR00(itr) = 0.D0
        s% d_extra_heat_dlnRp1(itr) = 0.D0
	ENDDO

	CALL store_hist(id,ierr)

	END SUBROUTINE DM_energy_transport

!!----------------------------
!!	grab star variables from MESA
!!	calculate vesc and V(k)
!!----------------------------
	SUBROUTINE get_star_variables(id,ierr)
	use const_def, only : mp, Rsun, standard_cgrav, amu
	! proton mass (g), solar radius (cm), Grav const (g^-1 cm^3 s^-2), amu (g)
	IMPLICIT NONE
	INTEGER, INTENT(IN) :: id
	INTEGER, INTENT(OUT) :: ierr
	INTEGER :: itr, j
	TYPE (star_info), pointer :: s ! pointer to star type
	ierr=0
	CALL GET_STAR_PTR(id, s, ierr)
	IF ( ierr /= 0 ) RETURN

	! star variables
	Age_star = s% star_age !! in years
	dttmp = s% dt !! in seconds
	maxT = 10.D0**(s% log_max_temperature) !! max temp in K
	M_star = s% mstar !! in grams
	R_star = (s% photosphere_r)* Rsun !! convert to cm
	vesc = SQRT(2.D0* standard_cgrav* M_star/ R_star)

	numspecies = s% species
	IF (numspecies .GT. maxspecies) THEN
		WRITE(*,*) '*** numspecies > maxspecies = ',maxspecies
		WRITE(*,*) '**** STOPPING RUN AT star_age = ',Age_star,' years'
		STOP
	ENDIF

	! copy cell variables
	kmax = s% nz
	IF ((kmax+1) .GT. maxcells) THEN
		WRITE(*,*) '**** kmax+1 > maxcells = ',maxcells
		WRITE(*,*) '**** STOPPING RUN AT star_age = ',Age_star,' years'
		STOP
	ENDIF

	DO itr = kmax,1, -1
		Xk(itr) = s% X(itr) !! mass fraction hydrogen
		Tk(itr) = s% T(itr) !! in K
		dmk(itr) = s% dm(itr)
		rhok(itr) = s% rho(itr) !! in g/cm^3
		npk(itr) = Xk(itr)*rhok(itr)/mp !! cm^-3
		rk(itr) = s% r(itr) !! in cm
		gravk(itr) = s% grav(itr) !! in cm/s^2
		! info on all elements in the net:
		IF (.NOT. spindep) THEN
			DO j = 1,numspecies
				chemj = s% chem_id(j) !! gives index of element j in chem_isos
				IF (itr .EQ. kmax) THEN
					mj(j) = chem_isos% W(chemj) * amu ! mass of element j (g)
					mGeVj(j) = mj(j)/gperGeV ! mass in GeV
					Aj(j) = chem_isos% Z_plus_N(chemj) ! mass number of element j
				ENDIF
				xajk(j,itr) = s% xa(j,itr) ! mass fraction of element j in cell k
				njk(j,itr) = xajk(j,itr)*rhok(itr)/mj(j) ! number fraction of element j in cell k
				IF (itr .EQ. kmax) THEN
!					WRITE(*,*) 'mod:  ',chem_isos% name(chemj), mGeVj(j), Aj(j), xajk(j,itr), njk(j,itr)
				ENDIF
			ENDDO
		ENDIF
	ENDDO

	! set central values (none set for elements beyond H1)
	Xk(kmax+1) = s% center_h1
	Tk(kmax+1) = 10.D0**(s% log_center_temperature)
	rhok(kmax+1) = 10.D0**(s% log_center_density)
	npk(kmax+1) = Xk(kmax+1)*rhok(kmax+1)/mp
	rk(kmax+1) = 0.D0
	gravk(kmax+1) = 0.D0
	Vk(kmax+1) = 0.D0

	! calculate V(k) from center to outer edge
	DO itr = kmax, 1, -1
		Vk(itr) = Vk(itr+1)+ 0.5D0*(gravk(itr)+ gravk(itr+1))* (rk(itr)- rk(itr+1)) ! cm^2/s^2
	ENDDO

	END SUBROUTINE get_star_variables

!!----------------------------
!!	set all DM properties
!!----------------------------
	SUBROUTINE set_DM_variables(id,ierr)
	IMPLICIT NONE
	DOUBLE PRECISION :: dNx
	INTEGER, INTENT(IN) :: id
	INTEGER, INTENT(OUT) :: ierr
	INTEGER :: itr, j
	TYPE (star_info), pointer :: s ! pointer to star type
	ierr=0
	CALL GET_STAR_PTR(id, s, ierr)
	IF ( ierr /= 0 ) RETURN

	mxGeV = s% X_CTRL(6)	! DM mass in GeV
	mx = mxGeV* gperGeV	! DM mass in grams

	sigmaxp = s% X_CTRL(7)	! DM-proton cross section, cm^2
	IF (.NOT.spindep) THEN
		DO j = 1,numspecies
			sigmaxj(j) = sigmaxp* (Aj(j)*mj(j)/mp)**2 * ((mx+mp)/(mx+mj(j)))**2
		ENDDO
	ENDIF

	Tx = calc_Tx(id,ierr)

!! in extras_finish_step (run_star_extras) s% xtra(1) = s% xtra(2)
!! so DM is not collected when step is not accepted
	dNx = calc_dNx()
	Nx = (s% xtra(1)) + dNx
	s% xtra(2) = Nx
	CALL calc_nxk()

	END SUBROUTINE set_DM_variables

!!----------------------------
!!	xheat is luminosity per unit mass (in cgs units) that DM gives up to baryons
!!	This is what needs to be given to MESA
!!	This is Spergel & Press 1985 equ 4.9 divided by rho(k)
!!----------------------------
	SUBROUTINE calc_xheat(Tx_given)
	USE const_def, only : pi, mp, kerg ! Boltzmann's constant (erg K^-1)
	IMPLICIT NONE
	INTEGER :: itr, j
	DOUBLE PRECISION :: mfact, dfact, Tfact, xheat_j, Tx_given

	IF (spindep) THEN
		mfact = 8.D0*SQRT(2.D0/pi)* sigmaxp* mx*mp/((mx+mp)**2) ! this is common to all cells

		DO itr = 1,kmax
			dfact = nxk(itr)*npk(itr)/rhok(itr)
			Tfact = SQRT((mp*kerg*Tx_given+ mx*kerg*Tk(itr))/(mx*mp)) * kerg*(Tx_given- Tk(itr))	! Tx_given-Tk gives correct sign
			xheat(itr) = mfact* dfact* Tfact

			! partial drvs for other_energy_implicit. all those not calculated here are zero.
			d_xheat_dlnd00(itr) = -xheat(itr)
			Tfact = mx*kerg/2.D0/(mp*kerg*Tx_given+ mx*kerg*Tk(itr)) - 1.D0/(Tx_given-Tk(itr))
			d_xheat_dlnT00(itr) = xheat(itr)*Tk(itr)*Tfact
		ENDDO

	ELSE
		DO itr = 1,kmax
			xheat(itr) = 0.D0
			d_xheat_dlnd00(itr) = 0.D0
			d_xheat_dlnT00(itr) = 0.D0
			DO j = 1,numspecies
				mfact = 8.D0*SQRT(2.D0/pi)* sigmaxj(j)* mx*mj(j)/((mx+mj(j))**2)
				dfact = nxk(itr)*njk(j,itr)/rhok(itr)
				Tfact = SQRT((mj(j)*kerg*Tx_given+ mx*kerg*Tk(itr))/(mx*mj(j))) * kerg*(Tx_given- Tk(itr))	! Tx_given-Tk gives correct sign
				xheat_j = mfact* dfact* Tfact
				xheat(itr) = xheat(itr)+ xheat_j
				! partial drvs for other_energy_implicit. all those not calculated here are zero.
				d_xheat_dlnd00(itr) = d_xheat_dlnd00(itr) - xheat_j
				Tfact = mx*kerg/2.D0/(mj(j)*kerg*Tx_given+ mx*kerg*Tk(itr)) - 1.D0/(Tx_given-Tk(itr))
				d_xheat_dlnT00(itr) = d_xheat_dlnT00(itr) + xheat_j*Tk(itr)*Tfact
			ENDDO
		ENDDO
	ENDIF

	END SUBROUTINE calc_xheat

!!----------------------------
!!	calculate DM number density for each cell
!!	assumes Nx has been calculated
!!----------------------------
	SUBROUTINE calc_nxk()
	USE const_def, only : pi, kerg ! Boltzmann's constant (erg K^-1)
	IMPLICIT NONE
	INTEGER :: itr
	DOUBLE PRECISION :: norm_integral, norm

	norm_integral = 0.D0
	DO itr=1,kmax!, 1, -1 ! integrate from r = 0 to Rstar
		norm_integral = norm_integral+ rk(itr+1)*rk(itr+1)* &
		EXP(-mx*Vk(itr)/ kerg/Tx)* (rk(itr)- rk(itr+1))
	ENDDO
	norm = Nx/ (4.D0*pi* norm_integral)

	nxk(kmax+1) = norm  ! this is central nx value since Vk(center) = 0
	DO itr = 1,kmax
		nxk(itr) = norm* EXP(-mx*Vk(itr)/ kerg/Tx)
	ENDDO

	END SUBROUTINE calc_nxk

!!----------------------------
!!	calculate number of captured DM particles in current time step
!!	crate is Zentner & Hearin 2011 equation 1
!!	with cboost subsuming rho and vbar factors
!!----------------------------
	FUNCTION calc_dNx()
	USE const_def, only : Msun ! solar mass (g)
	IMPLICIT NONE
	DOUBLE PRECISION :: calc_dNx, crate, Cfact

	IF (spindep) THEN
		Cfact = 5.D21* 5.D0/mxGeV ! s^-1
		crate = Cfact* cboost* sigmaxp/1.D-43* (vesc/6.18D7)**2* M_star/Msun
	ELSE
		Cfact = 7.D22 ! s^-1
		crate = Cfact* cboost* sigmaxp/1.D-43* (vesc/6.18D7)**2* M_star/Msun
	ENDIF

	calc_dNx = crate* dttmp

	END FUNCTION calc_dNx

!!----------------------------
!!	calculate DM temp using Spergel & Press 1985 one-zone model.
!!	this is the root of Spergel & Press 1985 equ 4.10
!!----------------------------
	FUNCTION calc_Tx(id,ierr)
	USE nrutil, ONLY : nrerror
	USE const_def, only : Lsun
	IMPLICIT NONE
	INTEGER, INTENT(IN) :: id
	INTEGER, INTENT(OUT) :: ierr
	INTEGER :: k
	DOUBLE PRECISION :: Txhigh, Txlow, tol
	DOUBLE PRECISION :: Ttmp, calc_Tx
	DOUBLE PRECISION :: Tarray(4)
	INTEGER :: tries, model_err=-1
	LOGICAL :: Tflag=.FALSE.
	PARAMETER ( tol = 1.D-4 )
	TYPE (star_info), pointer :: s ! pointer to star type
	ierr=0
	CALL GET_STAR_PTR(id, s, ierr)
	IF ( ierr /= 0 ) RETURN

	IF ((model_err.EQ. s% model_number) .AND. (.NOT.Tflag)) THEN
		call nrerror('Txlow > Txhigh or root must be bracketed')
	ENDIF

	Txhigh = maxT*1.1
	Txlow = maxT/25.0
	Tflag=.FALSE.
	tries=0

	DO WHILE ( .NOT. Tflag )
		tries = tries+1

		Tarray = zbrent(emoment, Txhigh, Txlow, tol) ! returns -1 if gets root must be bracketed error
		Ttmp = Tarray(1)
		IF (Txlow.GT.maxT) THEN ! treat as root must be bracketed error. Tx close to Tmax and slope is shallow (most likely)
			! WRITE(*,*) 'Txlow > maxT. Treating as root must be bracketed error. Model =', s% model_number
			Ttmp = -1.0
		ENDIF

		IF (Ttmp.GT.0.0) THEN
			Tflag = is_slope_negative(Ttmp)
			Txlow = 1.05*Txlow
		ELSE IF (tries.EQ.1) THEN ! expand the range and try again
			! WRITE(*,*) 'root must be bracketed. Txhigh=', Txhigh, 'Txlow=', Txlow, 'tries=', tries, 'Model =', s% model_number
			Txlow = Txlow/2.0
			Txhigh = Txhigh*1.25
		ELSE ! go back a step, find root, make sure slope is negative
			! WRITE(*,*) 'root must be bracketed. Txhigh=', Txhigh, 'Txlow=', Txlow, 'tries=', tries, 'Model =', s% model_number
			Txlow = Txlow/1.05
			Tarray = zbrent(emoment, Txhigh, Txlow, tol)
			Ttmp = Tarray(1)
			Tflag = is_slope_negative(Ttmp)
			IF (.NOT.Tflag) THEN
				! If the code gets here, a suitable root cannot be found.
				! In most cases there are other problems with the star at this
				! timestep and MESA will re-do the step.
				! The model_err variable (below) keeps track of this and will
				! kill the run if the step is not re-done, since this means
				! a suitable DM temperature was never found.
				Ttmp= 10.**(s% log_center_temperature)
				! WRITE(*,*) '**** Tx set to center_T ****', Ttmp, 'problem model ', s% model_number
				model_err= s% model_number +1
				EXIT ! terminate the DO WHILE loop
			ENDIF

		ENDIF

	ENDDO

	Tx_array = Tarray ! set DM_vars Tx_array
	calc_Tx = Ttmp
	END FUNCTION calc_Tx

!!----------------------------
!!	if slope < 0, return .TRUE.
!!	This helps to find the correct root of the energy equation
	FUNCTION is_slope_negative(Tx_given)
		IMPLICIT NONE
		DOUBLE PRECISION, INTENT(IN) :: Tx_given
		DOUBLE PRECISION :: Tl, Tx_int, Tl_int, slope
		LOGICAL :: is_slope_negative

		Tx_int = emoment(Tx_given)
		Tl = 0.999*Tx_given ! use narrower range to avoid function maximum
		Tl_int = emoment(Tl)
		slope = (Tx_int-Tl_int)/(Tx_given-Tl)

		IF (slope.LT.0.0) THEN
			is_slope_negative=.TRUE.
			! WRITE(*,*) 'is_slope_negative returns true. slope=',slope
		ELSE
			is_slope_negative=.FALSE.
			! WRITE(*,*) 'is_slope_negative returns false. slope=',slope
		ENDIF

	END FUNCTION is_slope_negative
!!----------------------------

!!----------------------------
!!	To avoid finding wrong (lower) root and inducing Tx oscillations
!!	calculate slope of Ttmp
!!	lower root has slope ~0, correct slope is usually very steep
!!	if abs(slope) < 10, return .FALSE.
!!
!!	This function turned out to be unnecessary and is not called from the code.
	FUNCTION is_slope_steep(Tx_given)
		IMPLICIT NONE
		DOUBLE PRECISION, INTENT(IN) :: Tx_given
		DOUBLE PRECISION :: Tl, Tx_int, Tl_int, slope
		LOGICAL :: is_slope_steep

		Tx_int = emoment(Tx_given)
		Tl = 0.99*Tx
		Tl_int = emoment(Tl)
		slope = ABS((Tx_int-Tl_int)/(Tx_given-Tl))
		! WRITE(*,*) 'slope=',slope

		IF (slope.GT.10.0) THEN
			is_slope_steep=.TRUE.
			! WRITE(*,*) 'is_slope_steep returns true. slope=',slope
		ELSE
			is_slope_steep=.FALSE.
			! WRITE(*,*) 'is_slope_steep returns false. slope=',slope
		ENDIF

	END FUNCTION is_slope_steep
!!----------------------------

!!----------------------------
!!	This is LHS of Spergel & Press 1985 equ 4.10 which implicitly defines Tx
!!	Called by zbrent() (which is called by calc_Tx())
!!	zbrent() finds Tx as the root of this equation
!!----------------------------
	FUNCTION emoment(Txtest)
	USE const_def, only : pi, Rsun, kerg ! Boltzmann's constant (erg K^-1)
	IMPLICIT NONE
	INTEGER :: itr, j
	DOUBLE PRECISION, INTENT(IN) :: Txtest
	DOUBLE PRECISION :: mpGeV, Tfact, efact, rfact, mfact, mjfact, sum, emoment
	DOUBLE PRECISION :: Tnorm, nnorm, Rnorm, m, npbar, Txbar, Tbar, rbar, drbar
	PARAMETER ( mpGeV=0.938272D0 ) ! Proton mass in GeV

	sum = 0.D0
	DO itr = kmax,1,-1 ! integrate over dm from r=0 to r_star
		mfact =  dmk(itr)/4./pi/rhok(itr) ! r^2 dr = dm / (4 pi rho)
		efact = EXP(-mx*Vk(itr)/ kerg/Txtest)
		IF (spindep) THEN
			Tfact = SQRT((mpGeV*Txtest+ mxGeV*Tk(itr))/(mxGeV*mpGeV))* (Tk(itr)-Txtest)
			sum = sum+ npk(itr)*Tfact*efact*mfact
		ENDIF
	ENDDO

	emoment = sum
	END FUNCTION emoment

!!----------------------------
!!	store data so it can be written to history.data
!!	THIS INFO MUST MATCH subroutine data_for_extra_history_columns
!!	IN run_star_extras !!!
!!----------------------------
	SUBROUTINE store_hist(id,ierr)
		IMPLICIT NONE
		INTEGER, INTENT(IN) :: id
		INTEGER, INTENT(OUT) :: ierr
		INTEGER :: j, k, idx
		TYPE (star_info), pointer :: s ! pointer to star type
		ierr=0
		CALL GET_STAR_PTR(id, s, ierr)
		IF ( ierr /= 0 ) RETURN

		! store info to write to extra_history_columns
		! indicies offset by one since cboost = s% X_CTRL(1)

		s% X_CTRL(2) = Tx ! names(1) = 'DM_temp'
		s% X_CTRL(3) = Nx ! names(2) = 'Nx_total'
		s% X_CTRL(4) = nxk((s% nz)+1) ! names(3) = 'center_nx'
		s% X_CTRL(5) = npk((s% nz)+1) ! names(4) = 'center_np'

		DO k = 1,kmax
			s% xtra1_array(k) = nxk(k)
			s% xtra2_array(k) = npk(k)
		ENDDO

	END SUBROUTINE store_hist

!!!***!!!*** TEST FUNCTIONS FOR THE Tx ROOT FINDING !!!***!!!***
! These functions are not currently used in the code, but can be used
! to understand the Tx root finding better.

!!----------------------------
!!	If extra energy is "too high", approximate emoment function with a
!!	straight line and return a root closer to actual zero
	FUNCTION linear_root(Tarray)
		DOUBLE PRECISION :: x1, x2, y1, y2
		DOUBLE PRECISION :: m, linear_root
		DOUBLE PRECISION :: Tarray(4)

		x1 = Tarray(1)
		y1 = Tarray(2)
		x2 = Tarray(3)
		y2 = Tarray(4)
		m = (y2-y1)/(x2-x1)
		linear_root = x1 - y1/m ! Tx

		! WRITE(*,*) "linear_root---***---"
		! WRITE(*,*) "linear_root---***--- slope = ", m
		! WRITE(*,*) "linear_root---***--- Tx OLD = ", x1
		! WRITE(*,*) "linear_root---***--- emoment OLD = ", y1
		! WRITE(*,*) "linear_root---***--- Tx NEW = ", linear_root
		! WRITE(*,*) "linear_root---***--- emoment NEW = ", emoment(linear_root)

	END FUNCTION linear_root

!!!!!!!!
! have Tx_array = [ Tx1, emom1, Tx2, emom2 ] from DM_num root finding.
! find linear_root Tx and emom with this Tx
! calc xenergy for all these Tx's
! write to separate file for plotting
!!!!!!!!
	SUBROUTINE energy_plots(id,ierr)
		IMPLICIT NONE
		INTEGER, INTENT(IN) :: id
		INTEGER, INTENT(OUT) :: ierr
		INTEGER :: i, k
		INTEGER, PARAMETER :: cols=853
		DOUBLE PRECISION :: matrx(3,cols)
		! matrx = 	[ Tx1, emom(Tx1), xEnergy(Tx1) xheat(Tx1) for each zone ]
		!			[ Tx2, emom(Tx2), xEnergy(Tx2) xheat(Tx2) for each zone ]
		!			[ Tx_linapprox, emom(Tx_linapprox), xEnergy(Tx_linapprox) xheat(Tx_linapprox) for each zone ]
		TYPE (star_info), pointer :: s ! pointer to star type
		ierr=0
		CALL GET_STAR_PTR(id, s, ierr)
		IF ( ierr /= 0 ) RETURN

		matrx(1,1) = Tx_array(1)
		matrx(1,2) = Tx_array(2)
		matrx(2,1) = Tx_array(3)
		matrx(2,2) = Tx_array(4)
		matrx(3,1) = linear_root(Tx_array)
		matrx(3,2) = emoment(matrx(3,1))

		DO i = 1,3
			CALL calc_xheat(matrx(i,1))
			matrx(i,3) = calc_xenrgy(id,ierr)
            DO k = 1,kmax
                matrx(i,k+3) = xheat(k)
            ENDDO
            DO k = kmax+1,850
                matrx(i,k+3) = 0.0D0
            ENDDO
		ENDDO

		OPEN(UNIT=9, FILE="/home/tjr63/mesaruns/LOGS/matrx.data", STATUS="NEW", ACTION="WRITE")
		DO i = 1,3
			WRITE(UNIT=9, FMT="(9999ES25.15)") matrx(i,:)
		ENDDO
		CLOSE(UNIT=9)

		STOP

	END SUBROUTINE energy_plots

!!! ONLY USED FOR energy_plots:
	FUNCTION calc_xenrgy(id, ierr)
		INTEGER, INTENT(IN) :: id
		INTEGER, INTENT(OUT) :: ierr
		real(dp) :: xe, calc_xenrgy
		integer :: k
		type (star_info), pointer :: s
		ierr = 0
		call GET_STAR_PTR(id, s, ierr)
		if (ierr /= 0) return

		xe = 0.d0
		DO k = 1, kmax
			xe = xe + xheat(k)* s% dm(k)* s% dt
		ENDDO

		calc_xenrgy = xe ! ergs
	END FUNCTION calc_xenrgy
!!! ONLY USED FOR energy_plots

!!!***!!!*** END TEST FUNCTIONS !!!***!!!***

END MODULE DM_module


text:
Living in Davis one is often exposed to English, Spanish and Greek.  To assist with deciphering Fraternity and Sororities Sorority names, or for reference in your math courses, use the following table:

Does this page contribute anything that isnt already found in the wiki:WikiPedia:Greek_alphabet Wikipedia article on the Greek alphabet?  ...especially, anything of specific relevance to Davis?
  On my first pass through I didnt add links from each letter to the organizations that refer to themselves as such, but that would help to make it Davis relevant.


text:
   program veradiff2
   use hdf5
   implicit none
!=======================================================================
!
!  Program to compare two VERA HDF output file and print summary
!
!  *** This version allows the comparison of two files with different numbers of exposure points ***
!
!  *** Use veradiff_testing for automated testing applications ****
!
!  Copyright (c) 2014-2017 Core Physics, Inc.
!
!  Distributed under the MIT license.
!  See the LICENSE file in the main directory for details.
!
!  Main directory: https://github.com/palmtag/VERAout-tools/
!
!=======================================================================
!
!  This program is intended to be used in regressing testing
!
!  Right now, we are just performing a simple diff on eigenvalues with a tolerance of 10 pcm.
!  This will be modifed in the future to check other parameters as well.
!
!  2013/09/06 - original file
!  2014/01/27 - compare boron if dataset is present on file
!  2014/04/26 - add comparison for multiple statepoints - see #3164
!  2014/11/01 - added pin power comparisons             - see #3430
!  2014/12/05 - added pin power tolerance               - see #3430
!  2016/02/28 - add statepoint summary
!  2017/05/01 - add ability to test files with different exposure points
!             - created new version "veradiff_testing" for automated testing applications
!
!--------------------------------------------------------------------------------

      character(len=250) :: fname1, fname2   ! HDF file names
      integer            :: iargs            ! number of command line arguments

      integer(hid_t)     :: file_id1, file_id2  ! HDF file id's
      integer            :: ierror
      integer            :: i, n, n2
      integer            :: nstate           ! statepoint number
      integer            :: nstate2          ! statepoint number

      integer, parameter :: maxstate=200

      real(8) :: xkeff1(maxstate)
      real(8) :: xkeff2(maxstate)
      real(8) :: xexpo1(maxstate)
      real(8) :: xexpo2(maxstate)
      real(8) :: xefpd1(maxstate)
      real(8) :: xefpd2(maxstate)
      real(8) :: xboron1(maxstate)
      real(8) :: xboron2(maxstate)
      real(8) :: dk, db
      real(8) :: xmax(maxstate)
      real(8) :: xrms(maxstate)

      logical            :: ifxst

      ierror=0
      fname1=' '
      fname2=' '

      nstate=0   ! number of statepoints on file 1
      nstate2=0  ! number of statepoints on file 2

      xexpo1=0.0d0
      xexpo2=0.0d0

!----------------------------------------------------------------------
!  Read in arguments from command line
!----------------------------------------------------------------------

      iargs = command_argument_count()
      if (iargs.lt.2) then
        write (*,*) 'usage:  veradiff.exe [file1] [file2]'
        goto 900
      endif

      call get_command_argument(1,fname1)
      call get_command_argument(2,fname2)

      write (*,'(/,1x,2a)') 'reading h5 file: ', trim(fname1)
      inquire(file=fname1, exist=ifxst)
      if (.not.ifxst) then
        write (*,*) 'ERROR: input file ',trim(fname1),' does not exist'
        goto 900
      endif

      write (*,'(1x,2a)') 'reading h5 file: ', trim(fname2)
      inquire(file=fname2, exist=ifxst)
      if (.not.ifxst) then
        write (*,*) 'ERROR: input file ',trim(fname2),' does not exist'
        goto 900
      endif

!--------------------------------------------------------------------------------
! Read HDF files
!--------------------------------------------------------------------------------

!--- initialize fortran interface

      call h5open_f(ierror)      ! NOTE: THIS IS REQUIRED

!--- open file1

      call h5fopen_f (fname1, H5F_ACC_RDONLY_F, file_id1, ierror)
      if (ierror<0) then
        write (*,'(3a)') 'ERROR: H5 input file ',trim(fname1),' could not be opened'
        goto 900
      endif

!--- open file2

      call h5fopen_f (fname2, H5F_ACC_RDONLY_F, file_id2, ierror)
      if (ierror<0) then
        write (*,'(3a)') 'ERROR: H5 input file ',trim(fname2),' could not be opened'
        goto 900
      endif

!--- read number of statepoints and eigenvalues

      call countstates(file_id1, maxstate, nstate, xexpo1, xefpd1, xkeff1, xboron1)

      call countstates(file_id2, maxstate, nstate2, xexpo2, xefpd2, xkeff2, xboron2)

!---- calculate pin difference values (experimental)

      do n=1, nstate
        n2=0
        do i=1, nstate2
          if (abs(xexpo1(n)-xexpo2(i)).lt.0.0001d0) n2=i
        enddo
        if (n2.eq.0) cycle  ! no match found

        call pin_compare(file_id1, file_id2, n, n2, xmax(n), xrms(n))

      enddo

!--------------------------------------------------------------------------------
! Print Results
!--------------------------------------------------------------------------------

      write (*,*)
      write (*,'(2a)') ' file1: ', trim(fname1)
      write (*,'(2a)') ' file2: ', trim(fname2)

      write (*,122)

      do n=1, nstate
        n2=0
        do i=1, nstate2
          if (abs(xexpo1(n)-xexpo2(i)).lt.0.0001d0) n2=i
        enddo
        if (n2.eq.0) cycle  ! no match found

        dk=(xkeff2(n2)-xkeff1(n))*1.0d5
        db= xboron2(n2)-xboron1(n)

        write (*,55) n, xexpo1(n), xefpd1(n), xkeff1(n), xkeff2(n2), dk, xboron1(n), xboron2(n2), db, xmax(n), xrms(n)

      enddo

  122 format ( &
            '  ==============================================================================', &
          /,'                         Statepoint Summary (file2-file1)', &
          /,'  ==============================================================================', &
          /,'     N     exp     EFPD    keff1     keff2      pcm     boron1   boron2     dbor   pinmax   pinrms')
   55 format (2x,i4, f9.4, f8.2, 2f10.6, f8.1, 2x, 3f9.2, 2f9.2)

!--- close HDF files

      call h5fclose_f(file_id1, ierror)
      call h5fclose_f(file_id2, ierror)

  900 continue

      end program

!=======================================================================
!
!  Subroutine to read two pin_power datasets and calculate RMS
!
!=======================================================================
      subroutine pin_compare(file_id1, file_id2, nstate1, nstate2, xmax, xrms)
      use hdf5
      use mod_hdftools, only : h5info, hdf5_read_double
      implicit none
      integer(hid_t),   intent(in)  :: file_id1, file_id2
      integer,          intent(in)  :: nstate1, nstate2
      real(8),          intent(out) :: xmax   ! return max difference (%)
      real(8),          intent(out) :: xrms   ! return RMS difference (%)

      integer :: ip, jp, k, na, np
      integer :: nassm, kd, npin
      integer :: ndim         ! number of array dimensions
      integer :: idim(10)     ! array dimension sizes
      integer :: itype        !
      integer :: kmax1(4)     ! location of max difference
      integer :: lmax1(4)     ! location of max power file 1
      integer :: lmax2(4)     ! location of max power file 2

      logical :: ifbad        ! flag for bad data

      real(8)              :: p1, p2, dd, dmax
      real(8)              :: pmax1, pave1
      real(8)              :: pmax2, pave2
      real(8)              :: rms
      real(8), allocatable :: power1(:,:,:,:)
      real(8), allocatable :: power2(:,:,:,:)

      character(len=24)  :: dataset1
      character(len=24)  :: dataset2

!--- dataset names

      dataset1='STATE_0001/pin_powers'
      dataset2='STATE_0001/pin_powers'

      write (dataset1(7:10),'(i4.4)') nstate1
      write (dataset2(7:10),'(i4.4)') nstate2

!d    write (*,*) 'debug: ', dataset1, nstate1
!d    write (*,*) 'debug: ', dataset2, nstate2

!--- find dimensions so we can allocate

      xmax=999999.0d0   ! set high in case reads fail
      ifbad=.false.

      call h5info(file_id1, dataset1, itype, ndim, idim)

      nassm=idim(1)    ! order (nassm, kd, npin, npin)
      kd   =idim(2)
      npin =idim(4)

!d    write (*,*) 'debug pin: ndim  = ', ndim
!d    write (*,*) 'debug pin: nassm = ', nassm
!d    write (*,*) 'debug pin: kd    = ', kd
!d    write (*,*) 'debug pin: npin  = ', npin

      if (ndim .ne.4) then
        ifbad=.true.
        write (*,*) 'invalid dimensions in pin data'
      endif
      if (nassm.eq.0) then
        ifbad=.true.
        write (*,*) 'invalid number of assemblies in pin data'
      endif
      if (npin .eq.0) then
        ifbad=.true.
        write (*,*) 'invalid number of pins in pin data'
      endif
      if (idim(3).ne.idim(4)) then
        ifbad=.true.
        write (*,*) 'invalid npin in pin data'
      endif
      if (ifbad) return   ! return with big xmax

      allocate (power1(nassm, kd, npin, npin))
      allocate (power2(nassm, kd, npin, npin))

      call hdf5_read_double(file_id1, dataset1, idim(1), idim(2), idim(3), idim(4), power1)
      call hdf5_read_double(file_id2, dataset2, idim(1), idim(2), idim(3), idim(4), power2)

      rms=0.0d0
      dmax=0.0d0
      np=0
      kmax1(:)=0    ! location of max difference
      lmax1(:)=0    ! location of max power file 1
      lmax2(:)=0    ! location of max power file 2
      pmax1=0.0d0   ! max power of distribution 1
      pmax2=0.0d0   ! max power of distribution 2
      pave1=0.0d0   ! ave power of distribution 1
      pave2=0.0d0   ! ave power of distribution 2
      do ip=1, npin
        do jp=1, npin
          do k=1, kd
            do na=1, nassm
              p1=power1(na,k,jp,ip)
              p2=power2(na,k,jp,ip)
              if (p1.gt.0.0d0) then
                np=np+1
                pave1=pave1+p1
                pave2=pave2+p2
                if (p1.gt.pmax1) then
                  pmax1=p1
                  lmax1(1)=ip
                  lmax1(2)=jp
                  lmax1(3)=k
                  lmax1(4)=na
                endif
                if (p2.gt.pmax2) then
                  pmax2=p2
                  lmax2(1)=ip
                  lmax2(2)=jp
                  lmax2(3)=k
                  lmax2(4)=na
                endif
                dd=abs(p2-p1)
                rms=rms+dd*dd
                if (dd.gt.dmax) then
                  dmax=dd
                  kmax1(1)=ip
                  kmax1(2)=jp
                  kmax1(3)=k
                  kmax1(4)=na
                endif
              endif
            enddo
          enddo
        enddo
      enddo

      if (np.gt.0) then
        rms=sqrt(rms/dble(np))
        pave1=pave1/dble(np)
        pave2=pave2/dble(np)
      endif

      xmax=dmax*100.0d0
      xrms=rms*100.0d0

      deallocate (power2)
      deallocate (power1)

      write (*,'(/1x,a)') 'Pin Power comparisons'
      write (*,120) np
      write (*,126) 1, pmax1, lmax1(:)
      write (*,126) 2, pmax2, lmax2(:)
      write (*,125) 1, pave1
      write (*,125) 2, pave2
      write (*,131) dmax*100.0d0, kmax1(:)
      write (*,130) rms*100.0d0
      write (*,138) (pmax2-pmax1)*100.0d0
 120  format (3x,'number values',2x,i10)
 125  format (3x,'ave power',i1,5x, f10.4,'   (does not include volume weighting or qtr sym)')
 126  format (3x,'max power',i1,5x, f10.4,'   at (ip,jp,k,na) ', 4i4)
 130  format (3x,'rms difference ', f10.4,' %')
 131  format (3x,'max difference ', f10.4,' % at (ip,jp,k,na) ', 4i4)
 138  format (3x,'difference in max pin ', f10.4,' %')

      return
      end subroutine pin_compare
!=======================================================================
!
!  Subroutine to count number of statepoints on a file and read some scalars
!
!=======================================================================
      subroutine countstates(file_id, maxstate, nstate, xexpo, xefpd, xkeff, xboron)
      use mod_hdftools, only : hdf5_read_double
      use hdf5
      implicit none
      integer(hid_t), intent(in)  :: file_id     ! HDF file id's
      integer,        intent(in)  :: maxstate    ! max number of statepoints on file
      integer,        intent(out) :: nstate      ! number of statepoints on file
      real(8)                     :: xexpo(maxstate)
      real(8)                     :: xefpd(maxstate)
      real(8)                     :: xkeff(maxstate)
      real(8)                     :: xboron(maxstate)

      logical :: ifxst
      integer :: ierror
      character(len=80) :: dataset
      character(len=10) :: statename

      nstate=0
      statename='STATE_0001'

      do
        write (statename(7:10),'(i4.4)') nstate+1
        call h5lexists_f(file_id, statename, ifxst, ierror)
        if (.not.ifxst) exit
        nstate=nstate+1
        if (nstate.gt.maxstate) stop 'maxstate exceeded'

        dataset=trim(statename)//'/exposure'
        call hdf5_read_double(file_id, dataset, xexpo(nstate))

        dataset=trim(statename)//'/exposure_efpd'
        call hdf5_read_double(file_id, dataset, xefpd(nstate))

        dataset=trim(statename)//'/keff'
        call hdf5_read_double(file_id, dataset, xkeff(nstate))

        dataset=trim(statename)//'/boron'
        call hdf5_read_double(file_id, dataset, xboron(nstate))

      enddo

      return
      end subroutine countstates

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

max_stars_repo_name: palmtag/VERAout-tools

max_issues_repo_name: palmtag/VERAout-tools

max_forks_repo_name: palmtag/VERAout-tools


text:
! wind_cdhf_alloc_def.for - wind/waves telemetry words, values are minor frame
! word numbers (0...255)
! Note: file wind_cdhf_rec_def.for must be included prior to including this
! file.

	integer*4	fixed_col_science_alloc( size_science_allocation)
	1		/17, 18, 24, 28, 32, 36, 40, 44, 48, 52, 56,
	1		 60, 64, 68, 72, 76, 80, 84, 88, 92, 96, 100, 104,
	1		108, 112, 116, 120, 124, 128, 132, 136, 140, 144,
	1		148, 152,
	1		156, 160, 164, 168, 172, 176, 180, 184, 188, 192/

	integer*4	fixed_col_maneuver_alloc(size_maneuver_allocation)
	1		/17, 18, 32, 33, 36, 37, 40, 41, 44, 45, 48, 49, 52,
	1		 53, 56, 57, 60, 61, 64, 65, 68, 69, 72, 73, 76, 77,
	1		 80, 81, 84, 85, 88, 89, 92, 93, 96, 97, 100, 101,
	1		104, 108, 112, 116, 120, 124, 128, 132, 136, 140, 144,
	1		148,
	1		152, 156, 160, 164, 168, 172, 176, 180, 184, 188, 192/

ext: for


text:
      subroutine dqc25f(f,a,b,omega,integr,nrmom,maxp1,ksave,result,
     *   abserr,neval,resabs,resasc,momcom,chebmo)
c***begin prologue  dqc25f
c***date written   810101   (yymmdd)
c***revision date  830518   (yymmdd)
c***category no.  h2a2a2
c***keywords  integration rules for functions with cos or sin
c             factor, clenshaw-curtis, gauss-kronrod
c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
c***purpose  to compute the integral i=integral of f(x) over (a,b)
c            where w(x) = cos(omega*x) or w(x)=sin(omega*x) and to
c            compute j = integral of abs(f) over (a,b). for small value
c            of omega or small intervals (a,b) the 15-point gauss-kronro
c            rule is used. otherwise a generalized clenshaw-curtis
c            method is used.
c***description
c
c        integration rules for functions with cos or sin factor
c        standard fortran subroutine
c        double precision version
c
c        parameters
c         on entry
c           f      - double precision
c                    function subprogram defining the integrand
c                    function f(x). the actual name for f needs to
c                    be declared e x t e r n a l in the calling program.
c
c           a      - double precision
c                    lower limit of integration
c
c           b      - double precision
c                    upper limit of integration
c
c           omega  - double precision
c                    parameter in the weight function
c
c           integr - integer
c                    indicates which weight function is to be used
c                       integr = 1   w(x) = cos(omega*x)
c                       integr = 2   w(x) = sin(omega*x)
c
c           nrmom  - integer
c                    the length of interval (a,b) is equal to the length
c                    of the original integration interval divided by
c                    2**nrmom (we suppose that the routine is used in an
c                    adaptive integration process, otherwise set
c                    nrmom = 0). nrmom must be zero at the first call.
c
c           maxp1  - integer
c                    gives an upper bound on the number of chebyshev
c                    moments which can be stored, i.e. for the
c                    intervals of lengths abs(bb-aa)*2**(-l),
c                    l = 0,1,2, ..., maxp1-2.
c
c           ksave  - integer
c                    key which is one when the moments for the
c                    current interval have been computed
c
c         on return
c           result - double precision
c                    approximation to the integral i
c
c           abserr - double precision
c                    estimate of the modulus of the absolute
c                    error, which should equal or exceed abs(i-result)
c
c           neval  - integer
c                    number of integrand evaluations
c
c           resabs - double precision
c                    approximation to the integral j
c
c           resasc - double precision
c                    approximation to the integral of abs(f-i/(b-a))
c
c         on entry and return
c           momcom - integer
c                    for each interval length we need to compute the
c                    chebyshev moments. momcom counts the number of
c                    intervals for which these moments have already been
c                    computed. if nrmom.lt.momcom or ksave = 1, the
c                    chebyshev moments for the interval (a,b) have
c                    already been computed and stored, otherwise we
c                    compute them and we increase momcom.
c
c           chebmo - double precision
c                    array of dimension at least (maxp1,25) containing
c                    the modified chebyshev moments for the first momcom
c                    momcom interval lengths
c
c ......................................................................
c***references  (none)
c***routines called  d1mach,dgtsl,dqcheb,dqk15w,dqwgtf
c***end prologue  dqc25f
c
      double precision a,abserr,ac,an,an2,as,asap,ass,b,centr,chebmo,
     *  cheb12,cheb24,conc,cons,cospar,d,dabs,dcos,dsin,dqwgtf,d1,
     *  d1mach,d2,estc,ests,f,fval,hlgth,oflow,omega,parint,par2,par22,
     *  p2,p3,p4,resabs,resasc,resc12,resc24,ress12,ress24,result,
     *  sinpar,v,x
      integer i,iers,integr,isym,j,k,ksave,m,momcom,neval,maxp1,
     *  noequ,noeq1,nrmom
c
      dimension chebmo(maxp1,25),cheb12(13),cheb24(25),d(25),d1(25),
     *  d2(25),fval(25),v(28),x(11)
c
      external f,dqwgtf
c
c           the vector x contains the values cos(k*pi/24)
c           k = 1, ...,11, to be used for the chebyshev expansion of f
c
      data x(1) / 0.9914448613 7381041114 4557526928 563d0 /
      data x(2) / 0.9659258262 8906828674 9743199728 897d0 /
      data x(3) / 0.9238795325 1128675612 8183189396 788d0 /
      data x(4) / 0.8660254037 8443864676 3723170752 936d0 /
      data x(5) / 0.7933533402 9123516457 9776961501 299d0 /
      data x(6) / 0.7071067811 8654752440 0844362104 849d0 /
      data x(7) / 0.6087614290 0872063941 6097542898 164d0 /
      data x(8) / 0.5000000000 0000000000 0000000000 000d0 /
      data x(9) / 0.3826834323 6508977172 8459984030 399d0 /
      data x(10) / 0.2588190451 0252076234 8898837624 048d0 /
      data x(11) / 0.1305261922 2005159154 8406227895 489d0 /
c
c           list of major variables
c           -----------------------
c
c           centr  - mid point of the integration interval
c           hlgth  - half-length of the integration interval
c           fval   - value of the function f at the points
c                    (b-a)*0.5*cos(k*pi/12) + (b+a)*0.5, k = 0, ..., 24
c           cheb12 - coefficients of the chebyshev series expansion
c                    of degree 12, for the function f, in the
c                    interval (a,b)
c           cheb24 - coefficients of the chebyshev series expansion
c                    of degree 24, for the function f, in the
c                    interval (a,b)
c           resc12 - approximation to the integral of
c                    cos(0.5*(b-a)*omega*x)*f(0.5*(b-a)*x+0.5*(b+a))
c                    over (-1,+1), using the chebyshev series
c                    expansion of degree 12
c           resc24 - approximation to the same integral, using the
c                    chebyshev series expansion of degree 24
c           ress12 - the analogue of resc12 for the sine
c           ress24 - the analogue of resc24 for the sine
c
c
c           machine dependent constant
c           --------------------------
c
c           oflow is the largest positive magnitude.
c
c***first executable statement  dqc25f
      oflow = d1mach(2)
c
      centr = 0.5d+00*(b+a)
      hlgth = 0.5d+00*(b-a)
      parint = omega*hlgth
c
c           compute the integral using the 15-point gauss-kronrod
c           formula if the value of the parameter in the integrand
c           is small.
c
      if(dabs(parint).gt.0.2d+01) go to 10
      call dqk15w(f,dqwgtf,omega,p2,p3,p4,integr,a,b,result,
     *  abserr,resabs,resasc)
      neval = 15
      go to 170
c
c           compute the integral using the generalized clenshaw-
c           curtis method.
c
   10 conc = hlgth*dcos(centr*omega)
      cons = hlgth*dsin(centr*omega)
      resasc = oflow
      neval = 25
c
c           check whether the chebyshev moments for this interval
c           have already been computed.
c
      if(nrmom.lt.momcom.or.ksave.eq.1) go to 120
c
c           compute a new set of chebyshev moments.
c
      m = momcom+1
      par2 = parint*parint
      par22 = par2+0.2d+01
      sinpar = dsin(parint)
      cospar = dcos(parint)
c
c           compute the chebyshev moments with respect to cosine.
c
      v(1) = 0.2d+01*sinpar/parint
      v(2) = (0.8d+01*cospar+(par2+par2-0.8d+01)*sinpar/parint)/par2
      v(3) = (0.32d+02*(par2-0.12d+02)*cospar+(0.2d+01*
     *  ((par2-0.80d+02)*par2+0.192d+03)*sinpar)/parint)/(par2*par2)
      ac = 0.8d+01*cospar
      as = 0.24d+02*parint*sinpar
      if(dabs(parint).gt.0.24d+02) go to 30
c
c           compute the chebyshev moments as the solutions of a
c           boundary value problem with 1 initial value (v(3)) and 1
c           end value (computed using an asymptotic formula).
c
      noequ = 25
      noeq1 = noequ-1
      an = 0.6d+01
      do 20 k = 1,noeq1
        an2 = an*an
        d(k) = -0.2d+01*(an2-0.4d+01)*(par22-an2-an2)
        d2(k) = (an-0.1d+01)*(an-0.2d+01)*par2
        d1(k+1) = (an+0.3d+01)*(an+0.4d+01)*par2
        v(k+3) = as-(an2-0.4d+01)*ac
        an = an+0.2d+01
   20 continue
      an2 = an*an
      d(noequ) = -0.2d+01*(an2-0.4d+01)*(par22-an2-an2)
      v(noequ+3) = as-(an2-0.4d+01)*ac
      v(4) = v(4)-0.56d+02*par2*v(3)
      ass = parint*sinpar
      asap = (((((0.210d+03*par2-0.1d+01)*cospar-(0.105d+03*par2
     *  -0.63d+02)*ass)/an2-(0.1d+01-0.15d+02*par2)*cospar
     *  +0.15d+02*ass)/an2-cospar+0.3d+01*ass)/an2-cospar)/an2
      v(noequ+3) = v(noequ+3)-0.2d+01*asap*par2*(an-0.1d+01)*
     *   (an-0.2d+01)
c
c           solve the tridiagonal system by means of gaussian
c           elimination with partial pivoting.
c
c***        call to dgtsl has been replaced by call to
c***        lapack routine dgtsv
c
c      call dgtsl(noequ,d1,d,d2,v(4),iers)
      call dgtsv(noequ,1,d1(2:),d,d2(:noequ-1),v(4),noequ,iers)
      go to 50
c
c           compute the chebyshev moments by means of forward
c           recursion.
c
   30 an = 0.4d+01
      do 40 i = 4,13
        an2 = an*an
        v(i) = ((an2-0.4d+01)*(0.2d+01*(par22-an2-an2)*v(i-1)-ac)
     *  +as-par2*(an+0.1d+01)*(an+0.2d+01)*v(i-2))/
     *  (par2*(an-0.1d+01)*(an-0.2d+01))
        an = an+0.2d+01
   40 continue
   50 do 60 j = 1,13
        chebmo(m,2*j-1) = v(j)
   60 continue
c
c           compute the chebyshev moments with respect to sine.
c
      v(1) = 0.2d+01*(sinpar-parint*cospar)/par2
      v(2) = (0.18d+02-0.48d+02/par2)*sinpar/par2
     *  +(-0.2d+01+0.48d+02/par2)*cospar/parint
      ac = -0.24d+02*parint*cospar
      as = -0.8d+01*sinpar
      if(dabs(parint).gt.0.24d+02) go to 80
c
c           compute the chebyshev moments as the solutions of a boundary
c           value problem with 1 initial value (v(2)) and 1 end value
c           (computed using an asymptotic formula).
c
      an = 0.5d+01
      do 70 k = 1,noeq1
        an2 = an*an
        d(k) = -0.2d+01*(an2-0.4d+01)*(par22-an2-an2)
        d2(k) = (an-0.1d+01)*(an-0.2d+01)*par2
        d1(k+1) = (an+0.3d+01)*(an+0.4d+01)*par2
        v(k+2) = ac+(an2-0.4d+01)*as
        an = an+0.2d+01
   70 continue
      an2 = an*an
      d(noequ) = -0.2d+01*(an2-0.4d+01)*(par22-an2-an2)
      v(noequ+2) = ac+(an2-0.4d+01)*as
      v(3) = v(3)-0.42d+02*par2*v(2)
      ass = parint*cospar
      asap = (((((0.105d+03*par2-0.63d+02)*ass+(0.210d+03*par2
     *  -0.1d+01)*sinpar)/an2+(0.15d+02*par2-0.1d+01)*sinpar-
     *  0.15d+02*ass)/an2-0.3d+01*ass-sinpar)/an2-sinpar)/an2
      v(noequ+2) = v(noequ+2)-0.2d+01*asap*par2*(an-0.1d+01)
     *  *(an-0.2d+01)
c
c           solve the tridiagonal system by means of gaussian
c           elimination with partial pivoting.
c
c***        call to dgtsl has been replaced by call to
c***        lapack routine dgtsv
c
c      call dgtsl(noequ,d1,d,d2,v(3),iers)
      call dgtsv(noequ,1,d1(2:),d,d2(:noequ-1),v(3),noequ,iers)
      go to 100
c
c           compute the chebyshev moments by means of forward recursion.
c
   80 an = 0.3d+01
      do 90 i = 3,12
        an2 = an*an
        v(i) = ((an2-0.4d+01)*(0.2d+01*(par22-an2-an2)*v(i-1)+as)
     *  +ac-par2*(an+0.1d+01)*(an+0.2d+01)*v(i-2))
     *  /(par2*(an-0.1d+01)*(an-0.2d+01))
        an = an+0.2d+01
   90 continue
  100 do 110 j = 1,12
        chebmo(m,2*j) = v(j)
  110 continue
  120 if (nrmom.lt.momcom) m = nrmom+1
       if (momcom.lt.(maxp1-1).and.nrmom.ge.momcom) momcom = momcom+1
c
c           compute the coefficients of the chebyshev expansions
c           of degrees 12 and 24 of the function f.
c
      fval(1) = 0.5d+00*f(centr+hlgth)
      fval(13) = f(centr)
      fval(25) = 0.5d+00*f(centr-hlgth)
      do 130 i = 2,12
        isym = 26-i
        fval(i) = f(hlgth*x(i-1)+centr)
        fval(isym) = f(centr-hlgth*x(i-1))
  130 continue
      call dqcheb(x,fval,cheb12,cheb24)
c
c           compute the integral and error estimates.
c
      resc12 = cheb12(13)*chebmo(m,13)
      ress12 = 0.0d+00
      k = 11
      do 140 j = 1,6
        resc12 = resc12+cheb12(k)*chebmo(m,k)
        ress12 = ress12+cheb12(k+1)*chebmo(m,k+1)
        k = k-2
  140 continue
      resc24 = cheb24(25)*chebmo(m,25)
      ress24 = 0.0d+00
      resabs = dabs(cheb24(25))
      k = 23
      do 150 j = 1,12
        resc24 = resc24+cheb24(k)*chebmo(m,k)
        ress24 = ress24+cheb24(k+1)*chebmo(m,k+1)
        resabs = dabs(cheb24(k))+dabs(cheb24(k+1))
        k = k-2
  150 continue
      estc = dabs(resc24-resc12)
      ests = dabs(ress24-ress12)
      resabs = resabs*dabs(hlgth)
      if(integr.eq.2) go to 160
      result = conc*resc24-cons*ress24
      abserr = dabs(conc*estc)+dabs(cons*ests)
      go to 170
  160 result = conc*ress24+cons*resc24
      abserr = dabs(conc*ests)+dabs(cons*estc)
  170 return
      end

max_stars_repo_name: jeanpauphilet/scipy

max_issues_repo_name: jeanpauphilet/scipy

max_forks_repo_name: jeanpauphilet/scipy


text:
!----------------------------------*-F90-*----------------------------------//
!!
!* \file   FortranChecks/f90sub/mpi_hw_ftest.f90
!* \author Allan Wollaber, Tom Evans, Kelly Thompson
!* \date   Fri Mar 15 11:25:00 2002
!* \brief  Infinite medium Draco shunt.
!* \note   Copyright (C) 2016-2019 Triad National Security, LLC.
!*         All rights reserved.
!*
!---------------------------------------------------------------------------//
! $Id$
!---------------------------------------------------------------------------//

subroutine tst_mpi_hw_f(nf) bind(C, name="tst_mpi_hw")

  use draco_mpi
  use draco_test

  use iso_c_binding, only : c_int, c_double

  implicit none

  integer(c_int), intent(inout) :: nf
  integer :: ierr

  ! ------------------------------------
  ! Initialize MPI
  ! ------------------------------------
  call f90_mpi_init(ierr)
  call check_fail(ierr, f90_rank)

  ! ------------------------------------
  ! Run the problem
  ! ------------------------------------

  if ( f90_rank < f90_num_ranks ) then
     call pass_msg( f90_rank, "MPI rank index ok" )
  else
     call it_fails( f90_rank, "MPI rank > max" )
  endif

  call f90_mpi_barrier(ierr)
  print *, "Hello from rank ", f90_rank, "/", f90_num_ranks

  ! ------------------------------------
  ! Finalize and clean up
  ! ------------------------------------

  call f90_mpi_finalize(ierr)
  call check_fail(ierr, f90_rank)

  ! Print the overall test result
  ! call test_report(f90_rank,nf)

  nf = nf + f90_num_failures

end subroutine tst_mpi_hw_f

max_stars_repo_name: rspavel/Draco

max_issues_repo_name: rspavel/Draco

max_forks_repo_name: rspavel/Draco