text:
!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!>  This module contains subroutines used in 3d eigen analysis for
module m_dynamic_mass

contains

  subroutine mass_C2(etype, nn, ecoord, gausses, sec_opt, thick, mass, lumped, temperature)
    use mMechGauss
    use m_MatMatrix
    use elementInfo
    implicit none
    type(tGaussStatus), intent(in) :: gausses(:)             !< status of qudrature points
    integer(kind=kint), intent(in) :: etype                  !< element type
    integer(kind=kint), intent(in) :: nn                     !< number of elemental nodes
    real(kind=kreal), intent(in)  :: ecoord(2,nn)           !< coordinates of elemental nodes
    real(kind=kreal), intent(out) :: mass(:,:)              !< mass matrix
    real(kind=kreal), intent(out) :: lumped(:)              !< mass matrix
    real(kind=kreal), intent(in), optional :: temperature(nn) !< temperature
    type(tMaterial), pointer :: matl !< material information
    integer(kind=kint), parameter :: ndof = 2
    integer(kind=kint) :: i, j, LX, sec_opt
    real(kind=kreal) :: naturalCoord(2)
    real(kind=kreal) :: func(nn), thick
    real(kind=kreal) :: det, wg, rho
    real(kind=kreal) :: D(2,2), N(2, nn*ndof), DN(2, nn*ndof)
    real(kind=kreal) :: gderiv(nn,2)
    logical :: is_lumped

    mass(:,:) = 0.0d0
    lumped = 0.0d0
    matl => gausses(1)%pMaterial

    if(sec_opt == 2) thick = 1.0d0

    do LX = 1, NumOfQuadPoints(etype)
      call getQuadPoint(etype, LX, naturalCoord)
      call getShapeFunc(etype, naturalCoord, func)
      call getGlobalDeriv(etype, nn, naturalcoord, ecoord, det, gderiv)

      if(present(temperature))then
        !ina(1) = temperature
        !call fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr, ina)
        !if(ierr)then
          rho = matl%variables(M_DENSITY)
        !else
        !  rho = outa(1)
        !endif
      else
        !call fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr)
        !if(ierr)then
          rho = matl%variables(M_DENSITY)
        !else
        !  rho = outa(1)
        !endif
      endif

      D = 0.0d0
      D(1,1) = rho*thick
      D(2,2) = rho*thick

      wg = getWeight(etype, LX)*det

      N = 0.0d0
      do i = 1, nn
        N(1,2*i-1) = func(i)
        N(2,2*i  ) = func(i)
      enddo

      DN(1:2, 1:nn*ndof) = matmul(D, N(1:2, 1:nn*ndof))
      forall(i = 1:nn*ndof,  j = 1:nn*ndof)
        mass(i,j) = mass(i,j) + dot_product(N(:,i), DN(:,j))*wg
      end forall
    enddo

    is_lumped = .true.
    if(is_lumped) call get_lumped_mass(nn, ndof, mass, lumped)
  end subroutine mass_C2

  subroutine mass_C3(etype, nn, ecoord, gausses, mass, lumped, temperature)
    use mMechGauss
    use m_MatMatrix
    use elementInfo
    implicit none
    type(tGaussStatus), intent(in) :: gausses(:)             !< status of qudrature points
    integer(kind=kint), intent(in) :: etype                  !< element type
    integer(kind=kint), intent(in) :: nn                     !< number of elemental nodes
    real(kind=kreal), intent(in)  :: ecoord(3,nn)           !< coordinates of elemental nodes
    real(kind=kreal), intent(out) :: mass(:,:)              !< mass matrix
    real(kind=kreal), intent(out) :: lumped(:)              !< mass matrix
    real(kind=kreal), intent(in), optional :: temperature(nn) !< temperature
    type(tMaterial), pointer :: matl !< material information
    integer(kind=kint), parameter :: ndof = 3
    integer(kind=kint) :: i, j, LX
    real(kind=kreal) :: naturalCoord(3)
    real(kind=kreal) :: func(nn)
    real(kind=kreal) :: det, wg, rho
    real(kind=kreal) :: D(3, 3), N(3, nn*ndof), DN(3, nn*ndof)
    real(kind=kreal) :: gderiv(nn, 3)
    logical :: is_lumped

    mass(:,:) = 0.0d0
    lumped = 0.0d0
    matl => gausses(1)%pMaterial

    do LX = 1, NumOfQuadPoints(etype)
      call getQuadPoint(etype, LX, naturalCoord)
      call getShapeFunc(etype, naturalCoord, func)
      call getGlobalDeriv(etype, nn, naturalcoord, ecoord, det, gderiv)

      if(present(temperature))then
        !ina(1) = temperature
        !call fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr, ina)
        !if(ierr)then
          rho = matl%variables(M_DENSITY)
        !else
        !  rho = outa(1)
        !endif
      else
        !call fetch_TableData(MC_ISOELASTIC, matl%dict, outa, ierr)
        !if(ierr)then
          rho = matl%variables(M_DENSITY)
        !else
        !  rho = outa(1)
        !endif
      endif

      D = 0.0d0
      D(1,1) = rho
      D(2,2) = rho
      D(3,3) = rho

      wg = getWeight(etype,LX)*det

      N = 0.0d0
      do i = 1, nn
        N(1,3*i-2) = func(i)
        N(2,3*i-1) = func(i)
        N(3,3*i  ) = func(i)
      enddo

      DN(1:3, 1:nn*ndof) = matmul(D, N(1:3, 1:nn*ndof))
      forall(i = 1:nn*ndof,  j = 1:nn*ndof)
        mass(i,j) = mass(i,j) + dot_product(N(:,i), DN(:,j))*wg
      end forall
    enddo

    is_lumped = .true.
    if(is_lumped) call get_lumped_mass(nn, ndof, mass, lumped)
  end subroutine mass_C3

  subroutine get_lumped_mass(nn, ndof, mass, lumped)
    use hecmw
    implicit none
    integer(kind=kint) :: i, j, nn, ndof
    real(kind=kreal) :: lumped(:), mass(:,:)
    real(kind=kreal) :: diag_mass, total_mass

    total_mass = 0.0d0
    do i = 1, nn*ndof, ndof
      do j = 1, nn*ndof, ndof
        total_mass = total_mass + mass(j,i)
      enddo
    enddo

    diag_mass = 0.0d0
    do i = 1, nn*ndof, ndof
      diag_mass = diag_mass + mass(i,i)
    enddo

    diag_mass = 1.0d0/diag_mass
    do i = 1, nn*ndof
      lumped(i) = lumped(i) + mass(i,i)*total_mass*diag_mass
    enddo

    mass = 0.0d0
    do i = 1, nn*ndof
      mass(i,i) = lumped(i)
    enddo
  end subroutine get_lumped_mass

  function get_length(ecoord)
    use hecmw
    implicit none
    real(kind=kreal) :: get_length, ecoord(3,20)

    get_length = dsqrt( &
      (ecoord(1,2) - ecoord(1,1))**2 + &
      (ecoord(2,2) - ecoord(2,1))**2 + &
      (ecoord(3,2) - ecoord(3,1))**2 )
  end function get_length

  function get_face3(ecoord)
    use hecmw
    implicit none
    real(kind=kreal) :: get_face3, ecoord(3,20)
    real(kind=kreal) :: a1, a2, a3
    real(kind=kreal) :: X(3), Y(3), Z(3)

    X(1) = ecoord(1,1); Y(1) = ecoord(2,1); Z(1) = ecoord(3,1)
    X(2) = ecoord(1,2); Y(2) = ecoord(2,2); Z(2) = ecoord(3,2)
    X(3) = ecoord(1,3); Y(3) = ecoord(2,3); Z(3) = ecoord(3,3)

    a1 = (X(2) - X(1))**2 + (Y(2) - Y(1))**2 + (Z(2) - Z(1))**2
    a2 = (X(1) - X(3))*(X(2) - X(1)) &
     & + (Y(1) - Y(3))*(Y(2) - Y(1)) &
     & + (Z(1) - Z(3))*(Z(2) - Z(1))
    a3 = (X(3) - X(1))**2 + (Y(3) - Y(1))**2 + (Z(3) - Z(1))**2

    get_face3 = 0.5d0*dsqrt(a1*a3 - a2*a2)
  end function get_face3

  function get_face4(ecoord)
    use hecmw
    implicit none
    integer(kind=kint) :: LX, LY
    real(kind=kreal) :: get_face4, ecoord(3,20)
    real(kind=kreal) :: XG(2), RI, SI, RP, SP, RM, SM, HR(4), HS(4)
    real(kind=kreal) :: XR, XS, YR, YS, ZR, ZS
    real(kind=kreal) :: X(4), Y(4), Z(4), det

    X(1) = ecoord(1,1); Y(1) = ecoord(2,1); Z(1) = ecoord(3,1)
    X(2) = ecoord(1,2); Y(2) = ecoord(2,2); Z(2) = ecoord(3,2)
    X(3) = ecoord(1,3); Y(3) = ecoord(2,3); Z(3) = ecoord(3,3)
    X(4) = ecoord(1,4); Y(4) = ecoord(2,4); Z(4) = ecoord(3,4)

    XG(1) = -0.5773502691896258D0
    XG(2) = -XG(1)
    get_face4 = 0.0d0

    do LX = 1, 2
      RI = XG(LX)
      do LY = 1, 2
        SI = XG(LY)
        RP = 1.0d0 + RI
        SP = 1.0d0 + SI
        RM = 1.0d0 - RI
        SM = 1.0d0 - SI

        !C*  FOR R-COORDINATE
        HR(1) =  0.25d0*SP
        HR(2) = -0.25d0*SP
        HR(3) = -0.25d0*SM
        HR(4) =  0.25d0*SM

        !C*  FOR S-COORDINATE
        HS(1) =  0.25d0*RP
        HS(2) =  0.25d0*RM
        HS(3) = -0.25d0*RM
        HS(4) = -0.25d0*RP

        !C*JACOBI MATRIX
        XR = HR(1)*X(1) + HR(2)*X(2) + HR(3)*X(3) + HR(4)*X(4)
        XS = HS(1)*X(1) + HS(2)*X(2) + HS(3)*X(3) + HS(4)*X(4)
        YR = HR(1)*Y(1) + HR(2)*Y(2) + HR(3)*Y(3) + HR(4)*Y(4)
        YS = HS(1)*Y(1) + HS(2)*Y(2) + HS(3)*Y(3) + HS(4)*Y(4)
        ZR = HR(1)*Z(1) + HR(2)*Z(2) + HR(3)*Z(3) + HR(4)*Z(4)
        ZS = HS(1)*Z(1) + HS(2)*Z(2) + HS(3)*Z(3) + HS(4)*Z(4)

        det = (YR*ZS - ZR*YS)**2 + (ZR*XS - XR*ZS)**2 + (XR*YS - YR*XS)**2
        det = dsqrt(det)

        get_face4 = get_face4 + det
      enddo
    enddo
  end function get_face4

end module m_dynamic_mass

max_stars_repo_name: fujitsu/FrontISTR

max_issues_repo_name: fujitsu/FrontISTR

max_forks_repo_name: fujitsu/FrontISTR


text:
**==dcstrg.spg  processed by SPAG 5.00N  at 11:12 on 23 Mar 1999
      SUBROUTINE DCSTRG(Str,Lstr,Ityp)
      IMPLICIT NONE
C     Note: Minor update by B. Melhus on 21.04.2021
C ======================================================================
C  PURPOSE:
C       Decode char. string according to criterion given by ITYP.
C       (See below).
C
C  METHOD:
C       - CALL MOVEIP to move line pointer IP to first non-blank
C         character in line buffer CHBUFF.
C       - Decode string from CHBUFF.
C       - After decoding of string, IP will be the first blank
C         after string.
C
C  ERROR HANDLING:
C       - If call to MOVEIP return error flag on, then error message
C         is printed.
C       - If invalid char. then error message is printed.
C       - If invalid char. IP is moved to the next blank character in
C         CHBUFF.
C
C  INPUT:
C       ITYP: Type of char. string to be decoded.
C             1: Decode any string until first blank char.
C             2: Decode any string between ' '.
C                Return STR and LSTR. If string enclosed by ' '
C                can not be decoded as next data-item, error-flag
C                is set to a negative value.
C             3: Decode string until first blank character.
C                Return STR and LSTR. Nearly equal to ITYP=1,
C                but call ALPHA(3) before any decoding.
C             4: Decode any string with length LSTR from first
C                non-blank character.
C                Return STR only.
C       LSTR: If ITYP = 1 - 4: Max length of string to be decoded is
C             verified.
C             (LEN() is used to calc. length of STR). If string to
C             be decoded is outside max limit: decoding is stopped.
C             IDCERR is set to a value < 0.
C
C
C  OUTPUT:
C       STR   - Decoded string.
C       LSTR  - Length of STR.
C
C  INTERNAL VARIABLES:
C
C  CALL TO:
C       MOVEIP
C
C  CALL FROM:
C       General purpose routine.
C
C  LIBRARY:
C       FFDLIB
C
C  MACHINE DEPENDENT PART:
C       INCLUDE statement
C
C  PROGRAMMED BY:
C       Per Ivar Loekstad / ETOE / Aker Engineering a.s.
C
C  DATE/VERSION:
C       19.02.92/1.0
C
C ======================================================================
C
      INTEGER iend , Ityp , maxwid
      INTEGER Lstr , istop , j , itell
      CHARACTER Str*(*) , ch
C
C Function
C
      LOGICAL ALPHA
C
      INCLUDE 'ffd.inc'
C
C ----------------------------------------------------------------------
C
C Initialize
C
      IF ( Ityp.NE.4 ) THEN
         maxwid = LEN(Str)
         Lstr = 0
      ENDIF
      Str = ' '
      istop = 0
C
C     New code due to run time error experienced by Intel Visual Fortran
C     Error: The variable 'J' is being used without being defined at line 193.
C     The parameter is now initialized to zero. B. Melhus on 21.04.2021
      j = 0
C
C Move IP to first non-blank character.
C
      CALL MOVEIP(iend)
      IF ( IDCerr.NE.istop ) THEN
C
         WRITE (LPU,99001)
         IDCerr = -1
      ELSEIF ( iend.EQ.0 ) THEN
C ----------------------------------------------------------------------
C Decode string
C
         IF ( Ityp.EQ.1 ) THEN
C
            Lstr = 1
            j = 1
 20         IF ( IP.LE.NCOl ) THEN
               ch = CHBuff(IP:IP)
               IF ( ch.EQ.' ' ) THEN
C
 25               IF ( IP.LE.NCOl ) THEN
                     IF ( CHBuff(IP:IP).NE.' ' ) THEN
                        IP = IP + 1
                        GOTO 25
                     ENDIF
                  ENDIF
               ELSE
                  IP = IP + 1
                  IF ( j.GT.maxwid ) GOTO 100
                  Str(j:j) = ch
                  j = j + 1
                  GOTO 20
               ENDIF
            ENDIF
            Lstr = j - 1
C
C ----------------------------------------------------------------------
C
            ELSEIF ( Ityp.EQ.2 ) THEN
C
            IF ( ALPHA(2) ) THEN
               j = 1
               itell = 0
 30            IF ( IP.LE.NCOl ) THEN
                  IF ( .NOT.CHBuff(IP:IP).EQ.'''' ) ch = CHBuff(IP:IP)
                  IF ( CHBuff(IP:IP).EQ.'''' ) THEN
                     itell = itell + 1
                     IF ( itell.EQ.1 ) Lstr = IP
                     IP = IP + 1
                     IF ( itell.EQ.1 ) GOTO 30
                  ENDIF
                  IP = IP + 1
                  IF ( itell.EQ.2 ) THEN
                     Lstr = IP - Lstr - 3
                     GOTO 200
                  ENDIF
                  IF ( j.GT.maxwid ) GOTO 100
                  Str(j:j) = ch
                  j = j + 1
                  GOTO 30
               ELSE
                  IDCerr = -1
                  WRITE (LPU,99003) NCOl
               ENDIF
            ELSE
               IDCerr = -1
               WRITE (LPU,99002)
            ENDIF
C
C ----------------------------------------------------------------------
C
         ELSEIF ( Ityp.EQ.3 ) THEN
C
C           Minor update: Error encountered with Visual Studio
C                         Use of uninitalized variable j
C                         B. Melhus by 21.04.2021
            j=0
C
            IF ( ALPHA(3) ) THEN
               Lstr = 1
               j = 1
 40            IF ( IP.LE.NCOl ) THEN
                  ch = CHBuff(IP:IP)
                  IF ( ch.EQ.' ' ) THEN
C
 42                  IF ( IP.LE.NCOl ) THEN
                        IF ( CHBuff(IP:IP).NE.' ' ) THEN
                           IP = IP + 1
                           GOTO 42
                        ENDIF
                     ENDIF
                  ELSE
                     IP = IP + 1
                     IF ( j.GT.maxwid ) GOTO 100
                     Str(j:j) = ch
                     j = j + 1
                     GOTO 40
                  ENDIF
               ENDIF
            ELSE
               IDCerr = -1
               WRITE (LPU,99002)
            ENDIF
            Lstr = j - 1
C
C ----------------------------------------------------------------------
C
         ELSEIF ( Ityp.EQ.4 ) THEN
C
            IF ( LEN(Str).LT.Lstr ) THEN
               IDCerr = -1
               WRITE (LPU,99004) LEN(Str)
               GOTO 200
            ENDIF
            IF ( IP.LE.NCOl ) THEN
               DO 50 j = 1 , Lstr
                  IF ( IP.LE.NCOl ) Str(j:j) = CHBuff(IP:IP)
                  IP = IP + 1
 50            CONTINUE
            ENDIF
C           LSTR = 1
C
C ----------------------------------------------------------------------
C
         ENDIF
      ENDIF
      GOTO 200
C
 100  IDCerr = -1
      WRITE (LPU,99004) maxwid
      IP = IP - 1
      Lstr = maxwid
C
 200  RETURN
C ----------------------------------------------------------------------
99001 FORMAT (/,' *** FFDLIB ERROR - Call to MOVEIP - (MODULE: DCSTRG)')
99002 FORMAT (/,' *** FFDLIB ERROR - Illegal character -',
     &        ' (MODULE: DCSTRG)')
99003 FORMAT (/,' *** FFDLIB ERROR - Try to read outside max.',
     &        ' line-width, NCOL: ',I3,' - (MODULE: DCSTRG)')
99004 FORMAT (/,' *** FFDLIB ERROR - Length of string exceed max.',
     &        ' length given as: ',I3,' - (MODULE: DCSTRG)')
C ----------------------------------------------------------------------
      END

ext: f


! Check that failing an image works.
! Image two is to fail, all others to continue.

program image_fail_test_1
  use iso_fortran_env , only : STAT_FAILED_IMAGE
  implicit none
  integer :: i, syncAllStat

  associate(np => num_images(), me => this_image())
    if (np < 3) error stop "I need at least 3 images to function."

    if (me == 2) fail image

    if (me == 2) print *, "Test failed."

    sync all (STAT=syncAllStat)
    if (me == 1) print *, "Test passed."
  end associate

end program image_fail_test_1


text:
      SUBROUTINE loop100_F90(N, a, b, c, d, x, u, v, w, z)
      INTEGER i, N
      REAL*8 a(N), b(N), c(N), d(N), x(N), u, v, w, z

      x=(1.0-c*c)/((4*w)*sin(1.0+c*c-2*v*c))*a*b*u*exp(-z*d)
      RETURN
      END

      SUBROUTINE loop100_F90Overhead(N, a, b, c, d, x, u, v, w, z)
      INTEGER i, N
      REAL*8 a(N), b(N), c(N), d(N), x(N), u, v, w, z

      RETURN
      END

max_stars_repo_name: fraguela/depspawn

max_issues_repo_name: MSV-Project/IBAMR

max_forks_repo_name: MSV-Project/IBAMR


text:
      subroutine qqbw2j_loop(i1,i2,i3,i4,i5,i6,
     .                       qqb_ijkk,qqb_iikl,qqb_ijkj,qqb_ijik,
     .                       qqb_ijii,qqb_ijjj,qqb_iiij,qqb_iiji)
      implicit none
      include 'constants.f'
      include 'zprods_com.f'
      integer i1,i2,i3,i4,i5,i6
      double complex a61_1234(2),a61_3412(2),a61_1432(2),a61_3214(2),
     .               a62_1234(2),a62_3412(2),a62_1432(2),a62_3214(2),
     .               atr_1234(2),atr_3412(2),atr_1432(2),atr_3214(2)
      double complex atrLLL,atrLRL,a61LLL,a61LRL,a62LLL,a62LRL
      double precision qqb_ijkk,qqb_iikl,qqb_ijkj,qqb_ijik,
     .                 qqb_ijii,qqb_ijjj,qqb_iiij,qqb_iiji
c--- this routine returns various pieces of the tree-loop interference
c--- for the process 0 -> q(1) qb(4) Q(3) Qb(2) + W (-> l(5) + lbar(6))
c--- the returned pieces are labelled according to the qqb -> qqb
c--- crossing and are calculated in the FORM program sq4q, with
c--- reference to the lowest order qqb_w2jetx.f

c--- there are 2 helicity amplitudes, we label LLL as 1, LRL as 2
      a61_1234(1)=a61LLL(i1,i2,i3,i4,i5,i6,za,zb)
      a61_1234(2)=a61LRL(i1,i2,i3,i4,i5,i6,za,zb)

      a61_3412(1)=a61LLL(i3,i4,i1,i2,i5,i6,za,zb)
      a61_3412(2)=a61LRL(i3,i4,i1,i2,i5,i6,za,zb)

      a61_1432(1)=a61LLL(i1,i4,i3,i2,i5,i6,za,zb)
      a61_1432(2)=a61LRL(i1,i4,i3,i2,i5,i6,za,zb)

      a61_3214(1)=a61LLL(i3,i2,i1,i4,i5,i6,za,zb)
      a61_3214(2)=a61LRL(i3,i2,i1,i4,i5,i6,za,zb)

      a62_1234(1)=a62LLL(i1,i2,i3,i4,i5,i6,za,zb)
      a62_1234(2)=a62LRL(i1,i2,i3,i4,i5,i6,za,zb)

      a62_3412(1)=a62LLL(i3,i4,i1,i2,i5,i6,za,zb)
      a62_3412(2)=a62LRL(i3,i4,i1,i2,i5,i6,za,zb)

      a62_1432(1)=a62LLL(i1,i4,i3,i2,i5,i6,za,zb)
      a62_1432(2)=a62LRL(i1,i4,i3,i2,i5,i6,za,zb)

      a62_3214(1)=a62LLL(i3,i2,i1,i4,i5,i6,za,zb)
      a62_3214(2)=a62LRL(i3,i2,i1,i4,i5,i6,za,zb)

      atr_1234(1)=Dconjg(atrLLL(i1,i2,i3,i4,i5,i6,za,zb))
      atr_1234(2)=Dconjg(atrLRL(i1,i2,i3,i4,i5,i6,za,zb))

      atr_3412(1)=Dconjg(atrLLL(i3,i4,i1,i2,i5,i6,za,zb))
      atr_3412(2)=Dconjg(atrLRL(i3,i4,i1,i2,i5,i6,za,zb))

      atr_1432(1)=Dconjg(atrLLL(i1,i4,i3,i2,i5,i6,za,zb))
      atr_1432(2)=Dconjg(atrLRL(i1,i4,i3,i2,i5,i6,za,zb))

      atr_3214(1)=Dconjg(atrLLL(i3,i2,i1,i4,i5,i6,za,zb))
      atr_3214(2)=Dconjg(atrLRL(i3,i2,i1,i4,i5,i6,za,zb))

* Note: in the expressions below, there are no LRL (2) interference
*       terms, since these would correspond to polarizations
*       that don't match, ie. (LL -> RR) x (LR -> LR)
      qqb_ijkk=dble(a61_1234(1)*atr_1234(1)+a61_1234(2)*atr_1234(2))
      qqb_iikl=dble(a61_3412(1)*atr_3412(1)+a61_3412(2)*atr_3412(2))
      qqb_ijkj=dble(a61_3214(1)*atr_3214(1)+a61_3214(2)*atr_3214(2))
      qqb_ijik=dble(a61_1432(1)*atr_1432(1)+a61_1432(2)*atr_1432(2))
      qqb_ijii=dble(a61_1234(1)*atr_1234(1)+a61_1234(2)*atr_1234(2)
     .             +a61_1432(1)*atr_1432(1)+a61_1432(2)*atr_1432(2)
     .            -(a62_1234(1)*atr_1432(1)+a62_1432(1)*atr_1234(1))/xn)
      qqb_ijjj=dble(a61_1234(1)*atr_1234(1)+a61_1234(2)*atr_1234(2)
     .             +a61_3214(1)*atr_3214(1)+a61_3214(2)*atr_3214(2)
     .            -(a62_1234(1)*atr_3214(1)+a62_3214(1)*atr_1234(1))/xn)
      qqb_iiij=dble(a61_1432(1)*atr_1432(1)+a61_1432(2)*atr_1432(2)
     .             +a61_3412(1)*atr_3412(1)+a61_3412(2)*atr_3412(2)
     .            -(a62_1432(1)*atr_3412(1)+a62_3412(1)*atr_1432(1))/xn)
      qqb_iiji=dble(a61_3214(1)*atr_3214(1)+a61_3214(2)*atr_3214(2)
     .             +a61_3412(1)*atr_3412(1)+a61_3412(2)*atr_3412(2)
     .            -(a62_3214(1)*atr_3412(1)+a62_3412(1)*atr_3214(1))/xn)

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text: #redirect Friends Urging Campus Kindness, Education and Diversity

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

!>    @brief OpenMP wrapper for "omp_neumann_temp"
!>    @param[out] neumann_max neumann criteria
!>    @param[in] ismpl local sample index
      SUBROUTINE neumann_temp(neumann_max,ismpl)
        use mod_genrl
        use mod_OMP_TOOLS
        IMPLICIT NONE
        integer :: ismpl

        INCLUDE 'OMP_TOOLS.inc'
        DOUBLE PRECISION neumann_max

#ifdef fOMP
!$OMP parallel num_threads(Tlevel_1)
!$      call omp_binding(ismpl)
#endif
        CALL omp_neumann_temp(neumann_max,ismpl)
#ifdef fOMP
!$OMP end parallel
#endif

        RETURN
      END

!>    @brief calculate grid neuman numbers (temp)
!>    @param[out] neumann_max maximal neuman number
!>    @param[in] ismpl local sample index
      SUBROUTINE omp_neumann_temp(neumann_max,ismpl)
        use arrays
        use mod_genrl
        use mod_genrlc
        use mod_temp
        use mod_time
        use mod_linfos
        IMPLICIT NONE
        integer :: ismpl
        integer :: i, j, k

        INTEGER c1, c2, c3
        DOUBLE PRECISION neumann_maxx, neumann_minx, neumann_avgx
        DOUBLE PRECISION neumann_maxy, neumann_miny, neumann_avgy
        DOUBLE PRECISION neumann_maxz, neumann_minz, neumann_avgz
        DOUBLE PRECISION val, neumann_max, delt, fac, davg
        DOUBLE PRECISION deltat
        EXTERNAL deltat
        DOUBLE PRECISION li, lj, lk, rhocf, por, rhoceff
        EXTERNAL li, lj, lk, rhocf, por, rhoceff

        delt = deltat(simtime(ismpl),ismpl)

        IF ( .NOT. (transient .AND. tr_switch(ismpl))) THEN
!$OMP master
          WRITE(*,*) ' neumann-temp: not defined for steady state'
!$OMP end master
          RETURN
        ELSE IF (linfos(3)>=2) THEN
!$OMP master
          WRITE(*,*)
          WRITE(*,'(A,1e16.8)') '  ... neumann-temp:  delt/tunit = ', &
            delt/tunit
          WRITE(*,*)
!$OMP end master
        END IF

! val in x
        c1 = 0
        neumann_maxx = small
        neumann_minx = big
        neumann_avgx = 0.0D0
!$OMP do schedule(static)
        DO k = 1, k0
          DO j = 1, j0
            DO i = 2, i0 - 1
              c1 = c1 + 1
              davg = 0.5D0*(delx(i)+delx(i+1))
              fac = delt/rhoceff(i,j,k,ismpl)
              val = fac*li(i,j,k,ismpl)/(davg*davg)
              IF (val>neumann_maxx) neumann_maxx = val
              IF (val<neumann_minx) neumann_minx = val
              neumann_avgx = neumann_avgx + val
            END DO
          END DO
        END DO
!$OMP end do nowait

! val in y
        c2 = 0
        neumann_maxy = small
        neumann_miny = big
        neumann_avgy = 0.0D0
!$OMP do schedule(static)
        DO k = 1, k0
          DO j = 2, j0 - 1
            DO i = 1, i0
              c2 = c2 + 1
              davg = 0.5D0*(dely(j)+dely(j+1))
              fac = delt/rhoceff(i,j,k,ismpl)
              val = fac*lj(i,j,k,ismpl)/(davg*davg)
              IF (val>neumann_maxy) neumann_maxy = val
              IF (val<neumann_miny) neumann_miny = val
              neumann_avgy = neumann_avgy + val
            END DO
          END DO
        END DO
!$OMP end do nowait

! val in z
        c3 = 0
        neumann_maxz = small
        neumann_minz = big
        neumann_avgz = 0.0D0
!$OMP do schedule(static)
        DO k = 2, k0 - 1
          DO j = 1, j0
            DO i = 1, i0
              c3 = c3 + 1
              davg = 0.5D0*(delz(k)+delz(k+1))
              fac = delt/rhoceff(i,j,k,ismpl)
              val = fac*lk(i,j,k,ismpl)/(davg*davg)
              IF (val>neumann_maxz) neumann_maxz = val
              IF (val<neumann_minz) neumann_minz = val
              neumann_avgz = neumann_avgz + val
            END DO
          END DO
        END DO
!$OMP end do nowait

!     compute global sum for all values
        CALL omp_summe(neumann_maxx,neumann_minx,neumann_avgx, &
          neumann_maxy,neumann_miny,neumann_avgy,neumann_maxz, &
          neumann_minz,neumann_avgz,c1,c2,c3,ismpl)

!$OMP master
        IF (i0>2) THEN
          neumann_avgx = neumann_avgx/dble(c1)
        ELSE
          neumann_maxx = 0.0D0
          neumann_minx = 0.0D0
          neumann_avgx = 0.0D0
        END IF
        IF (j0>2) THEN
          neumann_avgy = neumann_avgy/dble(c2)
        ELSE
          neumann_maxy = 0.0D0
          neumann_miny = 0.0D0
          neumann_avgy = 0.0D0
        END IF
        IF (k0>2) THEN
          neumann_avgz = neumann_avgz/dble(c3)
        ELSE
          neumann_maxz = 0.0D0
          neumann_minz = 0.0D0
          neumann_avgz = 0.0D0
        END IF

        neumann_max = max(neumann_maxx,neumann_maxy,neumann_maxz)

        IF (linfos(3)>=2) THEN
          WRITE(*,*) 'neumann number for temperature in x,y,z:'
          WRITE(*,'(a,1e10.3,a,1e10.3,a,1e10.3)') '  max. : ', &
            neumann_maxx, ',    ', neumann_maxy, ',    ', neumann_maxz
          WRITE(*,'(a,1e10.3,a,1e10.3,a,1e10.3)') '  min. : ', &
            neumann_minx, ',    ', neumann_miny, ',    ', neumann_minz
          WRITE(*,'(a,1e10.3,a,1e10.3,a,1e10.3)') '  avg. : ', &
            neumann_avgx, ',    ', neumann_avgy, ',    ', neumann_avgz
        END IF

        IF (linfos(3)>=1 .AND. neumann_max>1.D0) THEN
          WRITE(*,'(a)') &
            '!!!: neumann temp number(s) greater than 1 :'
          WRITE(*,'(a,1e12.3,a,1e10.3,a,1e10.3)') 'x: ', &
            neumann_maxx, 'y: ', neumann_maxy, 'z: ', neumann_maxz
          WRITE(*,*)
        END IF
!$OMP end master

        RETURN
      END


subroutine BALL (PES, RHA)
  use simsphere_mod
  implicit none

!        INCLUDE 'modvars.h'

!/ 24th March 1995 -- RST is always unscaled in this subroutine.

!       A BRAND NEW ROUTINE TO GET FCO2 USING PHOTO MODEL
!
!       PAR      >> PAR IN MOL UNITS
!       RRTOT    >> TOTAL RESISTANCE FOR CO2 DIFFUSION IN MOL UNITS
!       RAIR     >> SURFACE LAYER RESISTANCE
!       FCO2     >> CO2 FLUX (EITHER SCALED PER LEAF OR TOTAL)
!       CII      >> INTERNAL LEAF CO2 IN PPM

!       THE FOLLOWING ARE FROM THE TWO FARQUHAR PAPERS AND THE
!       "NOTE ON CARBON FIXATION MODELING" CONTAINS THE CURRENT
!       PARAMETERIZATIONS FOR CCOMP, KO, KC, WR.  MOST OF THE
!       PARAMETERS ARE NOW SET ACCORDING TO COLLATZ, BALL, BERRY.
!
!       AKC >> MM CONSTANT FOR CO2
!       AKO >> MM COMNSTANT FOR O2
!       OXY >> OXYGEN CONCENTRATION IN CHLOROPLASTS
!       VCMAX >> MAXIMUM CARBOXYLATION RATE
!       ALPI >> QUANTUM YEILD (PHOTONS PER CO2)
!       CCOMP >> CO2 COMPENSATION POINT
!       FCO2 >> CARBON EXCHANGE RATE
!       RESP >> DARK RESPIRATION
!       CII >> INTERNAL CO2 CONCENTRATION
!       WE >> RUBISCO LIMITED PHOTOSYNTHESIS
!       WR >> RUBP REGENERATION LIMITED PHOTOSYNTHESIS
!       VC >> CARBON EXCAHNGE IN ABSENSE OF PHOTORESPIRATION
!       PARQ >> PAR WHEN WR = WE FOR GIVEN CONDITIONS
!       TAU >> CO2/O2 SPECIFICITY RATIO
!
!       CSS >>>    LEAF SFC CO2 CONCENTRATION
!       GS  >>>    STOMATAL CONDUCTANCE

!
!       CARBON FLUX AND RESISTANCES SHOULD BE PER UNIT LEAF AREA
!       TO GET PSHYSIOLOGICAL CII
!
  character(len=7) :: LIMIT

  real :: PES, RHA, PAR, RAIR, RRTOT, AKC, TAU, AKO, VCMAX, RESP, CCOMP
  real :: AKC0, AKO0, OXY, VCMAX0, ALPI, RESP0, TAU0, AJPAR, AJMAX, AMPAR
  real :: CII, WE, AJ, WR, WS, VC, PARQ
  real :: AAA, BBB, CCC, GS, QBND, RHL, CSS

  RST = RST * XLAI / PES
  RAF = RAF * XLAI / PES  ! Unscaled
  FCO2 = FCO2 * PES / XLAI / FRVEG  ! Unscaled

!
!       PAR IN MOL UNITS
!
  PAR = ((SOL/2)*4.57)*1E-6      ! IN MOL UNITS

  RAIR = RHA + RZASCR
  RRTOT = 1.32 * RAF + 1.66*RST + RAIR
  RRTOT = RRTOT / 40.             ! IN MOL UNITS

!	SET SOME PARAMTERS  PhotoSyn Model (NOW READ IN START)

  AKC0 = 3.0E-4
  AKO0 = 0.300
  OXY = 0.209
  VCMAX0 =  2.0E-4 !COLLATZ  !7.5E-5    DATA 1992
  ALPI = 12
  RESP0 = 3.0*1E-6
  TAU0 = 2600
  AJPAR = 380E-6

  AMPAR = 9.0   ! Conductance Model

!
!       CORRECT KINEMATIC PROPERTIES TO TEMPERATURE
!	SEE COLLATZ, BALL, BERRY, FARQUHAR
!
  AKC = AKC0*(2.1**((TF-298)/10))
  TAU = TAU0*(0.57**((TF-298)/10))
  AKO = AKO0*(1.2**((TF-298)/10))
  VCMAX = VCMAX0*(2.4**((TF-298)/10))
  RESP = RESP0*(2.0**((TF-298)/10))
  CCOMP = OXY/(2*TAU)
  AJMAX = AJPAR * (1. + 0.0409*(TF-303.) - 1.54E-3 * ((TF-303.)**2) - 9.42E-5 * ((TF-303.)**3))

!
!	CORRECT RESP INHIBITION AT HIGH TEMPERATURES
!       FROMM COLLATZ, BALL, BERRY
!

  RESP = RESP*1/(1+EXP((1.3*(TF-328))))

!
!       DO INTERATIONS TO GET CII STABLE WHEN CONDITIONS CHANGING
!       RAPIDLY.
!

20      CONTINUE

!
!	GUESS AT CI
!

  CI = CO - FCO2 * RRTOT
  CII = CO - FCO2 * RRTOT  ! Previous Value

  if (CI < 0) CI = 220*1E-6

!
!       CALCULATE NEW RUBISCO AND RUBP LIMITED RATES
!

  WE = VCMAX*CI/(CI + AKC*(1+(OXY/AKO)))
  AJ = AJMAX * PAR / (PAR + 2.1 * AJMAX)
  WR = AJ * CI / (4.5 * CI + 10.5 * CCOMP)
  WS = VCMAX / 2

!
!       DETERMINE THE FACTOR WHICH IS LIMITIING
!

  IF (WE <= WR .and. WE <= WS) THEN
    LIMIT = 'RUBISCO'
  ELSE IF (WR < WE .and. WR < WS) THEN
    LIMIT = 'E TRANS'
  ELSE
    LIMIT = 'SINK'
  END IF

!
!       SET CARBOXYLATION RATE TO THE LIMITING FACTOR
!

  VC = MIN(WE,WR,WS)
  PARQ = WE * (ALPI*(1+(2*CCOMP/CI)))

!
!       CO2 FLUX
!

  FCO2 = VC*(1 - CCOMP/CI) - RESP
  IF (CCOMP > CI) FCO2 = - RESP

!
!       BALL BERRY STOMATAL CONDUCTANCE MODEL
!

!         gs = AMPAR * fco2 * (rhl)/css + 0.01  !IN MOL UNITS

!
!       ESITIMATE OF Gs ---- RECHECK THIS: THERE MAY BE A PROBLEM
!       USING 1.32 FOR WATER VAPOR.
!

  AAA = CI
  BBB = CO / (RAF * 1.32 + RAIR) - AMPAR * FCO2
  CCC = -AMPAR * FCO2 * QAF / (QSTF * RAF)

!
!       SOLUTION IF NO RELATIVE HUMIDITY DEPENDANCE
!	Following 2 statements ....

!	 BBB = CO/(RAF * 1.32 + RAIR)
!        CCC = - AMPAR * FCO2

  IF (BBB**2 < 4*AAA*CCC) THEN
    RST = RCUT
  ELSE
    GS = (-BBB + SQRT(BBB**2 - 4*AAA*CCC))/(2*AAA)
    GS = GS / 40.0    !MKS UNITS
!    IF (.not. eq(GS,0.0)) RS = 1. / GS  ! Checks
!    IF (eq(GS,0.0)) RS = 500.0     ! -"-
    if (eq(gs,0.0)) then
      rs = 500.0
    else
      rs = 1.0 / gs
    end if
    RST = RS * RCUT / ( RS + RCUT )
    RRTOT = 1.32 * RAF + 1.66 * RST + RAIR
    RRTOT = RRTOT / 40.0      ! Mole units

    CI = CO - RRTOT * FCO2

!
!	CHECK IF CI IS CONVERGING
!

    IF (ABS(1 - CI/CII) > 0.05) GOTO 20

  END IF

  IF (FCO2 < 0.0) RST = RCUT

!
!	SCALE THE RESISTANCE AND CO2 FLUX
!	Not Used in any equations -- just for output

  QBND = (QSTF/RST + QAF/RAF) * (1/RAF + 1/RST)**(-1)
! Humidity at leaf surface or leaf boundary layer
  RHL = QBND/QSTF  !  RH at Leaf Surface

  CCAN = (CI/(RST*1.66 + RAF*1.32) + CO/(RAIR))*(1/(RST*1.66 + RAF*1.32) + 1/RAIR)**(-1) ! CO2 in canopy
  CSS = FCO2 * RAF * 1.66 + CCAN    ! CO2 at leaf surface

  RAF = RAF * PES / XLAI   ! Scaling Back
  RST = RST * PES / XLAI
  FCO2 = FCO2 * XLAI / PES * FRVEG

  return
end


! ============================================================================
! "Read_Archive_Map" reads a sky map from a binary table in a FITS file.  The
! map consists of two columns in the table; the first contains temperature
! and the second contains N_Obs.
!
! Arguments:
!   FileName  Character  The FITS filename.
!   Temp      Real*4     The array of temperatures.
!   NObs      Real*4     The array of observation numbers.
!   NPix      Integer    The number of elements in the output arrays.
!   Status    Integer    A status code: 0=success.
!
! Written by Michael R. Greason, SSAI, 07 February 2006.
! ============================================================================
Subroutine Read_Archive_Map (FileName, Temp, NObs, NPix, Status)
!
Character(*),                Intent(In)  :: FileName
Integer (Kind=4),            Intent(Out) :: NPix
Integer (Kind=4),            Intent(Out) :: Status
Real (Kind=4), Dimension(:), Intent(Out) :: Temp
Real (Kind=4), Dimension(:), Intent(Out) :: NObs
!
Character (80)   :: str
Integer (Kind=4) :: i, n, Unit, blk, typ
Logical          :: anyf
! ----------------------------------------------------------------------------
!           Open the FITS file.
!
NPix   = 0
Status = 0
i      = 0
Call FTGIOU (Unit, Status)
Call FTOPEN (Unit, FileName, i, blk, Status)
!
!           Select the HDU.
!
i   = 1
typ = 0
Do While ((typ .NE. 2) .AND. (Status .EQ. 0))
  i = i + 1
  Call FTMAHD (Unit, i, typ, Status)
End Do
If ((typ .NE. 2) .OR. (Status .NE. 0)) Then
  Print 1, trim(FileName)
  Return
End If
!
!           Extract the number of pixels in the map.
!
Call FTGNRW (Unit, i, Status)
If (Status .NE. 0) Return
NPix = i
!
!           Extract the map.
!
Call FTGCVE (Unit, 1, 1, 1, NPix, 0.0E0, Temp, anyf, Status)
Call FTGCVE (Unit, 2, 1, 1, NPix, 0.0E0, NObs, anyf, Status)
!
!           Close the FITS file.
!
Call FTCLOS (Unit, Status)
If (Unit .GE. 50) Call FTFIOU (Unit, Status)
!
Return
! ----------------------------------------------------------------------------
1 Format ('Read_Archive_Map: There is no binary FITS table extension in ',A)
!
End Subroutine Read_Archive_Map


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

Howdy mysterious stranger, and Welcome to the Wiki!   You might want to check out the importance of using your RealName, just so we can get to know you.  Or not: its your choice, but people were pretty friendly here, and this is a community of local people.  My names Evan, pleased to meet you!  You know, it might be better to just edit the Renaissance Park Apartments directly to add some of that kind of stuff to the main part of the page.  You might want to check with management to confirm some of the upgrades first. Users/JabberWokky Evan JabberWokky Edwards


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine precfd(ne,ipkon,kon,lakon,ipnei,neifa,neiel,ipoface,
     &  nodface,ielfa,nkonnei,nface,ifaext,nfaext,
     &  isolidsurf,nsolidsurf,set,nset,istartset,iendset,ialset,
     &  vel,vold,mi,neij,nef,nactdoh,ipkonf,lakonf,ielmatf,ielmat,
     &  ielorienf,ielorien,norien)
!
!     gathering topological information for computational fluid
!     dynamics calculations
!
      implicit none
!
      character*8 lakon(*),lakonf(*)
      character*81 set(*),noset
!
      integer ne,ipkon(*),ipnei(*),ipoface(*),nodface(5,*),neifa(*),
     &  ielfa(4,*),nkonnei,nface,i,j,k,index,indexe,neiel(*),ithree,
     &  nfaext,ifaext(*),isolidsurf(*),nsolidsurf,indexf,
     &  nset,istartset(*),iendset(*),ialset(*),iaux,kflag,ifour,iel2,
     &  ifaceq(8,6),ifacet(7,4),ifacew(8,5),kon(*),nodes(4),iel1,j2,
     &  indexold,ifree,ifreenew,ifreenei,mi(*),neij(*),ifreenei2,nef,
     &  nactdoh(*),ipkonf(*),ielmatf(mi(3),*),ielmat(mi(3),*),nf(5),
     &  nope,ielorien(mi(3),*),ielorienf(mi(3),*),norien,jopposite8(6),
     &  jopposite6(5)
!
      real*8 vel(nef,0:5),vold(0:mi(2),*)
!
!     nodes belonging to the element faces
!
      data ifaceq /4,3,2,1,11,10,9,12,
     &            5,6,7,8,13,14,15,16,
     &            1,2,6,5,9,18,13,17,
     &            2,3,7,6,10,19,14,18,
     &            3,4,8,7,11,20,15,19,
     &            4,1,5,8,12,17,16,20/
      data ifacet /1,3,2,7,6,5,11,
     &             1,2,4,5,9,8,12,
     &             2,3,4,6,10,9,13,
     &             1,4,3,8,10,7,14/
      data ifacew /1,3,2,9,8,7,0,0,
     &             4,5,6,10,11,12,0,0,
     &             1,2,5,4,7,14,10,13,
     &             2,3,6,5,8,15,11,14,
     &             4,6,3,1,12,15,9,13/
      data jopposite6 /2,1,0,0,0/
      data jopposite8 /2,1,5,6,3,4/
      data nf /3,3,4,4,4/
!
!     transfer element information from ipkon, lakon, ielmat
!     and ielorien to ipkonf, lakonf, ielmatf and ielorienf.
!     The latter fields were created
!     for CFD applications in which the elements are numerated in
!     ascending order (without gaps)
!
      do i=1,ne
c         write(*,*) 'precfd ',i,nactdoh(i)
         if(nactdoh(i).ne.0) then
            ipkonf(nactdoh(i))=ipkon(i)
            lakonf(nactdoh(i))=lakon(i)
            do j=1,mi(3)
               ielmatf(j,nactdoh(i))=ielmat(j,i)
               if(norien.gt.0) ielorienf(j,nactdoh(i))=ielorien(j,i)
            enddo
         endif
      enddo
!
      kflag=1
      ithree=3
      ifour=4
!
!     determining the external element faces of the fluid mesh
!     the faces are catalogued by the three lowes nodes numbers
!     in ascending order. ipoface(i) points to a face for which
!     node i is the lowest node and nodface(1,ipoface(i)) and
!     nodface(2,ipoface(i)) are the next lower ones.
!     nodface(3,ipoface(i)) contains the element number,
!     nodface(4,ipoface(i)) the face number and nodface(5,ipoface(i))
!     is a pointer to the next surface for which node i is the
!     lowest node; if there are no more such surfaces the pointer
!     has the value zero
!     An external element face is one which belongs to one element
!     only
!
      ifree=1
      ifreenei=0
!
      do i=1,6*nef
         nodface(5,i)=i+1
      enddo
      do i=1,nef
         indexe=ipkonf(i)
         if(lakonf(i)(4:4).eq.'8') then
!
!           hex element
!
            ipnei(i)=ifreenei
            do j=1,6
               do k=1,4
                  nodes(k)=kon(indexe+ifaceq(k,j))
               enddo
               call isortii(nodes,iaux,ifour,kflag)
               indexold=0
               index=ipoface(nodes(1))
               do
!
!                 adding a surface which has not been
!                 catalogued so far
!
                  if(index.eq.0) then
                     ifreenew=nodface(5,ifree)
                     nodface(1,ifree)=nodes(2)
                     nodface(2,ifree)=nodes(3)
                     nodface(3,ifree)=i
                     nodface(4,ifree)=j
                     nodface(5,ifree)=ipoface(nodes(1))
                     ipoface(nodes(1))=ifree
                     ifreenei=ifreenei+1
                     neifa(ifreenei)=ifree
                     ielfa(1,ifree)=i
                     ielfa(4,ifree)=j
                     ifree=ifreenew
                     exit
                  endif
!
!                 removing a surface which has already
!                 been catalogued
!
                  if((nodface(1,index).eq.nodes(2)).and.
     &               (nodface(2,index).eq.nodes(3))) then
                     ifreenei=ifreenei+1
!
!                    completing the facial info in neifa
!
                     neifa(ifreenei)=index
                     ielfa(2,index)=i
!
!                    the neighboring elements to the face are i and iel2
!                    with local face number for the face of j and j2
!
                     iel2=ielfa(1,index)
                     j2=ielfa(4,index)
!
!                    completing the neighboring info for (element i,side j)
!
                     neiel(ifreenei)=iel2
                     neij(ifreenei)=j2
!
!                    completing the neighboring info for (element iel2,side j2)
!
                     ifreenei2=ipnei(iel2)+j2
                     neiel(ifreenei2)=i
                     neij(ifreenei2)=j
                     exit
                  endif
                  indexold=index
                  index=nodface(5,index)
               enddo
            enddo
         else if(lakonf(i)(4:4).eq.'6') then
!
!           wedge element
!
            ipnei(i)=ifreenei
            do j=1,5
               do k=1,nf(j)
                  nodes(k)=kon(indexe+ifacew(k,j))
               enddo
               call isortii(nodes,iaux,nf(j),kflag)
               indexold=0
               index=ipoface(nodes(1))
               do
!
!                 adding a surface which has not been
!                 catalogued so far
!
                  if(index.eq.0) then
                     ifreenew=nodface(5,ifree)
                     nodface(1,ifree)=nodes(2)
                     nodface(2,ifree)=nodes(3)
                     nodface(3,ifree)=i
                     nodface(4,ifree)=j
                     nodface(5,ifree)=ipoface(nodes(1))
                     ipoface(nodes(1))=ifree
                     ifreenei=ifreenei+1
                     neifa(ifreenei)=ifree
                     ielfa(1,ifree)=i
                     ielfa(4,ifree)=j
                     ifree=ifreenew
                     exit
                  endif
!
!                 removing a surface which has already
!                 been catalogued
!
                  if((nodface(1,index).eq.nodes(2)).and.
     &               (nodface(2,index).eq.nodes(3))) then
                     ifreenei=ifreenei+1
!
!                    completing the facial info in neifa
!
                     neifa(ifreenei)=index
                     ielfa(2,index)=i
!
!                    the neighboring elements to the face are i and iel2
!                    with local face number for the face of j and j2
!
                     iel2=ielfa(1,index)
                     j2=ielfa(4,index)
!
!                    completing the neighboring info for (element i,side j)
!
                     neiel(ifreenei)=iel2
                     neij(ifreenei)=j2
!
!                    completing the neighboring info for (element iel2,side j2)
!
                     ifreenei2=ipnei(iel2)+j2
                     neiel(ifreenei2)=i
                     neij(ifreenei2)=j
                     exit
                  endif
                  indexold=index
                  index=nodface(5,index)
               enddo
            enddo
         else
!
!           tet element
!
            ipnei(i)=ifreenei
            do j=1,4
               do k=1,3
                  nodes(k)=kon(indexe+ifacet(k,j))
               enddo
               call isortii(nodes,iaux,ithree,kflag)
               indexold=0
               index=ipoface(nodes(1))
               do
!
!                 adding a surface which has not been
!                 catalogued so far
!
                  if(index.eq.0) then
                     ifreenew=nodface(5,ifree)
                     nodface(1,ifree)=nodes(2)
                     nodface(2,ifree)=nodes(3)
                     nodface(3,ifree)=i
                     nodface(4,ifree)=j
                     nodface(5,ifree)=ipoface(nodes(1))
                     ipoface(nodes(1))=ifree
                     ifreenei=ifreenei+1
                     neifa(ifreenei)=ifree
                     ielfa(1,ifree)=i
                     ielfa(4,ifree)=j
                     ifree=ifreenew
                     exit
                  endif
!
!                 removing a surface which has already
!                 been catalogued
!
                  if((nodface(1,index).eq.nodes(2)).and.
     &               (nodface(2,index).eq.nodes(3))) then
                     ifreenei=ifreenei+1
!
!                    completing the facial info in neifa
!
                     neifa(ifreenei)=index
                     ielfa(2,index)=i
!
!                    the neighboring elements to the face are i and iel2
!                    with local face number for the face of j and j2
!
                     iel2=ielfa(1,index)
                     j2=ielfa(4,index)
!
!                    completing the neighboring info for (element i,side j)
!
                     neiel(ifreenei)=iel2
                     neij(ifreenei)=j2
!
!                    completing the neighboring info for (element iel2,side j2)
!
                     ifreenei2=ipnei(iel2)+j2
                     neiel(ifreenei2)=i
                     neij(ifreenei2)=j
                     exit
                  endif
                  indexold=index
                  index=nodface(5,index)
               enddo
            enddo
         endif
      enddo
!
      nkonnei=ifreenei
      nface=ifree-1
!
!     catalogueing external faces
!
      nfaext=0
!
      do i=1,nface
         if(ielfa(2,i).ne.0) cycle
         nfaext=nfaext+1
         ifaext(nfaext)=i
c         j=ielfa(4,i)
c         iel1=ielfa(1,i)
c         indexf=ipnei(iel1)
c         j=jopposite(j)
         iel1=ielfa(1,i)
         if(lakonf(iel1)(4:4).eq.'8') then
            j=ielfa(4,i)
            j=jopposite8(j)
         elseif(lakonf(iel1)(4:4).eq.'6') then
            j=ielfa(4,i)
            j=jopposite6(j)
            if(j.eq.0) cycle
         else
            cycle
         endif
         indexf=ipnei(iel1)
         ielfa(3,i)=neiel(indexf+j)
      enddo
!
c      do i=1,ne
c         k=nactdoh(i)
c         if(k.eq.0) cycle
c         write(*,*)'precfd neiel',i,k,ipnei(k),(neiel(ipnei(k)+j),j=1,6)
c         write(*,*)'precfd neij',i,k,ipnei(k),(neij(ipnei(k)+j),j=1,6)
c         write(*,*)'precfd neifa',i,k,ipnei(k),(neifa(ipnei(k)+j),j=1,6)
c      enddo
c      do k=1,nface
c         write(*,*) 'precfd ielfa',k,(ielfa(j,k),j=1,4)
c      enddo
c      do k=1,nfaext
c         write(*,*) 'precfd ifaext',k,ifaext(k)
c      enddo
c      write(*,*)
!
!     faces belonging to solid surfaces
!
      nsolidsurf=0
      noset(1:13)='SOLIDSURFACET'
      do i=1,nset
         if(set(i)(1:13).eq.noset(1:13)) exit
      enddo
      if(i.gt.nset) then
         write(*,*) '*WARNING in precfd: facial surface SOLID SURFACE '
         write(*,*) '         has not been defined.'
         write(*,*)
      else
c         nsolidsurf=0
         do j=istartset(i),iendset(i)
            nsolidsurf=nsolidsurf+1
            isolidsurf(nsolidsurf)=ialset(j)
         enddo
         call isortii(isolidsurf,iaux,nsolidsurf,kflag)
      endif
c      do i=1,nsolidsurf
c         write(*,*) 'precfd solidsurf ',i,isolidsurf(i)
c      enddo
!
!     initial conditions: element values
!     vel was initialized at allocation
!
      do i=1,nef
         indexe=ipkonf(i)
         if(lakonf(i)(4:4).eq.'8') then
            nope=8
         elseif(lakonf(i)(4:4).eq.'6') then
            nope=6
         else
            nope=4
         endif
         do j=0,4
            do k=1,nope
               vel(i,j)=vel(i,j)+vold(j,kon(indexe+k))
            enddo
            vel(i,j)=vel(i,j)/nope
         enddo
      enddo
c      do i=1,nef
c         write(*,*) 'precfd vel ',i,(vel(j,i),j=1,3)
c      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
A favorite amongst westside Arboretum frequenters, Monsieur Tig McJRoolems resides within a largeish bush beside the UCD Equestrian Center.  Described as dignified, enigmatic, and a relentless cuddlywhumper, Tig takes to young female strangers, enticing them with his soulful eyes, lustrous coat, and feline charm.  But unlike any other cat you might meet, Tig has a rather captivating secret...  If you take the time to acquaint yourself with himgiven his unpretentious presenceand treat him well, he might just bestow upon you, a startlingly affectionate gift...

Tig McJRoolems was recently adopted by a wonderful family to try and give him a warm loving home. Upon bringing him into the veterinarian for a standard health check it was discovered that he was positive for FIV (feline HIV). FIV is a feline only disease and is contagious to any unvaccinated cats. It is transferred just like the human strain of HIV/AIDS and is just as fatal. Unfortunately Tig had to be humanely euthanized due to the suffering the disease causes. He will be greatly missed. For more information on FIV please visit http://www.vet.cornell.edu/fhc/brochures/fiv.html


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      LOGICAL FUNCTION LXSCNP(DELIM,STR,NS,CH)
      IMPLICIT INTEGER (A-Z)
      CHARACTER*504 ILINE
      COMMON /LXCOM1/ILINE
      COMMON /LXCOM2/JLINE,LXINIT
      CHARACTER*(*) DELIM
      CHARACTER*(*) STR
      INTEGER NS
      CHARACTER CH
      LOGICAL QFLAG
      CHARACTER QCH

      NS = 0
      J = JLINE
      LP = 0
      QFLAG = .FALSE.
 2600 CONTINUE
      CH = ILINE(J:J)
      IF (LP.EQ.0) THEN
         ID = INDEX(DELIM,CH)
         IF (ID.GT.0) THEN
            GO TO 2620

         END IF

      END IF

      IF (.NOT.QFLAG) THEN
         IF (LP.EQ.0) THEN
            ID = INDEX(DELIM,CH)
            IF (ID.GT.0) THEN
               GO TO 2620

            END IF

         END IF

         IF (CH.EQ.CHAR(0)) THEN
            GO TO 2620

         ELSE IF (CH.EQ.'(' .OR. CH.EQ.'[' .OR. CH.EQ.'{') THEN
            LP = LP + 1

         ELSE IF (CH.EQ.')' .OR. CH.EQ.']' .OR. CH.EQ.'}') THEN
            LP = LP - 1

         ELSE IF (CH.EQ.'''' .OR. CH.EQ.'"' .OR. CH.EQ.CHAR(96)) THEN
            QFLAG = .TRUE.
            QCH = CH
         END IF

         IF (LP.LT.0) THEN
            GO TO 2620

         END IF

         NS = NS + 1
         STR(NS:NS) = CH

      ELSE
         NS = NS + 1
         STR(NS:NS) = CH
         IF (CH.EQ.CHAR(0)) THEN
            GO TO 2620

         ELSE IF (CH.EQ.QCH) THEN
            IF (CH.EQ.ILINE(J+1:J+1)) THEN
               NS = NS + 1
               STR(NS:NS) = CH
               J = J + 2
               GO TO 2610

            ELSE
               QFLAG = .FALSE.
            END IF

         END IF

      END IF

      J = J + 1
 2610 GO TO 2600

 2620 CONTINUE
      JLINE = J
      LXSCNP = ((LP.EQ.0) .AND. .NOT.QFLAG)
      RETURN

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
!
!     (c) 2019 Guide Star Engineering, LLC
!     This Software was developed for the US Nuclear Regulatory Commission (US NRC)
!     under contract "Multi-Dimensional Physics Implementation into Fuel Analysis under
!     Steady-state and Transients (FAST)", contract # NRC-HQ-60-17-C-0007
!
!
!    NEMO - Numerical Engine (for) Multiphysics Operators
! Copyright (c) 2007, Stefano Toninel
!                     Gian Marco Bianchi  University of Bologna
!              David P. Schmidt    University of Massachusetts - Amherst
!              Salvatore Filippone University of Rome Tor Vergata
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification,
! are permitted provided that the following conditions are met:
!
!     1. Redistributions of source code must retain the above copyright notice,
!        this list of conditions and the following disclaimer.
!     2. Redistributions in binary form must reproduce the above copyright notice,
!        this list of conditions and the following disclaimer in the documentation
!        and/or other materials provided with the distribution.
!     3. Neither the name of the NEMO project nor the names of its contributors
!        may be used to endorse or promote products derived from this software
!        without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
! ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
! WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
! ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!---------------------------------------------------------------------------------
!
! $Id: rd_inp_bc_math.f90 3093 2008-04-22 14:51:09Z sfilippo $
!
! Description:
!    To be added...
!
SUBMODULE(tools_bc) rd_inp_bc_math_implementation
    IMPLICIT NONE

    CONTAINS

        MODULE PROCEDURE rd_inp_bc_math
        USE class_psblas
        USE json_module, ONLY : json_file
        USE tools_bc, ONLY: bc_dirichlet_, bc_neumann_, bc_robin_, &
            & bc_dirichlet_map_, bc_neumann_map_, bc_robin_map_
        USE tools_input

        IMPLICIT NONE
        !
        LOGICAL, PARAMETER :: debug = .FALSE.
        !
        TYPE(json_file) :: nemo_json
        INTEGER :: mypnum, icontxt
        INTEGER :: info, inp, na, nb, nc
        CHARACTER(len=15) :: par
        REAL(psb_dpk_) :: abc(3)

        icontxt = icontxt_()
        mypnum  = mypnum_()

        IF(mypnum == 0) THEN
            CALL open_file(input_file,nemo_json)
            CALL find_section(sec,nemo_json)

            WRITE(*,*) '- Reading ',TRIM(sec),' section: type MATH'

            READ(inp,'()')

            ! Reads BC section and the corresponding ID
            READ(inp,*) par, id

            ! Broadcasts the ID
            CALL psb_bcast(icontxt,id)

            ! Reads ABC parameters
            SELECT CASE(id)
            CASE(bc_dirichlet_)
                abc(1) = 1.d0
                abc(2) = 0.d0
                READ(inp,*) par, abc(3)
            CASE(bc_neumann_)
                abc(1) = 0.d0
                abc(2) = 1.d0
                READ(inp,*) par, abc(3)
            CASE(bc_robin_)
                READ(inp,*) par, abc(1), abc(2), abc(3)
            CASE DEFAULT
                abc(:) = 0.d0
            END SELECT

            ! Broadcasts ABC parameters (only uniform BCs)
            CALL psb_bcast(icontxt,abc)

            CLOSE(inp)
        ELSE
            CALL psb_bcast(icontxt,id)
            CALL psb_bcast(icontxt,abc)
        END IF

        ! Allocates A, B, C pointers depending on the BC type
        SELECT CASE(id)
        CASE(  bc_dirichlet_, &
            & bc_neumann_,   &
            & bc_robin_ )
            na = 1; nb = 1; nc = 1
        CASE(  bc_dirichlet_map_, &
            & bc_neumann_map_)
            na = 1; nb = 1; nc = nbf
        CASE(  bc_robin_map_)
            na = nbf; nb = nbf; nc = nbf
        CASE DEFAULT
            WRITE(*,100)
            CALL abort_psblas
        END SELECT
        ALLOCATE(a(na),b(nb),c(nc),stat=info)
        IF(info /= 0) THEN
            WRITE(*,200)
            CALL abort_psblas
        END IF

        ! Sets BC values
        SELECT CASE(id)
        CASE(bc_dirichlet_, bc_neumann_, bc_robin_)
            a(1) = abc(1)
            b(1) = abc(2)
            c(1) = abc(3)
        CASE(bc_dirichlet_map_)
            a(1) = 1.d0
            b(1) = 0.d0
            c(:) = 0.d0 ! To be actually set at the first mapping
        CASE(bc_neumann_map_)
            a(1) = 0.d0
            b(1) = 1.d0
            c(:) = 0.d0 ! To be actually set at the first mapping
        CASE(bc_robin_map_)
            a(:) = 0.d0 ! To be actually set at the first mapping
            b(:) = 0.d0 ! To be actually set at the first mapping
            c(:) = 0.d0 ! To be actually set at the first mapping
        CASE DEFAULT
            WRITE(*,100)
            CALL abort_psblas
        END SELECT

        ! Debug
        IF(debug) THEN
            WRITE(*,*)
            WRITE(*,300) mypnum
            WRITE(*,400) TRIM(sec),' - Type: Math'
            WRITE(*,500) '  BC%id     = ', id
            WRITE(*,600) '  BC%a = ', a
            WRITE(*,600) '  BC%b = ', b
            WRITE(*,600) '  BC%c = ', c
            WRITE(*,*)
        END IF

100     FORMAT(' ERROR! Unsupported BC in RD_INP_BC_MATH')
200     FORMAT(' ERROR! Memory allocation failure in RD_INP_BC_MATH')
300     FORMAT(' ----- Process ID = ',i2,' -----')
400     FORMAT(1x,a,a)
500     FORMAT(1x,a,i2)
600     FORMAT(1x,a,es10.3)

        END PROCEDURE rd_inp_bc_math

END SUBMODULE rd_inp_bc_math_implementation

max_stars_repo_name: sourceryinstitute/MORFEUS-Source

max_issues_repo_name: sourceryinstitute/MORFEUS-Source

max_forks_repo_name: sourceryinstitute/MORFEUS-Source


text:
Image(personpicture.jpg, thumbnail)  He always wore a hat.  You know the one.

Name: George Seamer

Phone Number: No, no.

Office: Borders

Personality: Funny, selfloathing, creative, impossible to contact, intelligent, selfabsorbed. Roomie of another former California Aggie personality, Doug Rutsch.


text:
Net Impact is a local chapter within the Graduate School of Management of a nationwide association of business students working to promote socially responsible business.  The Davis Graduate School of Management has a very active chapter that frequently does site visits to companies such as clif bar, PRIDE industries and Cisco systems.  The chapter is also involved in promoting organic farming, microfinance and environmental protection.

Each quarter Net Impact holds green bag lunch seminars where a guest is invisted to speak on a subject relevant to socially responsible business.   Last years events included:
  Sarah Null, PhD candidate at Davis and Tom Philip, Pulitzer Prize winning columnist for the Sacramento Bee speaking on restoring the Hetch Hetchy Reservoir.
  Josh Dorf, President & CEO of JOG Distribution, Inc. (Producer of Organic Flower)
  UCD Professor Brad Barber, speaking on Socially Responsible Investing


text:
module cc_mg_tower_smoother_module

  use multifab_module
  use cc_stencil_module
  use mg_tower_module
  use bl_timer_module

  implicit none

contains

  subroutine cc_mg_tower_smoother(mgt, lev, ss, uu, ff, mm)

    use omp_module
    use bl_prof_module
    use cc_smoothers_module, only: gs_line_solve_1d, gs_rb_smoother_1d, &
         jac_smoother_2d, jac_smoother_3d, &
         gs_rb_smoother_2d,  gs_rb_smoother_3d, &
         fourth_order_smoother_2d, fourth_order_smoother_3d
    use itsol_module, only: itsol_bicgstab_solve

    integer        , intent(in   ) :: lev
    type( mg_tower), intent(inout) :: mgt
    type( multifab), intent(inout) :: uu
    type( multifab), intent(in   ) :: ff
    type( multifab), intent(in   ) :: ss
    type(imultifab), intent(in   ) :: mm

    real(kind=dp_t), pointer :: fp(:,:,:,:)
    real(kind=dp_t), pointer :: up(:,:,:,:)
    real(kind=dp_t), pointer :: sp(:,:,:,:)
    integer        , pointer :: mp(:,:,:,:)
    integer :: i, n, ng, nn, stat, npts
    integer :: lo(mgt%dim)
    type(bl_prof_timer), save :: bpt
    logical :: pmask(mgt%dim), singular_test
    real(kind=dp_t) :: local_eps

    call bl_proffortfuncstart("cc_mg_tower_smoother")

    singular_test = ( mgt%bottom_singular .and. mgt%coeffs_sum_to_zero )

    pmask = get_pmask(get_layout(uu))
    !
    ! Make sure to define this here so we don't assume a certain number of ghost cells for uu.
    !
    ng = nghost(uu)

    call build(bpt, "cc_mgt_smoother")

    if ( mgt%skewed_not_set(lev) ) then
       !$OMP PARALLEL DO PRIVATE(i,mp)
       do i = 1, nfabs(mm)
          mp => dataptr(mm, i)
          mgt%skewed(lev,i) = skewed_q(mp)
       end do
       !$OMP END PARALLEL DO
       mgt%skewed_not_set(lev) = .false.
    end if

    if ( (mgt%dim .eq. 1) .and. (nboxes(uu%la) .eq. 1) ) then

       call fill_boundary(uu, cross = mgt%lcross)
       !
       ! We do these line solves as a preconditioner.
       !
       do i = 1, nfabs(ff)
          up => dataptr(uu, i)
          fp => dataptr(ff, i)
          sp => dataptr(ss, i)
          mp => dataptr(mm, i)
          lo =  lwb(get_box(ss, i))
          do n = 1, mgt%nc
             call gs_line_solve_1d(sp(:,:,1,1), up(:,1,1,n), fp(:,1,1,n), &
                  mp(:,1,1,1), lo, ng)
          end do
       end do

       call fill_boundary(uu, cross = mgt%lcross)

       local_eps = 0.001
       npts = multifab_volume(ff)
       call itsol_bicgstab_solve(ss, uu, ff, mm, &
            local_eps, npts, mgt%cg_verbose, &
            mgt%stencil_type, mgt%lcross, &
            stat = stat, &
            singular_in = singular_test, &
            uniform_dh = mgt%uniform_dh)

       call fill_boundary(uu, cross = mgt%lcross)

       if ( stat /= 0 ) then
          if ( parallel_IOProcessor() ) call bl_error("BREAKDOWN in 1d CG solve")
       end if

    else
       !
       ! Cell-centered stencils.
       !
       select case ( mgt%smoother )

       case ( MG_SMOOTHER_GS_RB )

          do nn = 0, 1
             call fill_boundary(uu, cross = mgt%lcross)
             do i = 1, nfabs(ff)
                up => dataptr(uu, i)
                fp => dataptr(ff, i)
                sp => dataptr(ss, i)
                mp => dataptr(mm, i)
                lo =  lwb(get_box(ss, i))
                do n = 1, mgt%nc
                   select case ( mgt%dim)
                   case (1)
                      call gs_rb_smoother_1d(sp(:,:,1,1), up(:,1,1,n), &
                           fp(:,1,1,n), mp(:,1,1,1), lo, ng, nn, &
                           mgt%skewed(lev,i))
                   case (2)
                      call gs_rb_smoother_2d(sp(:,:,:,1), up(:,:,1,n), &
                           fp(:,:,1,n), mp(:,:,1,1), lo, ng, nn, &
                           mgt%skewed(lev,i))
                   case (3)
                      call gs_rb_smoother_3d(sp(:,:,:,:), up(:,:,:,n), &
                           fp(:,:,:,n), mp(:,:,:,1), lo, ng, nn, &
                           mgt%skewed(lev,i))
                   end select
                end do
             end do
          end do

       case ( MG_SMOOTHER_EFF_RB )

          call fill_boundary(uu, cross = mgt%lcross)
          do i = 1, nfabs(ff)
             up => dataptr(uu, i)
             fp => dataptr(ff, i)
             sp => dataptr(ss, i)
             mp => dataptr(mm, i)
             lo =  lwb(get_box(ss, i))
             do n = 1, mgt%nc
                select case ( mgt%dim)
                case (1)
                   call gs_rb_smoother_1d(sp(:,:,1,1), up(:,1,1,n), &
                        fp(:,1,1,n), mp(:,1,1,1), lo, ng, 0, &
                        mgt%skewed(lev,i))
                   call gs_rb_smoother_1d(sp(:,:,1,1), up(:,1,1,n), &
                        fp(:,1,1,n), mp(:,1,1,1), lo, ng, 1, &
                        mgt%skewed(lev,i))
                case (2)
                   call gs_rb_smoother_2d(sp(:,:,:,1), up(:,:,1,n), &
                        fp(:,:,1,n), mp(:,:,1,1), lo, ng, 0, &
                        mgt%skewed(lev,i))
                   call gs_rb_smoother_2d(sp(:,:,:,1), up(:,:,1,n), &
                        fp(:,:,1,n), mp(:,:,1,1), lo, ng, 1, &
                        mgt%skewed(lev,i))
                case (3)
                   call gs_rb_smoother_3d(sp(:,:,:,:), up(:,:,:,n), &
                        fp(:,:,:,n), mp(:,:,:,1), lo, ng, 0, &
                        mgt%skewed(lev,i))
                   call gs_rb_smoother_3d(sp(:,:,:,:), up(:,:,:,n), &
                        fp(:,:,:,n), mp(:,:,:,1), lo, ng, 1, &
                        mgt%skewed(lev,i))
                end select
             end do
          end do

       case ( MG_SMOOTHER_MINION_CROSS )

          do nn = 0, 1
             if ( (nn == 1) .and. (mgt%dim == 3) ) &
                  !
                  ! Only the 2D fourth order stencil does red-black GS sweeps.
                  !
                  exit
             call fill_boundary(uu, cross = mgt%lcross)
             do i = 1, nfabs(ff)
                up => dataptr(uu, i)
                fp => dataptr(ff, i)
                sp => dataptr(ss, i)
                mp => dataptr(mm, i)
                lo =  lwb(get_box(ss, i))
                do n = 1, mgt%nc
                   select case (mgt%dim)
                   case (2)
                      call fourth_order_smoother_2d(sp(:,:,:,1), up(:,:,1,1), &
                           fp(:,:,1,1), lo, ng, mgt%stencil_type, nn)
                   case (3)
                      call fourth_order_smoother_3d(sp(:,:,:,:), up(:,:,:,1), &
                           fp(:,:,:,1), lo, ng, mgt%stencil_type)
                   end select
                end do
             end do
          end do

       case ( MG_SMOOTHER_MINION_FULL )

          call fill_boundary(uu, cross = mgt%lcross)
          do i = 1, nfabs(ff)
             up => dataptr(uu, i)
             fp => dataptr(ff, i)
             sp => dataptr(ss, i)
             mp => dataptr(mm, i)
             lo =  lwb(get_box(ss, i))
             do n = 1, mgt%nc
                select case (mgt%dim)
                case (2)
                   call fourth_order_smoother_2d(sp(:,:,:,1), up(:,:,1,1), &
                        fp(:,:,1,1), lo, ng, mgt%stencil_type, 0)
                case (3)
                   call fourth_order_smoother_3d(sp(:,:,:,:), up(:,:,:,1), &
                        fp(:,:,:,1), lo, ng, mgt%stencil_type)
                end select
             end do
          end do

       case ( MG_SMOOTHER_JACOBI )

          call fill_boundary(uu, cross = mgt%lcross)
          do i = 1, nfabs(ff)
             up => dataptr(uu, i)
             fp => dataptr(ff, i)
             sp => dataptr(ss, i)
             mp => dataptr(mm, i)
             lo =  lwb(get_box(ss, i))
             do n = 1, mgt%nc
                select case ( mgt%dim)
                case (2)
                   call jac_smoother_2d(sp(:,:,:,1), up(:,:,1,n), fp(:,:,1,n), &
                        mp(:,:,1,1), ng)
                case (3)
                   call jac_smoother_3d(sp(:,:,:,:), up(:,:,:,n), fp(:,:,:,n), &
                        mp(:,:,:,1), ng)
                end select
             end do
          end do

       case default
          call bl_error("MG_TOWER_SMOOTHER: no such smoother")
       end select

    end if ! if (mgt%dim > 1)

    call destroy(bpt)
    call bl_proffortfuncstop("cc_mg_tower_smoother")

  end subroutine cc_mg_tower_smoother

  subroutine mg_jacobi_smoother(mgt, lev, ss, uu, ff, mm)

    use bl_prof_module
    use cc_smoothers_module, only: jac_smoother_1d, jac_smoother_2d, jac_smoother_3d

    type(mg_tower), intent(inout) :: mgt
    type(multifab), intent(inout) :: uu
    type(multifab), intent(in) :: ff
    type(multifab), intent(in) :: ss
    type(imultifab), intent(in) :: mm
    integer, intent(in) :: lev
    real(kind=dp_t), pointer :: fp(:,:,:,:)
    real(kind=dp_t), pointer :: up(:,:,:,:)
    real(kind=dp_t), pointer :: sp(:,:,:,:)
    integer        , pointer :: mp(:,:,:,:)
    integer :: i, n, ng, iter
    type(bl_prof_timer), save :: bpt

    call bl_proffortfuncstart("mg_jacobi_smoother")
    call build(bpt, "mgt_jacobi_smoother")

    ng = nghost(uu)

    if ( mgt%skewed_not_set(lev) ) then
       do i = 1, nfabs(mm)
          mp => dataptr(mm, i)
          mgt%skewed(lev,i) = skewed_q(mp)
       end do
       mgt%skewed_not_set(lev) = .false.
    end if

    if ( .not. cell_centered_q(uu) ) then
       call bl_error('MG_JACOBI_SMOOTHER ONLY DESIGNED FOR CELL-CENTERED RIGHT NOW')
    end if

    do iter = 1, mgt%nu1

       call fill_boundary(uu, cross = mgt%lcross)

       do i = 1, nfabs(ff)
          up => dataptr(uu, i)
          fp => dataptr(ff, i)
          sp => dataptr(ss, i)
          mp => dataptr(mm, i)
          do n = 1, mgt%nc
             select case ( mgt%dim)
             case (1)
                call jac_smoother_1d(sp(:,:,1,1), up(:,1,1,n), fp(:,1,1,n), &
                                     mp(:,1,1,1), ng)
             case (2)
                call jac_smoother_2d(sp(:,:,:,1), up(:,:,1,n), fp(:,:,1,n), &
                                     mp(:,:,1,1), ng)
             case (3)
                call jac_smoother_3d(sp(:,:,:,:), up(:,:,:,n), fp(:,:,:,n), &
                                     mp(:,:,:,1), ng)
             end select
          end do
       end do
    end do

    call destroy(bpt)
    call bl_proffortfuncstop("mg_jacobi_smoother")

  end subroutine mg_jacobi_smoother

end module cc_mg_tower_smoother_module

max_stars_repo_name: memmett/BoxLib

max_issues_repo_name: memmett/BoxLib

max_forks_repo_name: memmett/BoxLib


text:
REAL FUNCTION sinc(x)
IMPLICIT NONE
  REAL, INTENT(IN) :: x
  REAL, PARAMETER :: EPSILON = 1.0E-30
  IF (ABS(x) > EPSILON) THEN
    sinc = SIN(x) / x
  ELSE
    sinc = 1
  END IF
END FUNCTION sinc

PROGRAM test_sinc
IMPLICIT NONE
  REAL :: x, sinc
  INTEGER :: i, num_results, offset, adj

  WRITE(*, *) "We will calculate sinc over a range of positive and negative values."
  WRITE(*, *) "How many results would you like to see?"
  READ(*, *) num_results

  offset = num_results / 2
  IF (MOD(num_results, 2) .eq. 0) THEN
    adj = 0
  ELSE
    adj = 1
  END IF

  DO i = 1, num_results
    WRITE(*, '(F9.6)') sinc(REAL(i - adj - offset))
  END DO
  !WRITE(*, 100) "sinc(", x, ") = ", sinc(x)
  !100 FORMAT (A, F12.6, A, F12.6)
END PROGRAM test_sinc

max_stars_repo_name: evanmacbride/fortran-practice

max_issues_repo_name: evanmacbride/fortran-practice

max_forks_repo_name: evanmacbride/fortran-practice


text:
Jonathon Leathers attended UC Davis from 2003 until graduating in 2007. He became President of the Davis College Green Party in his first year at UC Davis and remained President until he graduated. Off campus he served as the Chair of the County Council of the Green Party of Yolo County while living in Davis.

Jonathon Leathers unsuccessfully ran for ASUCD ASUCD Senate Senate twice. His first campaign was in the Fall 2004 ASUCD Election on the SOSSS slate where he fell 42 votes short of being elected. He ran again, this time as an Independent independent, in the Winter 2005 ASUCD Election. Following the Winter election he was appointed as the Chair of the ASUCD Elections Committee and held that position for two consecutive years: 20052006 and 20062007. Jonathon was a supporter of Choice Voting in ASUCD elections and even wrote the Countback Amendment to utilize election data for filing Senate vacancies. During his time in ASUCD, Jonathon was involved in a number of controversial situations. These are just a few:
  Fall 2005 ASUCD Election/Unqualified Candidates Scandal
  Fall 2005 ASUCD Election/Turnout Controversy
  Election Statement Controversy

In 2007, Jonathon Leathers left Davis and moved to Portland, Oregon. He currently lives in Vancouver, B.C.

 More information
  Users/JonathonLeathers Jons personal page
  http://www.jtleathers.com/asucd/elections.html Jons ASUCD Election History Website


program main
    implicit none
    real(8), parameter :: pi = acos(-1.d0)
    integer :: clock, n, i
    integer, allocatable :: seed(:)

    integer :: n_sample = 10000000, n_accept = 0
    real(8), parameter :: Xmax = 5.d0
    real(8) :: x = 0.d0, r1, r2, A, Integral = 0.d0

    integer, allocatable :: bin(:)
    integer :: Nbin = 1000
    real(8) :: binwidth

    real(8) :: mu, theta

    !real(8) ::

    ! seed initialization
    call SYSTEM_CLOCK(clock)
    call RANDOM_SEED(size = n)
    allocate(seed(n))
    do i = 1, n
        seed(i) = clock + 37 * i
    end do
    call RANDOM_SEED(put = seed)

    ! bin initialization
    binwidth = 2.d0 * Xmax / dble(Nbin)
    allocate(bin(Nbin))
    bin = 0

    ! Metropolis
    do i = 1, n_sample
        call RANDOM_NUMBER(r1)
        r1 = (r1 - .5d0) * (2.d0 * Xmax)
        A = exp(-(r1**2 - x**2) / 2.d0)
        call RANDOM_NUMBER(r2)
        if (r2 <= A) then
            x = r1
            bin(floor((x + Xmax) / binwidth)) = bin(floor((x + Xmax) / binwidth)) + 1
            n_accept = n_accept + 1
        end if
        integral = integral + x**2
    end do
    integral = sqrt(2.d0 * pi) * integral / dble(n_sample)
    write(*,'("Metropolis : ",f10.5," Acceptance rate : ",f10.5)') integral, dble(n_accept) / dble(n_sample)
    open(unit = 1, file = 'samples.txt', status = 'unknown')
    do i = 1, Nbin
        write(1,'(f10.5,i8)') -Xmax + (i - .5d0) * binwidth, bin(i)
    end do
    close(1)

    ! Importance sampling - inverse transform
    integral = 0.d0
    do i = 1, n_sample
        call RANDOM_NUMBER(mu)
        call RANDOM_NUMBER(theta)
        theta = 2.d0 * pi * theta
        x = sqrt(-2.d0 * log(mu)) * cos(theta)
        integral = integral + x**2
    end do
    integral = sqrt(2 * pi) * integral / n_sample
    write(*,'("Importance sampling (inverse transform) : ",f10.5)') integral

    ! Importance sampling - acceptance & rejection
    integral = 0.d0
    n_accept = 0
    do i = 1,n_sample
        call RANDOM_NUMBER(x)
        x = (x - .5d0) * (2 * Xmax)
        call WEIGHT(x,r1)
        call random_NUMBER(r2)
        if (r2 <= r1) then
            integral = integral + x**2
            n_accept = n_accept + 1
        end if
    end do
    integral = sqrt(2.d0 * pi) * integral / dble(n_accept)
    write(*,'("Important sampling (acceptance & rejection) : ",f10.5," Acceptance rate : ",f10.5)')&
        integral, dble(n_accept) / dble(n_sample)
end program main

subroutine WEIGHT(x,p)
    implicit none
    real(8), parameter :: pi = acos(-1.d0)
    real(8), intent(in) :: x
    real(8), intent(out) :: p

    p = exp(-x**2 / 2.d0) / sqrt(2.d0 * pi)
end subroutine WEIGHT


text:
C    @(#)chglod.f	20.5 11/12/98
      subroutine chglod (nwrd, word, pctp, pctq, pltoto, qltoto,
     1                   pltotn, qltotn, typout, error)

      integer nwrd, error
      character word (*)*(*), typout * 1
      real pctp, pltotn, pctq, qltotn, pltoto, qltoto

C     process %LOAD_CHANGE within /CHANGE_PARAMETERS
C
C     Input: NWRD - counter of WORD().
C            WORD - character array from SCAN of BUF.
C
C     Output: NWRD - updated counter.
C             WORD - stripped character array.
C             PCTP - percent P specified.
C             PCTQ - percent Q specified.
C             PLTOTO - total quantity in MW before percentage change
C                     applied.
C             QLTOTO - total quantity in MVAR before percentage change
C                     applied.
C             PLTOTN - total quantity in MW afterpercentage change
C                      applied.
C             QLTOTN - total quantity in MVAR after percentage change
C                      applied.
C             TYPOUT - 'P' or 'Q', the disired plotted output quantity.
C             ERROR - return status (0=normal/1=error)

      include 'ipfinc/parametr.inc'

      include 'ipfinc/alpha.inc'
      include 'ipfinc/alpha2.inc'
      include 'ipfinc/arcntl.inc'
      include 'ipfinc/area.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/cbus.inc'
      include 'ipfinc/ecvar.inc'
      include 'ipfinc/lfiles.inc'
      include 'ipfinc/ordsta.inc'
      include 'ipfinc/prt.inc'
      include 'ipfinc/slnopt.inc'
      include 'ipfinc/tbx.inc'

      common /amtrx/ ikk(5,MAXBUS), ownlis(100)
      character ownlis*3

      common /is_batch / is_batch

      double precision  dpltot, dqltot

      dimension pctchg(4), ipctch(2)
      character cbown * 3, cbtyp * 1, cbkyr * 2, owntmp * 3,  own * 3
      logical found

      error = 0
      pctp = 0.0
      pctq = 0.0
      pltotn = 0.0
      qltotn = 0.0
      pltoto = 0.0
      qltoto = 0.0
      typout = ' '

C     > CHANGE_PARAMETERS, BUS = <name> <base>, QY=? -
C                         %LOAD_CHANGE, %PX = <##>, %QY = <##>, -
C                                               AREAS = NORTHWEST, -
C                                               ZONES = NB, -
C                                               OWNERS = PSP, -
C                                               TYPES = N <
C
C     Initialize IKK array:
C
C     (1,*) = 0 :  normal.
C             1 :  Bus has its load modified.
C     (2,*) = 0 :  bus is not eligible for systematic % changes.
C             1 :  bus is eligible for systematic % changes.
C     (3,*) = I    (cross index to TBX array)
C     (4,*) = 0    (not used)
C     (5,*) = 0    (not used)
C
      do nb = 1, ntot
         ikk(1,nb) = 0
         ikk(2,nb) = 0
         ikk(3,nb) = 0
         ikk(4,nb) = 0
         ikk(5,nb) = 0
      enddo

      do i = 1, ntotb
         ltyp = tbx(1,i)
         if (ltyp .lt. 10) then
            nb = tbx(2,i)
            if (ordtbx .eq. 2) nb = opt2inp(nb)
            ikk(3,nb) = i
         endif
      enddo
C                                            Search for %LOAD.
      do i = 1, nwrd
         if (word(i)(1:1) .eq. '%') then
            if (word(i)(1:2) .eq. '%P') then
               next = i + 1
               last = index(word(next), '%') - 1
               if (last .le. 0) last = lastch (word(next))
               pctp = ftn_atof (word(next)(1:last))
               if (typout .eq. ' ') typout = 'P'
               word(next) = ' '
            else if (word(i)(1:2) .eq. '%Q') then
               next = i + 1
               last = index(word(next), '%') - 1
               if (last .le. 0) last = lastch (word(next))
               pctq = ftn_atof (word(next)(1:last))
               if (typout .eq. ' ') typout = 'Q'
               word(next) = ' '
            endif
            word(i) = ' '
         endif

      enddo

      if (pctp .eq. 0.0 .and. pctq .eq. 0.0) then

         write (errbuf(1),190)
  190    format ('No %P or %Q percentages specified.')
         call prterx ('W', 1)
         error = 1
         go to 900

      endif

C     Search for OWNERS = <owner_1,...>

      kntown = 0

      do i = 1, nwrd
         if (word(i)(1:4) .eq. 'OWNE') then
            word(i) = ' '
            last = i
            do j = i+1, nwrd
               if (word(j)(1:4) .eq. 'AREA') then
                  go to 220
               else if (word(j)(1:4) .eq. 'ZONE') then
                  go to 220
               else
                  kntown = kntown + 1
                  ownlis(kntown) = word(j)
                  word(j) = ' '
                  last = j
               endif
            enddo
         endif
  210    continue
      enddo
  220 continue
C
C     Search for AREAS = <area_1,...> or
C                ZONES = <zone_1,...>
C
      do i = 1, nwrd
         if (word(i)(1:4) .eq. 'AREA') then
            word(i) = ' '
            do nb = 1, ntot
               ikk(2,nb) = 0
            enddo
            do j = i+1, nwrd
               do k = 1, ntotc
                  if (arcnam(k) .eq. word(j)) then
                     do nb = 1, ntot
                        if (jarzn(nb) .eq. k) ikk(2,nb) = 1
                     enddo
                     go to 270
                  endif
               enddo
               last = lastch (word(j))
               write (errbuf(1),260) word(j)(1:last)
  260          format ('Interchange area (', a,
     1                 ') is not in system')
               call prterx ('W', 1)

  270          word(j) = ' '
            enddo
            go to 340
         else if (word(i)(1:4) .eq. 'ZONE') then
            word(i) = ' '
            do nb = 1, ntot
               ikk(2,nb) = 0
            enddo
            do j = i+1, nwrd
               found = .false.
               do nb = 1, ntot
                  if (zone(nb) .eq. word(j)) then
                     found = .true.
                     ikk(2,nb) = 1
                  endif
               enddo
               if (.not.found) then
                  last = lastch (word(j))
                  write (errbuf(1), 300) word(j)(1:last)
  300             format ('Zone (', a, ') is not in system')
                  call prterx ('W', 1)
               endif
  310          word(j) = ' '
            enddo
            go to 340
         endif
      enddo
C
C     Neither ZONES nor AREAS have been specified.  The entire
C     system is eligible.
C
      do i = 1, ntot
         ikk(2,i) = 1
      enddo

  340 continue
C                                    Sort OWNLIS
      if (kntown .eq. 0) then
         kntown = 1
         ownlis(kntown) = '###'
      else
         do i = 1, kntown-1
            do j = i+1, kntown
               komp = kompr(ownlis(i), ownlis(i), komp)
               if (komp .gt. 0) then
                  owntmp = ownlis(i)
                  ownlis(j) = ownlis(i)
                  ownlis(i) = owntmp
               endif
            enddo
         enddo

      endif
C                                         Print out results
      pctppu = 0.01 * pctp
      pctqpu = 0.01 * pctq

      write (outbuf, 370) pctp, pctq
  370 format ('0 %LOAD_CHANGE, PL = ', f6.1, '%, QL = ', f6.1, '%')
      call prtout (1)

      if (idswb .ne. 0) then
         write (outbuf, 380)
  380    format ('0            BUS             CARD ZONE OWNER   ---- ',
     &    'ORIGINAL QUANTITIES -----  -------- FINAL QUANTITIES ----')
         call prtout (1)
         write (outbuf, 390)
  390    format ('                                               ---- ',
     1   'LOAD ---  - GENERATION --   ---- LOAD ----  - GENERATION -')
         call prtout (1)
         write (outbuf, 400)
  400    format ('                                                ',
     &   '(MW) (MVAR)   (MW)   (MVAR)    (MW)   (MVAR)   (MW)   (MVAR)')
         call prtout (1)
      endif

      ipctch(1) = 0
      ipctch(2) = 0
      do i = 1, 4
         pctchg(i) = 0.0
      enddo
      pltoto = 0.0
      qltoto = 0.0
      pltotn = 0.0
      qltotn = 0.0

      do nb = 1, ntot
         if (ikk(2,nb) .eq. 0) go to 520
         dpltot = 0.0
         dqltot = 0.0
         ntyp = kbsdta(1,nb)
         if (ntyp .eq. 5 .or. ntyp .eq. 12) go to  520
         own = owner(nb)
         do l = 1, kntown
            if (own .eq. ownlis(l) .or. ownlis(l) .eq. '###') go to 430
         enddo
         go to  450

  430    continue
         plo = busdta(3,nb)
         qlo = busdta(4,nb)
         pctchg(1) = pctchg(1) + plo
         pctchg(2) = pctchg(2) + qlo
         dp = plo * pctppu
         dq = qlo * pctqpu
         if (dp .ne. 0.0 .or. dq .ne. 0.0) then
            dpltot = dpltot + dp / bmva
            dqltot = dqltot + dq / bmva

            pltoto = pltoto + plo
            qltoto = qltoto + qlo
            pltotn = pltotn + plo + dp
            qltotn = qltotn + qlo + dq

            busdta(3,nb) = plo + dp
            busdta(4,nb) = qlo + dq
            ipctch(1) = ipctch(1) + 1
            if (idswb .ne. 0) then
               write (outbuf, 440) '  ', '  ', bus(nb), base(nb),
     1            zone(nb), owner(nb), plo, qlo, busdta(3,nb),
     2            busdta(4,nb)
  440          format (3x, 2x, a2, 2x, a3, 1x, a8, f6.1, 2x,
     1            'B ', 3x, a2, 3x, a3, 2f9.2, 17x, 2f9.2)
               call prtout(1)
            endif
         endif
C                                 Perform %changes on customer buses.
  450    icb = kbsdta(15,nb)
         do while (icb .gt. 0)
C
C           Exclude type "A" code year "  " and "00" continuation bus
C           data (equivalent or cut system data)
C
            call getchr(1,cbtyp,kbctbl(8,icb))
            call getchr(2,cbkyr,kbctbl(9,icb))
            call getchr(3,cbown,kbctbl(10,icb))
            if ((cbtyp .eq. 'A') .and.
     1          (cbkyr .ne. '01' .and. cbkyr .ne. '02')) go to 500
            do l = 1, kntown
               if (cbown .eq. ownlis(l) .or. ownlis(l) .eq. '###')
     1            go to  480
            enddo
            go to  500
C
C           Found customer bus match
C
  480       plo = bctbl(2,icb)
            qlo = bctbl(3,icb)

            pctchg(1) = pctchg(1) + plo
            pctchg(2) = pctchg(2) + qlo
            dp = plo * pctppu
            dq = qlo * pctqpu

            if (dp .ne. 0.0 .or. dq .ne. 0.0) then
               dpltot = dpltot + dp / bmva
               dqltot = dqltot + dq / bmva
               pltoto = pltoto + plo
               qltoto = qltoto + qlo
               pltotn = pltotn + plo + dp
               qltotn = qltotn + qlo + dq

               bctbl(2,icb) = plo + dp
               bctbl(3,icb) = qlo + dq
               ipctch(2) = ipctch(2) + 1
               k = kbctbl(1,icb)
               call getchr(3,cbown,kbctbl(10,icb))
               if (idswb .ne. 0) then
                  write (outbuf, 490) bus(k), base(k), cbtyp, cbown,
     1               plo, qlo, bctbl(2,icb), bctbl(3,icb)
  490             format (13x, a8, f6.1, 2x, '+', a1, 8x, a3,
     1               2f9.2, 17x, 2f9.2)
                  call prtout (1)
               endif
            endif
  500       icb = bctbl_nxt(icb)
         enddo
  510    continue
         if (dpltot .ne. 0.0 .or. dqltot .ne. 0.0) then
            kt = inp2opt(nb)
            pnetu(kt) = pnetu(kt) - dpltot
            qnetu(kt) = qnetu(kt) - dqltot
            ploadu(kt) = ploadu(kt) + dpltot
            qloadu(kt) = qloadu(kt) + dqltot
            ntbx = ikk(3,nb)
            if (ntbx .gt. 0) then
               ltyp = tbx(1,ntbx)
               if (ltyp .eq. 2 .or. ltyp .eq. 3 .or. ltyp .eq. 5) then
                  tbx(3,ntbx) = tbx(3,ntbx) - dqltot
                  tbx(4,ntbx) = tbx(4,ntbx) - dqltot
               endif
            endif
         endif

  520    continue
      enddo

      write (outbuf, 530) pltoto, qltoto, pltotn, qltotn
  530 format ('0', 2x, 'TOTAL', t28, 'ORIGINAL LOAD', 2f9.2,
     1   7x, 'FINAL LOAD', 2f9.2)
      call prtout (1)
      write (outbuf, 540) ipctch(1), ipctch(2)
  540 format ('0', 2x, i5, 'BUSES', i5, '+ BUSES AFFECTED')
      call prtout (1)
C
C     Squeeze out blanked out words from WORD().
C
  900 jwrd = nwrd
      nwrd = 0
      do i = 1, jwrd
         if (word(i) .ne. ' ') then
            nwrd = nwrd + 1
            word(nwrd) = word(i)
         endif
      enddo
C
C     Check proper control card upon program exit
C
      if (index ('[(/HSC', card) .eq. 0) then
C
C        Meaningless control record
C
         write (errbuf(1), 920) buf(1:80)
  920    format ('Unrecognized command: (', a80, ')')
         if (is_batch .eq. 0) then
            call prterx ('E',1)
         else
            call prterx ('F',1)
         endif
         error = 1
         return
      endif

      outbuf = '0End of %LOAD_CHANGE '
      call prtout(1)
      call forbtm

      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
!***********************************************************************
!                                                                      *
      SUBROUTINE QQSORT(NFILE, NUMBER, KSTART, NAME, KAMAX)
!                                                                      *
!     The list of unique integrals (j,i) is formed in the order of     *
!     increasing symmetry, i.e. with j .le. i.                         *
!     With each integral (INTGRL) there is a pointer (INTPTR)          *
!     that points to the last element in the array of coefficients     *
!     for that integral.                                               *
!                                                                      *
!     On Exit                                                          *
!                                                                      *
!     INTGRL -- array of integrals (packed form of j,i)                *
!     NINTG  -- number of integrals                                    *
!     INTPTR -- array of pointers to last element in list of coeff.    *
!     CNN    -- array of coefficients                                  *
!     NCOEFF -- number of coefficients                                 *
!     JANN   -- row of matrix                                          *
!     JBNN   -- column of matrix                                       *
!                                                                      *
!     Written by Per Jonsson                                           *
!                                                                      *
!***********************************************************************
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:38:04   1/ 6/07
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY: DOUBLE
      USE parameter_def,   ONLY: KEYORB
      USE memory_man
      USE orb_C
      USE default_C
      USE mcpdata_C
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER, INTENT(IN) :: NFILE
      INTEGER, INTENT(IN) :: NUMBER
      INTEGER  :: KSTART
      INTEGER, INTENT(IN) :: KAMAX
      CHARACTER  :: NAME*24
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
      INTEGER, PARAMETER :: KEY = KEYORB
      INTEGER, PARAMETER :: NF = 200
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: I, L, IR, JA, JB, INT, J, INTT
      REAL(DOUBLE) :: CN
!-----------------------------------------------
!
      REWIND (NFILE + 80)

      NCOEFF = NUMBER
!
!   Sort the list
!
!   Allocate storage for all required arrays these arrays are then deallocated
!   in mcp
!
      CALL ALLOC (JANN,   NCOEFF, 'JANN',   'QQSORT')
      CALL ALLOC (JBNN,   NCOEFF, 'JBNN',   'QQSORT')
      CALL ALLOC (INTGRL, NCOEFF, 'INTGRL', 'QQSORT')
      CALL ALLOC (CNN,    NCOEFF, 'CNN',    'QQSORT')
      CALL ALLOC (INTPTR, NCOEFF, 'INTPTR', 'QQSORT')
!
!   Read arrays into memory from NFILE
!
      DO I = 1, NCOEFF
         READ (NFILE + 80) JANN(I), JBNN(I), INTGRL(I), CNN(I)
      END DO
!
!   Sort INTGRL into ascending order using the heapsort algorithm;
!   (Numerical recepies page 231.) move the associated members of
!   of JANN and JBNN in the same
!   manner;
!
      IF (NCOEFF > 1) THEN
!
         L = NCOEFF/2 + 1
         IR = NCOEFF
    2    CONTINUE
         IF (L > 1) THEN
            L = L - 1
            JA = JANN(L)
            JB = JBNN(L)
            INT = INTGRL(L)
            CN = CNN(L)
         ELSE
            JA = JANN(IR)
            JB = JBNN(IR)
            INT = INTGRL(IR)
            CN = CNN(IR)
            JANN(IR) = JANN(1)
            JBNN(IR) = JBNN(1)
            INTGRL(IR) = INTGRL(1)
            CNN(IR) = CNN(1)
            IR = IR - 1
            IF (IR == 1) THEN
               JANN(1) = JA
               JBNN(1) = JB
               INTGRL(1) = INT
               CNN(1) = CN
               GO TO 4
            ENDIF
         ENDIF
         I = L
         J = L + L
    3    CONTINUE
         IF (J <= IR) THEN
            IF (J < IR) THEN
               IF (INTGRL(J) < INTGRL(J+1)) J = J + 1
            ENDIF
            IF (INT < INTGRL(J)) THEN
               JANN(I) = JANN(J)
               JBNN(I) = JBNN(J)
               INTGRL(I) = INTGRL(J)
               CNN(I) = CNN(J)
               I = J
               J = J + J
            ELSE
               J = IR + 1
            ENDIF
            GO TO 3
         ENDIF
         JANN(I) = JA
         JBNN(I) = JB
         INTGRL(I) = INT
         CNN(I) = CN
         GO TO 2
      ENDIF
!

!   Sorting complete; close the file
!
!ww    4 CLOSE (80+NFILE)
    4 CONTINUE
      NINTG = 1
      INTT = INTGRL(1)
!
      DO I = 1, NCOEFF
         IF (INTGRL(I) == INTT) CYCLE
         INTPTR(NINTG) = I - 1
         NINTG = NINTG + 1
         INTT = INTGRL(I)
      END DO

      INTPTR(NINTG) = NCOEFF
!
      DO I = 1, NINTG
         INTGRL(I) = INTGRL(INTPTR(I))
      END DO
!
!  If output option is set dump the data on file
!
      IF (NDUMP == 1) THEN
!
!  If first set of data open the file and print
!  some data to later be able to identify the file
!
!        IF (NFILE.EQ.1.AND.IBLK.EQ.1) THEN
!         print *, ' open sorted ang. file .TB(NF)', NF
!          J = INDEX(NAME,' ')
!          OPEN (UNIT = NF,FILE=NAME(1:J-1)//'.TB',
!     :          STATUS='UNKNOWN',FORM='UNFORMATTED')
!          REWIND (NF)
!        ENDIF
         IF (NFILE == 1) WRITE (NF) NCF, NW, KAMAX
!
!  Print out angular data for this kappa
!
         WRITE (NF) NINTG, NCOEFF
         DO I = 1, NINTG
            WRITE (NF) INTGRL(I), INTPTR(I)
         END DO
         DO I = 1, NCOEFF
            WRITE (NF) CNN(I), JANN(I), JBNN(I)
         END DO
      ENDIF
      CALL DALLOC (JANN,   'JANN',   'QQSORT')
      CALL DALLOC (JBNN,   'JBNN',   'QQSORT')
      CALL DALLOC (INTGRL, 'INTGRL', 'QQSORT')
      CALL DALLOC (CNN,    'CNN',    'QQSORT')
      CALL DALLOC (INTPTR, 'INTPTR', 'QQSORT')
!
!  Has all data been processed? If so close
!  the file
!
!     IF (NFILE.EQ.KAMAX) CLOSE (NF)

!  Debug output
!
!      IF (KSTART.EQ.1) THEN
!        NNNN = 120
!      ELSE
!        NNNN = 140
!      ENDIF
!      WRITE(NNNN+NFILE,*) ' Data for '
!      WRITE(NNNN+NFILE,*) NINTG, ' Integrals'
!      DO I = 1,NINTG
!        IA = INTGRL(I)/KEY
!        IB = MOD(INTGRL(I),KEY)
!        WRITE(NNNN+NFILE,302) NP(IA),NH(IA),NP(IB),NH(IB),INTPTR(I)
!      ENDDO
!      WRITE(NNNN+NFILE,*) NCOEFF,' Coefficients'
!      WRITE(NNNN+NFILE,'(F12.8,2I6)')
!     :     (CNN(I),JANN(I),JBNN(I),I=1,NCOEFF)
!
      RETURN
!
  301 FORMAT(' T_[',1I2,',',1I2,']',' (',1I2,1A2,',',1I2,1A2,') = ',1P,D19.12)
  302 FORMAT(' (',1I2,1A2,',',1I2,1A2,')',I6)
      RETURN
!
      END SUBROUTINE QQSORT

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
      program main
        integer*2 i2_1
        parameter (i2_1 = '1'x)
        ! pads to '0001'x and resolves to Int 1

        integer*2 i2_2
        parameter (i2_2 = '10'x)
        ! pads to '0010'x and resolves to Int 16

        integer*2 i2_3
        parameter (i2_3 = '100'x)
        ! pads to '0100'x and resolves to Int 256

        integer*2 i2_4
        parameter (i2_4 = '1000'x)
        ! doesn't pad and resolves to Int 4096

        integer*2 i2_5
        parameter (i2_5 = '7fff'x)
        ! doesn't pad and resolves to Int 32767
        ! the largest number with a positive sign
        !   (most significant bit = 0) for kind 2
        ! in other words, the highest number INTEGER*2 can take

        integer*2 i2_6
        parameter (i2_6 = '8000'x)
        ! doesn't pad and resolves to Int -32768
        ! the largest number with a negative sign
        !   (most significant bit = 1) for kind 2
        ! in other words, the lowest number INTEGER*2 can take

        integer*2 i2_7
        parameter (i2_7 = '8001'x)
        ! doesn't pad and resolves to Int -32767

        integer*2 i2_8
        parameter (i2_8 = '10000'x)
        ! truncates to '0000'x and resolves to Int 0

        integer*2 i2_9
        parameter (i2_9 = '10010'x)
        ! trunactes to '0010'x and resolves to Int 16

        integer*2 i2_10
        parameter (i2_10 = '00010'x)
        ! truncates to '0010'x and resolves to Int 16
      end

max_stars_repo_name: camfort/fortran-vars

max_issues_repo_name: camfort/fortran-vars

max_forks_repo_name: camfort/fortran-vars