text:
c ==============================================================================
c =                                                                            =
c = Author: Sylwester Wysocki <sw143@wp.pl>                                    =
c = Created on: 2008                                                           =
c =                                                                            =
c = This is free and unencumbered software released into the public domain.    =
c =                                                                            =
c = Anyone is free to copy, modify, publish, use, compile, sell, or            =
c = distribute this software, either in source code form or as a compiled      =
c = binary, for any purpose, commercial or non-commercial, and by any          =
c = means.                                                                     =
c =                                                                            =
c = In jurisdictions that recognize copyright laws, the author or authors      =
c = of this software dedicate any and all copyright interest in the            =
c = software to the public domain. We make this dedication for the benefit     =
c = of the public at large and to the detriment of our heirs and               =
c = successors. We intend this dedication to be an overt act of                =
c = relinquishment in perpetuity of all present and future rights to this      =
c = software under copyright law.                                              =
c =                                                                            =
c = THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,            =
c = EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF         =
c = MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.     =
c = IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR          =
c = OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,      =
c = ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR      =
c = OTHER DEALINGS IN THE SOFTWARE.                                            =
c =                                                                            =
c = For more information, please refer to <https://unlicense.org>              =
c =                                                                            =
c ==============================================================================

c================================================
c==   Wczytuje geometrie molekuly z pliku      ==
c==                                            ==
c==     u_in = uchwyt do pliku (WE)            ==
c==      Nat = ile atomow (WY)                 ==
c==   GEOM = wspolrzedne atomow i l.at. (WY)   ==
c================================================

      subroutine ReadGeometry (u_in)
        IMPLICIT NONE

        integer u_in                          ! uchwyt do pliku
        integer Nat
        real*8 GEOM(1000,4)                   ! 3 wsp. i liczba at.

        character*10 tag
        real*8 x,y,z
        character*2  sym                      ! symbol pierwiastka
        integer i
        integer SymToZ
        integer occupied                      ! ile zajetych orbitali

        COMMON /GEOM/ GEOM,Nat
        COMMON /OCCUPIED/ occupied

c       ======================
c       = sprawdz tag 'GEOM' =
c       ======================

        read (u_in,'(A4)') tag

        if (tag.ne.'GEOM') then
          write (*,*) 'BLAD! Niepoprawna geometria!'
          stop
        endif

c       ====================
c       = czytaj geometrie =
c       ====================

        Nat = 0                                    ! zliczaj atomy
        occupied = 0                               ! licz orbitale
  100   read (u_in,'(A8)') tag
        if (tag.eq.'END_GEOM') goto 200

        backspace (u_in)
        read (u_in,*) sym, x,y,z

        Nat = Nat + 1

        GEOM(Nat,1) = x/0.5292
        GEOM(Nat,2) = y/0.5292
        GEOM(Nat,3) = z/0.5292
        GEOM(Nat,4) = SymToZ (sym)
        occupied = occupied + GEOM(Nat,4)          ! dodaj liczbe atom.
        goto 100                                   ! nastepny atom
  200   continue

c       ==========================
c       = sprwadz liczbe elektr. =
c       ==========================

        if (MOD(occupied,2).eq.1) then
          write (*,*) 'Blad! Nieparzysta liczba elektronow'
          stop
        endif
        occupied = occupied / 2
      end subroutine

c================================================
c==   Wczytuje baze funkcyjna z pliku          ==
c==                                            ==
c==     u_in = uchwyt do pliku z baza  (WE)    ==
c==     Base = baza funkcyjna w postaci: (WY)  ==
c==            (wyk,C)                         ==
c==              C = wsp. w rozwinieciu        ==
c==            wyk = wykladnik                 ==
c== BaseIndex = zakres funkcji w bazie dla (WY)==
c==             zadanej liczby atomowej        ==
c================================================

      subroutine ReadBase (BaseSymbol)
        IMPLICIT NONE

        real*8 Base(1000,2)
        integer BaseIndex(100)
        character*10 BaseSymbol

        integer u_in
        parameter (u_in=35)

        character bufor                            ! bufor do czytania
        integer La                                 ! liczba atomowa
        character*2 Sym                            ! symbol
        integer charge

        character*2 typ                            ! typ orbitalu
        integer l                                  ! poboczna licz. kw.
        integer IleFun                             ! ile funkcji
        real*8 cos

        integer IleOrb                             ! ile orbitali

        real*8 wyk(20), c(20),c2(20)               ! prym. GTO
        integer index                              ! index w bazie
        integer SymToZ                             ! symbol to l.at.
        integer i,k

        character*10 symbol                        ! symbol bazy
        character*10 plik                          ! nazwa pliku

        COMMON /BASE/ Base, BaseIndex

c       ============================
c       = szukaj bazy w bibliotece =
c       ============================

        open (u_in, file='basis.lib')
  500   read (u_in,*,end=510) symbol, plik
        if (symbol.eq.BaseSymbol) goto 520
        goto 500

  510   close (u_in)
        write (*,*) 'Blad! Brak bazy w biliotece!'
        stop

  520   close (u_in)
        open (u_in, file=plik)                      ! otworz plik z baza

        index = 1                                   ! index w bazie

c       ========================================
c       ====== kolejny pierwiastek  ============
c       ========================================

  100   read (u_in, '(A)',end=200) bufor            ! czytaj 1 znak
        if (bufor.ne.'!') then                      ! czy to komentarz?

          backspace (u_in)                          ! wroc na pocz. linijki
          read (u_in,*) sym, charge

          La = SymToZ (sym)                         ! pobierz liczbe at.
          BaseIndex (La) = index                    ! zachowaj polozenie
                                                    ! pierw. w bazie
          IleOrb = 0                                ! zliczaj orbitale
          index = index + 1                         ! zostaw miejsce
                                                    ! na liczbe orbitali
c         ====================
c         = kolejny orbital  =
c         ====================

  300     read (u_in,*) bufor
          if (bufor.eq.'*') then                    ! czy nastepny pierw?
            k = BaseIndex (La)                      ! poczatek pierwiastka
            Base(k,1) = La                          ! zapisz liczb. atm.
            Base(k,2) = IleOrb                      ! oraz ilosc orbiali
            goto 100                                ! czytaj nas. pierw.
          endif

          backspace (u_in)
          read (u_in,*) typ, IleFun, cos            ! kolejny orbital

c         ========================
c         = orbitale S i P razem =
c         ========================

          if (typ.eq.'SP') then
            IleOrb = IleOrb + 2
            do i=1,IleFun
              read (u_in,*) wyk(i), c(i),c2(i)
            enddo

            Base(index,1) = 0
            Base(index,2) = IleFun
            index = index + 1

            do i=1,IleFun
              Base(index,1) = wyk(i)
              Base(index,2) = c(i)
              index = index + 1
            enddo

            Base(index,1) = 1
            Base(index,2) = IleFun
            index = index + 1

            do i=1,IleFun
              Base(index,1) = wyk(i)
              Base(index,2) = c2(i)
              index = index + 1
            enddo

          else

c         =======================
c         = pozostale typy orb. =
c         =======================

            IleOrb = IleOrb + 1
            do i=1,IleFun
              read (u_in,*) wyk(i), c(i)
            enddo

            if (typ.eq.'S') l=0
            if (typ.eq.'P') l=1
            if (typ.eq.'D') l=2
            if (typ.eq.'F') l=3

            Base(index,1) = l
            Base(index,2) = IleFun
            index = index + 1

            do i=1,IleFun
              Base(index,1) = wyk(i)
              Base(index,2) = c(i)
              index = index + 1
            enddo
          endif
          goto 300                                 ! nastepny orbital
        endif

        goto 100
  200   continue

        close (u_in)
      end subroutine


text:
program main
	USE Plotter
!	Функции для графиков
	USE plotter_functions
	implicit none
!	Массив, в котором хранятся указатели.
!	Тип нужен потому что Фортран не умеет делать массивы указателей,
!	А вот указатели на массивы - умеет
	type(functions),dimension(:),allocatable::multiple_graphic
!	Мы собираемся отправить 2 функции поверхности
	allocate(multiple_graphic(2))
	multiple_graphic(1)%function_surface=>first_surface
	multiple_graphic(2)%function_surface=>second_surface
!       Шаг по осям
        call Set_Step_X(1.)
        call Set_Step_Y(1.)
!       Передаём в функцию Размер окна, ограничения по осям Х и У и указатели с их количеством
	call DrawGraph(960,720,-5.,5.,-5.,5.,multiple_graphic,2)
end

max_stars_repo_name: LibPlotter/Plotter

max_issues_repo_name: LibPlotter/Plotter

max_forks_repo_name: LibPlotter/Plotter


text:
SUBROUTINE bigvoi
USE PARACLUSTER
USE CLUSTER
USE POTENTIAL
USE ENFORCE

IMPLICIT NONE
INTEGER :: i,j
REAL*8 :: xij,yij,zij,dij2,dc55
!
      dc55 = cutoff_end * cutoff_end
!
nv4(1:natom)=0.d0
!
      DO i=1,natom-1
        DO j=i+1,natom
           xij=x(j)+u(j)-x(i)-u(i)
           yij=y(j)+v(j)-y(i)-v(i)
           zij=z(j)+w(j)-z(i)-w(i)

!------------------------------------------------------------------------------
! PBC for no_cluster system
!------------------------------------------------------------------------------
IF ( wires ) THEN
!
          IF (ABS(zij+pbcz)<ABS(zij)) zij=zij+pbcz
          IF (ABS(zij-pbcz)<ABS(zij)) zij=zij-pbcz
ENDIF
!
IF( surface ) THEN
           if (abs(xij+pbcx)<abs(xij)) xij=xij+pbcx
           if (abs(xij-pbcx)<abs(xij)) xij=xij-pbcx
           if (abs(yij+pbcy)<abs(yij)) yij=yij+pbcy
           if (abs(yij-pbcy)<abs(yij)) yij=yij-pbcy
ENDIF
!
IF ( bulk ) THEN
           if (abs(xij+pbcx)<abs(xij)) xij=xij+pbcx
           if (abs(xij-pbcx)<abs(xij)) xij=xij-pbcx
           if (abs(yij+pbcy)<abs(yij)) yij=yij+pbcy
           if (abs(yij-pbcy)<abs(yij)) yij=yij-pbcy
           IF (ABS(zij+pbcz)<ABS(zij)) zij=zij+pbcz
           IF (ABS(zij-pbcz)<ABS(zij)) zij=zij-pbcz
ENDIF
!------------------------------------------------------------------------------
!
! Verlet cage system
! dc55 > cutz(itype(i), itype(j))
!
!nv4 == number of neighbours for atom i in a radius dc55
!iv4 == label to identify the number of the neighbour
!
!
           dij2=xij*xij+yij*yij+zij*zij
!
           IF (dij2.LT.dc55) THEN
              nv4(i)=nv4(i)+1
              iv4(nv4(i),i)=j
           ENDIF
!
        ENDDO
!
        IF (nv4(i).GE.nvsiz) THEN
           WRITE (*,'(i4,1x,f6.3,a18)')i,nv4(i),'Too many NN in bigvoi'
           STOP
        ENDIF

      ENDDO
!
End SUBROUTINE bigvoi

max_stars_repo_name: kcl-tscm/LoDiS

max_issues_repo_name: kcl-tscm/LoDiS

max_forks_repo_name: kcl-tscm/LoDiS


text:
      SUBROUTINE PBDMATADD( ICONTXT, MODE, M, N, ALPHA, A, LDA, BETA, B,
     $                      LDB )
*
*  -- PB-BLAS routine (version 2.1) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory.
*     April 28, 1996
*
*     .. Scalar Arguments ..
      CHARACTER*1        MODE
      INTEGER            ICONTXT, LDA, LDB, M, N
      DOUBLE PRECISION   ALPHA, BETA
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
*     ..
*
*  Purpose
*  =======
*
*  PBDMATADD performs the matrix add operation B := alpha*A + beta*B,
*  where alpha and beta are scalars, and A and B are m-by-n
*  upper/lower trapezoidal matrices, or rectangular matrices.
*
*  Arguments
*  =========
*
*  ICONTXT (input) INTEGER
*          ICONTXT is the BLACS mechanism for partitioning communication
*          space.  A defining property of a context is that a message in
*          a context cannot be sent or received in another context.  The
*          BLACS context includes the definition of a grid, and each
*          process' coordinates in it.
*
*  MODE    (input) CHARACTER*1
*          Specifies the part of the matrix A, or (conjugate) transposed
*          matrix A to be added to the matrix B,
*          = 'U':  Upper triangular part
*                  up(B) = alpha*up(A) + beta*up(B)
*          = 'L':  Lower triangular part
*                  lo(B) = alpha*lo(A) + beta*lo(B)
*          = 'T':  Transposed matrix A
*                  B = alpha*A**T + beta*B
*          = 'C':  Conjugate transposed matrix A
*                  B = alpha*A**H + beta*B
*          Otherwise: B = alpha*A + beta*B
*              if M = LDA = LDB: use one BLAS loop
*              if MODE = 'V' :   columnwise copy using BLAS if possible
*              else :            use double loops
*
*  M       (input) INTEGER
*          M specifies the number of columns of the matrix A if
*          MODE != 'T'/'C', and it specifies the number of rows of the
*          matrix A otherwise.  It also specifies the number of rows of
*          the matrix B. M >= 0.
*
*  N       (input) INTEGER
*          N specifies the number of rows of the matrix A if
*          MODE != 'T'/'C', and it specifies the number of columns of
*          the matrix A otherwise. It also specifies the number of
*          columns of the matrix B. N >= 0.
*
*  ALPHA   (input) DOUBLE PRECISION
*          ALPHA specifies the scalar alpha.
*
*  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
*          The m by n matrix A if MODE != 'T'/'C'.
*          If MODE = 'U', only the upper triangle or trapezoid is
*          accessed; if MODE = 'L', only the lower triangle or
*          trapezoid is accessed.  Otherwise all m-by-n data matrix
*          is accessed.
*          And the n by m matrix A if MODE = 'T'/'C'.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,M) if
*          MODE != 'T'/'C'.  And LDA >= max(1,N) if MODE = 'T'/'C'.
*
*  BETA    (input) DOUBLE PRECISION
*          BETA specifies the scalar beta.
*
*  B       (input) DOUBLE PRECISION array, dimension (LDB,N)
*          On exit, B = alpha*A + beta*B
*
*  LDB     (input) INTEGER
*          The leading dimension of the array B.  LDB >= max(1,M).
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO,          ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0)
*     ..
*     .. Local Scalars ..
      INTEGER            I, J
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. External Subroutines ..
      EXTERNAL           DSCAL, DCOPY, DAXPY
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MIN
*     ..
*     .. Executable Statements ..
*
      IF( M.LE.0 .OR. N.LE.0 .OR. ( ALPHA.EQ.ZERO.AND.BETA.EQ.ONE ) )
     $   RETURN
*
*     A is upper triangular or upper trapezoidal,
*
      IF( LSAME( MODE, 'U' ) ) THEN
         IF( ALPHA.EQ.ZERO ) THEN
            IF( BETA.EQ.ZERO ) THEN
               DO 20 J = 1, N
                  DO 10 I = 1, MIN( J, M )
                     B( I, J ) = ZERO
   10             CONTINUE
   20          CONTINUE
            ELSE
               DO 40 J = 1, N
                  DO 30 I = 1, MIN( J, M )
                     B( I, J ) = BETA * B( I, J )
   30             CONTINUE
   40          CONTINUE
            END IF
*
         ELSE IF( ALPHA.EQ.ONE ) THEN
            IF( BETA.EQ.ZERO ) THEN
               DO 60 J = 1, N
                  DO 50 I = 1, MIN( J, M )
                     B( I, J ) = A( I, J )
   50             CONTINUE
   60          CONTINUE
            ELSE IF( BETA.EQ.ONE ) THEN
               DO 80 J = 1, N
                  DO 70 I = 1, MIN( J, M )
                     B( I, J ) = A( I, J ) + B( I, J )
   70             CONTINUE
   80          CONTINUE
            ELSE
               DO 100 J = 1, N
                  DO 90 I = 1, MIN( J, M )
                     B( I, J ) = A( I, J ) + BETA * B( I, J )
   90             CONTINUE
  100          CONTINUE
            END IF
*
         ELSE
            IF( BETA.EQ.ZERO ) THEN
               DO 120 J = 1, N
                  DO 110 I = 1, MIN( J, M )
                     B( I, J ) = ALPHA * A( I, J )
  110             CONTINUE
  120          CONTINUE
            ELSE IF( BETA.EQ.ONE ) THEN
               DO 140 J = 1, N
                  DO 130 I = 1, MIN( J, M )
                     B( I, J ) = ALPHA * A( I, J ) + B( I, J )
  130             CONTINUE
  140          CONTINUE
            ELSE
               DO 160 J = 1, N
                  DO 150 I = 1, MIN( J, M )
                     B( I, J ) = ALPHA * A( I, J ) + BETA * B( I, J )
  150             CONTINUE
  160          CONTINUE
            END IF
         END IF
*
*     A is lower triangular or upper trapezoidal,
*
      ELSE IF( LSAME( MODE, 'L' ) ) THEN
         IF( ALPHA.EQ.ZERO ) THEN
            IF( BETA.EQ.ZERO ) THEN
               DO 180 J = 1, N
                  DO 170 I = J, M
                     B( I, J ) = ZERO
  170             CONTINUE
  180          CONTINUE
            ELSE
               DO 200 J = 1, N
                  DO 190 I = J, M
                     B( I, J ) = BETA * B( I, J )
  190             CONTINUE
  200          CONTINUE
            END IF
*
         ELSE IF( ALPHA.EQ.ONE ) THEN
            IF( BETA.EQ.ZERO ) THEN
               DO 220 J = 1, N
                  DO 210 I = J, M
                     B( I, J ) = A( I, J )
  210             CONTINUE
  220          CONTINUE
            ELSE IF( BETA.EQ.ONE ) THEN
               DO 240 J = 1, N
                  DO 230 I = J, M
                     B( I, J ) = A( I, J ) + B( I, J )
  230             CONTINUE
  240          CONTINUE
            ELSE
               DO 260 J = 1, N
                  DO 250 I = J, M
                     B( I, J ) = A( I, J ) + BETA * B( I, J )
  250             CONTINUE
  260          CONTINUE
            END IF
*
         ELSE
            IF( BETA.EQ.ZERO ) THEN
               DO 280 J = 1, N
                  DO 270 I = J, M
                     B( I, J ) = ALPHA * A( I, J )
  270             CONTINUE
  280          CONTINUE
            ELSE IF( BETA.EQ.ONE ) THEN
               DO 300 J = 1, N
                  DO 290 I = J, M
                     B( I, J ) = ALPHA * A( I, J ) + B( I, J )
  290             CONTINUE
  300          CONTINUE
            ELSE
               DO 320 J = 1, N
                  DO 310 I = J, M
                     B( I, J ) = ALPHA * A( I, J ) + BETA * B( I, J )
  310             CONTINUE
  320          CONTINUE
            END IF
         END IF
*
*     If MODE = 'Transpose'/'Conjugate'
*
      ELSE IF( LSAME( MODE, 'T' ) .OR. LSAME( MODE, 'C' ) ) THEN
         IF( ALPHA.EQ.ZERO ) THEN
            IF( BETA.EQ.ZERO ) THEN
               DO 340 J = 1, N
                  DO 330 I = 1, M
                     B( I, J ) = ZERO
  330             CONTINUE
  340          CONTINUE
            ELSE
               DO 360 J = 1, N
                  DO 350 I = 1, M
                     B( I, J ) = BETA * B( I, J )
  350             CONTINUE
  360          CONTINUE
            END IF
*
         ELSE IF( ALPHA.EQ.ONE ) THEN
            IF( BETA.EQ.ZERO ) THEN
               DO 380 J = 1, N
                  DO 370 I = 1, M
                     B( I, J ) = A( J, I )
  370             CONTINUE
  380          CONTINUE
            ELSE IF( BETA.EQ.ONE ) THEN
               DO 400 J = 1, N
                  DO 390 I = 1, M
                     B( I, J ) = A( J, I ) + B( I, J )
  390             CONTINUE
  400          CONTINUE
            ELSE
               DO 420 J = 1, N
                  DO 410 I = 1, M
                     B( I, J ) = A( J, I ) + BETA * B( I, J )
  410             CONTINUE
  420          CONTINUE
            END IF
*
         ELSE
            IF( BETA.EQ.ZERO ) THEN
               DO 440 J = 1, N
                  DO 430 I = 1, M
                     B( I, J ) = ALPHA * A( J, I )
  430             CONTINUE
  440          CONTINUE
            ELSE IF( BETA.EQ.ONE ) THEN
               DO 460 J = 1, N
                  DO 450 I = 1, M
                     B( I, J ) = ALPHA * A( J, I ) + B( I, J )
  450             CONTINUE
  460          CONTINUE
            ELSE
               DO 480 J = 1, N
                  DO 470 I = 1, M
                     B( I, J ) = ALPHA * A( J, I ) + BETA * B( I, J )
  470             CONTINUE
  480          CONTINUE
            END IF
         END IF
*
*     Other cases (for genral matrix additions)
*
      ELSE
         IF( ALPHA.EQ.ZERO ) THEN
            IF( BETA.EQ.ZERO ) THEN
               DO 500 J = 1, N
                  DO 490 I = 1, M
                     B( I, J ) = ZERO
  490             CONTINUE
  500          CONTINUE
*
            ELSE
               IF( M.EQ.LDB ) THEN
                  CALL DSCAL( M*N, BETA, B( 1, 1 ), 1 )
               ELSE IF( LSAME( MODE, 'V' ) ) THEN
                  DO 510 J = 1, N
                     CALL DSCAL( M, BETA, B( 1, J ), 1 )
  510             CONTINUE
               ELSE
                  DO 530 J = 1, N
                     DO 520 I = 1, M
                        B( I, J ) = BETA * B( I, J )
  520                CONTINUE
  530             CONTINUE
               END IF
            END IF
*
         ELSE IF( ALPHA.EQ.ONE ) THEN
            IF( BETA.EQ.ZERO ) THEN
               IF( M.EQ.LDA .AND. M.EQ.LDB ) THEN
                  CALL DCOPY( M*N, A( 1, 1 ), 1, B( 1, 1 ), 1 )
               ELSE IF( LSAME( MODE, 'V' ) ) THEN
                  DO 540 J = 1, N
                     CALL DCOPY( M, A( 1, J ), 1, B( 1, J ), 1 )
  540             CONTINUE
               ELSE
                  DO 560 J = 1, N
                     DO 550 I = 1, M
                        B( I, J ) = A( I, J )
  550                CONTINUE
  560             CONTINUE
               END IF
*
            ELSE IF( BETA.EQ.ONE ) THEN
               DO 580 J = 1, N
                  DO 570 I = 1, M
                     B( I, J ) = A( I, J ) + B( I, J )
  570             CONTINUE
  580          CONTINUE
*
            ELSE
               DO 600 J = 1, N
                  DO 590 I = 1, M
                     B( I, J ) = A( I, J ) + BETA * B( I, J )
  590             CONTINUE
  600          CONTINUE
            END IF
*
         ELSE
            IF( BETA.EQ.ZERO ) THEN
               DO 620 J = 1, N
                  DO 610 I = 1, M
                     B( I, J ) = ALPHA * A( I, J )
  610             CONTINUE
  620          CONTINUE
*
            ELSE IF( BETA.EQ.ONE ) THEN
               IF( M.EQ.LDA .AND. M.EQ.LDB ) THEN
                  CALL DAXPY( M*N, ALPHA, A( 1, 1 ), 1, B( 1, 1 ), 1 )
               ELSE IF( LSAME( MODE, 'V' ) ) THEN
                  DO 630 J = 1, N
                     CALL DAXPY( M, ALPHA, A( 1, J ), 1, B( 1, J ), 1 )
  630             CONTINUE
               ELSE
                  DO 650 J = 1, N
                     DO 640 I = 1, M
                        B( I, J ) = ALPHA * A( I, J ) + B( I, J )
  640                CONTINUE
  650             CONTINUE
               END IF
*
            ELSE
               DO 670 J = 1, N
                  DO 660 I = 1, M
                     B( I, J ) = ALPHA * A( I, J ) + BETA * B( I, J )
  660             CONTINUE
  670          CONTINUE
            END IF
         END IF
      END IF
*
      RETURN
*
*     End of PBDMATADD
*
      END

ext: f

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack


! { dg-do run }
program nint_1
  if (int(anint(8388609.0)) /= 8388609) STOP 1
  if (int(anint(0.49999997)) /= 0) STOP 2
  if (nint(8388609.0) /= 8388609) STOP 3
  if (nint(0.49999997) /= 0) STOP 4
  if (int(dnint(4503599627370497.0d0),8) /= 4503599627370497_8) STOP 5
  if (int(dnint(0.49999999999999994d0)) /= 0) STOP 6
  if (int(anint(-8388609.0)) /= -8388609) STOP 7
  if (int(anint(-0.49999997)) /= 0) STOP 8
  if (nint(-8388609.0) /= -8388609) STOP 9
  if (nint(-0.49999997) /= 0) STOP 10
  if (int(dnint(-4503599627370497.0d0),8) /= -4503599627370497_8) STOP 11
  if (int(dnint(-0.49999999999999994d0)) /= 0) STOP 12
end program nint_1


module GwfInterfaceModelModule
  use KindModule, only: I4B, DP
  use ConstantsModule, only: DZERO
  use MemoryManagerModule, only: mem_allocate
  use MemoryHelperModule, only: create_mem_path
  use NumericalModelModule, only: NumericalModelType, GetNumericalModelFromList
  use GwfModule, only: GwfModelType, CastAsGwfModel
  use Xt3dModule, only: xt3d_cr
  use GwfBuyModule, only: buy_cr
  use GridConnectionModule
  use BaseDisModule
  use GwfDisuModule
  use GwfNpfModule
  use GwfNpfOptionsModule
  use GwfNpfGridDataModule
  use GwfBuyInputDataModule
  use GwfOcModule
  implicit none
  private

  !> The GWF Interface Model is a utility to calculate the solution's
  !! exchange coefficients from the interface between a GWF model and
  !! its GWF neighbors. The interface model itself will not be part
  !! of the solution, it is not being solved.
  !! Patching (a part of the) discretizations of two GWF models in a
  !! general way, e.g. DIS+DIS with refinement, requires the resulting
  !< discretization to be of type DISU.
  type, public, extends(GwfModelType) :: GwfInterfaceModelType
    class(GridConnectionType), pointer    :: gridConnection => null() !< The grid connection class will provide the interface grid
    class(GwfModelType), private, pointer :: owner => null()          !< the real GWF model for which the exchange coefficients
                                                                      !! are calculated with this interface model
  contains
    procedure, pass(this) :: gwfifm_cr
    procedure :: model_df => gwfifm_df
    procedure :: model_ar => gwfifm_ar
    procedure :: model_da => gwfifm_da

    ! private
    procedure, private, pass(this) :: setNpfOptions
    procedure, private, pass(this) :: setNpfGridData
    procedure, private, pass(this) :: setBuyData
  end type

contains

  !> @brief set up the interface model, analogously to what
  !< happens in gwf_cr
  subroutine gwfifm_cr(this, name, iout, gridConn)
    class(GwfInterfaceModelType) :: this        !< the GWF interface model
    character(len=*), intent(in)  :: name                      !< the interface model's name
    integer(I4B), intent(in) :: iout                           !< the output unit
    class(GridConnectionType), pointer, intent(in) :: gridConn !< the grid connection for creating a DISU
    ! local
    class(*), pointer :: modPtr

    this%memoryPath = create_mem_path(name)
    call this%allocate_scalars(name)

    this%iout = iout

    this%gridConnection => gridConn
    modPtr => this%gridConnection%model
    this%owner => CastAsGwfModel(modPtr)

    this%innpf = huge(1_I4B)
    this%inewton = this%owner%inewton
    this%inewtonur = this%owner%inewtonur

    if (this%owner%inbuy > 0) then
      this%inbuy = huge(1_I4B)
    end if

    ! create discretization and packages
    call disu_cr(this%dis, this%name, -1, this%iout)
    call npf_cr(this%npf, this%name, this%innpf, this%iout)
    call xt3d_cr(this%xt3d, this%name, this%innpf, this%iout)
    call buy_cr(this%buy, this%name, this%inbuy, this%iout)

  end subroutine gwfifm_cr

  !> @brief Define, mostly DISU and the NPF package
  !< for this interface model
  subroutine gwfifm_df(this)
    class(GwfInterfaceModelType) :: this !< the GWF interface model
    ! local
    type(GwfNpfOptionsType) :: npfOptions
    type(GwfBuyInputDataType) :: buyData
    class(*), pointer :: disPtr

    this%moffset = 0

    ! define DISU
    disPtr => this%dis
    call this%gridConnection%getDiscretization(CastAsDisuType(disPtr))

    ! define NPF package
    call npfOptions%construct()
    call this%setNpfOptions(npfOptions)
    call this%npf%npf_df(this%dis, this%xt3d, 0, npfOptions)
    call npfOptions%destroy()

    ! define BUY package
    if (this%inbuy > 0) then
      call buyData%construct(this%owner%buy%nrhospecies)
      call this%setBuyData(buyData)
      call this%buy%buy_df(this%dis, buyData)
      call buyData%destruct()
    end if

    this%neq = this%dis%nodes
    this%nja = this%dis%nja
    this%ia  => this%dis%con%ia
    this%ja  => this%dis%con%ja

    call this%allocate_arrays()

  end subroutine gwfifm_df

  !> @brief allocate and read the packages
  !<
  subroutine gwfifm_ar(this)
    class(GwfInterfaceModelType) :: this !< the GWF interface model
    ! local
    type(GwfNpfGridDataType) :: npfGridData

    call npfGridData%construct(this%dis%nodes)
    call this%setNpfGridData(npfGridData)
    call this%npf%npf_ar(this%ic, this%ibound, this%x, npfGridData)
    call npfGridData%destroy()

    if (this%inbuy > 0) call this%buy%buy_ar(this%npf, this%ibound)

  end subroutine gwfifm_ar

  !> @brief Clean up
  !<
  subroutine gwfifm_da(this)
  use MemoryManagerModule, only: mem_deallocate
    class(GwfInterfaceModelType) :: this !< the GWF interface model

    ! -- Internal flow packages deallocate
    call this%dis%dis_da()
    call this%npf%npf_da()
    call this%xt3d%xt3d_da()
    call this%buy%buy_da()
    !
    ! -- Internal package objects
    deallocate(this%dis)
    deallocate(this%npf)
    deallocate(this%xt3d)
    !
    ! -- Scalars
    call mem_deallocate(this%inic)
    call mem_deallocate(this%inoc)
    call mem_deallocate(this%inobs)
    call mem_deallocate(this%innpf)
    call mem_deallocate(this%inbuy)
    call mem_deallocate(this%insto)
    call mem_deallocate(this%incsub)
    call mem_deallocate(this%inmvr)
    call mem_deallocate(this%inhfb)
    call mem_deallocate(this%ingnc)
    call mem_deallocate(this%iss)
    call mem_deallocate(this%inewtonur)
    !
    ! -- NumericalModelType
    call this%NumericalModelType%model_da()

  end subroutine

  !> @brief Copy NPF options from the model owning
  !! the interface to the data structure
  !<
  subroutine setNpfOptions(this, npfOptions)
    class(GwfInterfaceModelType) :: this  !< the GWF interface model
    type(GwfNpfOptionsType) :: npfOptions !< the options data to be filled

    ! for now, assuming full homogeneity, so just take
    ! the options from the owning model's npf package
    npfOptions%icellavg = this%owner%npf%icellavg
    npfOptions%ithickstrt = this%owner%npf%ithickstrt
    npfOptions%iperched = this%owner%npf%iperched
    npfOptions%ivarcv = this%owner%npf%ivarcv
    npfOptions%idewatcv = this%owner%npf%idewatcv
    npfOptions%irewet = this%owner%npf%irewet
    npfOptions%wetfct = this%owner%npf%wetfct
    npfOptions%iwetit = this%owner%npf%iwetit
    npfOptions%ihdwet = this%owner%npf%ihdwet

  end subroutine setNpfOptions

  !> @brief Loop over the interface grid and fill the structure
  !! with NPF grid data, copied from the models that participate
  !! in this interface
  !<
  subroutine setNpfGridData(this, npfGridData)
    class(GwfInterfaceModelType) :: this    !< the interface model
    type(GwfNpfGridDataType) :: npfGridData !< grid data to be set
    ! local
    integer(I4B) :: icell, idx
    class(*), pointer :: modelPtr
    class(GwfModelType), pointer :: gwfModel

    ! TODO_MJR: deal with inhomogeneity, for now, we assume
    ! that we can just take the owning model's settings...
    npfGridData%ik22 = this%owner%npf%ik22
    npfGridData%ik33 =  this%owner%npf%ik33
    npfGridData%iwetdry = this%owner%npf%iwetdry
    npfGridData%iangle1 = this%owner%npf%iangle1
    npfGridData%iangle2 = this%owner%npf%iangle2
    npfGridData%iangle3 = this%owner%npf%iangle3
    if (this%npf%ixt3d > 0) then
      npfGridData%iangle1 = 1
      npfGridData%iangle2 = 1
      npfGridData%iangle3 = 1
    end if

    do icell = 1, this%gridConnection%nrOfCells
      idx = this%gridConnection%idxToGlobal(icell)%index
      modelPtr => this%gridConnection%idxToGlobal(icell)%model
      gwfModel => CastAsGwfModel(modelPtr)

      npfGridData%icelltype(icell) = gwfModel%npf%icelltype(idx)
      npfGridData%k11(icell) = gwfModel%npf%k11(idx)
      npfGridData%k22(icell) = gwfModel%npf%k22(idx)
      npfGridData%k33(icell) = gwfModel%npf%k33(idx)

      ! the K rotation angles, or default (0.0)
      if (npfGridData%iangle1 == 1) then
        if (gwfModel%npf%iangle1 == 1) then
          npfGridData%angle1(icell) = gwfModel%npf%angle1(idx)
        else
          npfGridData%angle1(icell) = DZERO
        end if
      end if
      if (npfGridData%iangle2 == 1) then
        if (gwfModel%npf%iangle2 == 1) then
          npfGridData%angle2(icell) = gwfModel%npf%angle2(idx)
        else
          npfGridData%angle2(icell) = DZERO
        end if
      end if
      if (npfGridData%iangle3 == 1) then
        if (gwfModel%npf%iangle3 == 1) then
          npfGridData%angle3(icell) = gwfModel%npf%angle3(idx)
        else
          npfGridData%angle3(icell) = DZERO
        end if
      end if

      ! wetdry parameter, TODO_MJR: where is it ever set to 1??
      if (npfGridData%iwetdry == 1) then
        if (gwfModel%npf%iwetdry == 1) then
          npfGridData%wetdry(icell) = gwfModel%npf%wetdry(idx)
        else
          npfGridData%wetdry(icell) = DZERO
        end if
      end if

    end do

  end subroutine setNpfGridData

  !> @brief Sets the BUY input data from the models that
  !! make up this interface. We adopt everything from the
  !! owning model, but during validation it should be
  !< checked that the models are compatible.
  subroutine setBuyData(this, buyData)
    class(GwfInterfaceModelType) :: this !< the interface model
    type(GwfBuyInputDataType) :: buyData !< the data for the buoyancy package
    ! local
    integer(I4B) :: i

    buyData%denseref = this%owner%buy%denseref
    buyData%iform = this%owner%buy%iform
    buyData%nrhospecies = this%owner%buy%nrhospecies

    do i = 1, buyData%nrhospecies
      buyData%drhodc(i) = this%owner%buy%drhodc(i)
      buyData%crhoref(i) = this%owner%buy%crhoref(i)
      buyData%cmodelname(i) = this%owner%buy%cmodelname(i)
      buyData%cauxspeciesname(i) = this%owner%buy%cauxspeciesname(i)
    end do

  end subroutine setBuyData

end module GwfInterfaceModelModule


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

20100423 19:51:55 nbsp Hello, Rich, and Welcome to the Wiki!   My names Evan, pleased to meet you! Thanks for creating the Verizon entry... are you a manager or owner there?  Its a great first entry, and I encourage you to poke around and edit the rest of the wiki.  Everything here was created by people like you.  Once again, welcome to the wiki!  Users/JabberWokky Evan JabberWokky Edwards

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine frcfac(nowcrp)
c
c     + + + PURPOSE + + +
c    Computes the soil grain friction factor (frcsol(iplane)), the total
c    rill friction factor (frctrl(iplane)), and the equivalent weighting
c    friction factor for rills (frcteq(iplane)).
c
c     Generates the interrill and rill friction factors for the
c     erosion component; computes the time-variant roughness
c     coefficients for overland flow routing.
c
c     Called from IRS.
c     Author(s):
c     Reference in User Guide:
c
c     Changes:
c           1) Common blocks: CLIYR, CRPVR1, & CRPVR5 were not accessed.
c              Reference to them was deleted.
c           2) Since INRCOV & RILCOV are just read in INPUT, used
c              in this routine, and not used anywhere else,  they
c              are no longer set to zero, if their original values
c              are negative.
c           3) RROUGH is handled in a similar fashion, for the same
c              reason.
c           4) Calculations for RILARE(MXPLAN) simplified, and
c              RILARE(MXPLAN) changed to local variable RILLAR.
c     ******************************************************************
c     *  NOTE:  "RILARE" SHOULD BE REMOVED FROM COMMON BLOCK "FFACT".  *
c     ******************************************************************
c           5) SAVE statement, which saves ALL local variables, removed.
c           6) Removed use of variables FRSPRV and FRTPRV because
c              variables FRCSOL and FRCTRL have now been dimensioned
c              to MXPLAN, and we can use the previous plane's value of
c              these in SR PARAM directly.  Also removed them from
c              common /ffact/.  Computation of FRSPRV and FRTPRV in
c              WEPP version 93.2 was in error.    dcf  10/14/93
c
c
c     Version: This module recoded from WEPP version 91.10.
c     Date recoded: 04/22/91 - 4/23/91.
c     Recoded by: Charles R. Meyer.
c
c     + + + KEYWORDS + + +
c
c     + + + ARGUMENT DECLARATIONS + + +
      integer nowcrp
c
c     + + + ARGUMENT DEFINITIONS + + +
c
c     + + + PARAMETERS + + +
      include 'pntype.inc'
      include 'pmxnsl.inc'
      include 'ptilty.inc'
      include 'pmxpln.inc'
      include 'pmxcrp.inc'
      include 'pmxelm.inc'
      include 'pmxres.inc'
      include 'pmxtls.inc'
      include 'pmxtil.inc'
      include 'pmxhil.inc'
      include 'pmxgrz.inc'
c
c     + + + COMMON BLOCKS + + +
      include 'cffact.inc'
c     modify: frcsol(iplane), frctrl(iplane)
c      write: frcteq(iplane)
c
      include 'ccontcv.inc'
c       read: tilseq(mxcrop,mxplan)
c
      include 'ccrpvr2.inc'
c       read: ar(ntype),flivmx(nowcrp(iplane))
c
      include 'ccrpvr3.inc'
c       read: hmax(ntype)
c
      include 'ccover.inc'
c       read: lanuse(mxplan), canhgt(mxplan), inrcov(mxplan),
c       rilcov(mxplan)
c
      include 'ccrpout.inc'
c       read: rescov(mxplan), rrc(mxplan)
c
      include 'ccrpprm.inc'
c       read: itype, itill
c
      include 'cends.inc'
c       read: rspace(mxplan), width(mxplan)
c
      include 'cfripas.inc'
c
      include 'crinpt1.inc'
c       read: wcf(mxplan)
c
c     include 'crinpt1a.inc'
c
      include 'crinpt6.inc'
c       read: rufcov(mxplan), rrough(mxplan)
c
      include 'cslinit.inc'
c       read: rrinit(mxplan)
c
      include 'csolvar.inc'
c       read: sand, clay
c
      include 'cstruc.inc'
c       read: iplane,imodel
c
      include 'ctillge.inc'
c       read: rro
c
      include 'cupdate.inc'
c       read: indxy(mxplan)
c
c     + + + LOCAL VARIABLES + + +
c     real inrfso, inrfo, inrrou, inrfro, inrfco, inrfto, rroinr,
c    1    rrrinr, frrraf, frrrok, frirok, rillar, frires
      real inrfso, inrfo, inrrou, inrfco, rroinr, rrrinr
      integer i1, i2, idplnt
c
c     + + + LOCAL DEFINITIONS + + +
c     inrfso - friction due to soil texture
c     inrfo  - Darcy friction coefficient for interrill surface
c     inrrou -
c     inrfro - friction due to form roughness
c     inrfco - Darcy friction coefficient for interrill surface
c              area
c     inrfto - total Darcy friction factor for the interrill area
c     rroinr - initial random roughness of a freshly tilled
c              surface
c     rrrinr -
c     frccov(iplane) - Darcy friction coefficient for rill surface cover
c     frican(iplane) - friction due to canopy cover
c     flivmx -
c     frlive(iplane) -
c     frrraf - friction due to form roughness
c     frrrok - friction (rill) due to rock and gravel
c     frirok - friction (interrill) due to rock and gravel
c     frrres(iplane) - friction (rill) due to residue
c     frires(iplane) - friction (interrill) due to residue
c     fribas(iplane) - friction due to basal cover and cryptogram
c     fritot - sum of fribas(iplane) and frican under rangeland
c     rillar - rill area
c     i1     -
c     i2     -
c     idplnt -
c
c     + + + END SPECIFICATIONS + + +
c
c      *** L0 IF ***
      if (lanuse(iplane).eq.1) then
c
c       hydraulic roughness for interrill areas
c
        inrfro = 0.0
c
c       friction coefficient for smooth bare soil on an
c       interrill area
c
        inrfso = 4.07
c
c       For the case of a continuous simulation
c
        if(imodel.eq.1)then
          i1 = indxy(iplane)
          i2 = tilseq(nowcrp,iplane)
          if (i2.gt.0) then
            if (indxy(iplane).eq.0)then
              rroinr = rrinit(iplane)
            elseif(resman(i1,i2).eq.10 .or. resman(i1,i2).eq.11) then
              rroinr = rrinit(iplane)
            else
              rroinr = rro(i1,i2)
            end if
          else
            rroinr = rrinit(iplane)
          end if
c
c       ELSE for the case of a single storm simulation
c
        else
          rroinr=rrinit(iplane)
        end if
        rrrinr = rrc(iplane) / rroinr
        if (rrrinr.gt.1.0) rrrinr = 1.0
        inrfo = exp(3.024-5.042*exp(-161.*rroinr))
c
c         friction coefficient for interrill surface roughness
c
        inrrou = 0.5 * (inrfo**1.128) * exp(-3.088*(1.0-rrrinr))
        if (inrrou.lt.inrfso) inrrou = inrfso
        inrfro = inrrou - inrfso
c
        if (inrcov(iplane).gt.0.0) then
c
c         friction coefficient for interrill surface cover
c
          inrfco = 14.5 * (inrcov(iplane)) ** 1.5544
        else
          inrcov(iplane) = 0.0
          inrfco = 0.0
        end if
c
c       New section added  5/16/90 by dcf to account for roughness
c       introduced by living grasses, small grains, and alfalfa
c
        idplnt = itype(nowcrp,iplane)
c
        if (hmax(idplnt).gt.0.0) then
          frlive(iplane) = (canhgt(iplane)/hmax(idplnt)) *
     1        flivmx(idplnt)
        else
          frlive(iplane) = 0.0
        end if
c
c       total friction factor for interrill area
c
        inrfto = inrfro + inrfco + inrfso + frlive(iplane)
c
c       hydraulic roughness for rills
c
        frcsol(iplane) = 1.11
c
c       Gilley's new equation - 5/16/90
c       Changed by Nearing - 12/94
c
        if (rilcov(iplane).gt.0.0) then
          frccov(iplane) = 4.5 * (rilcov(iplane)) ** 1.5544
        else
          rilcov(iplane) = 0.0
          frccov(iplane) = 0.0
        end if
c
c------- total rill friction
        frctrl(iplane) = frccov(iplane) + frlive(iplane) +
     1      frcsol(iplane)
c
c
c     *** L0 ELSEIF ***
      else if (lanuse(iplane).eq.2) then
c
c-------rill friction calculations
c-------friction due to soil texture
c
c       frcsol(iplane)=(3.42**clay(1,iplane)/12.42**sand(1,iplane))**0.5
c
c       Use Gilley's value for bare soil frcsol(iplane) instead of present
c       range equation.
c
        frcsol(iplane) = 1.11
        if (rrough(iplane).gt.0.006) then
          frrraf = 42.76 * (1.0-exp(-77.3*rrough(iplane)))
        else
          frrraf = 42.76 * (1.0-exp(-77.3*0.006))
        end if
c
c       Changes from Mary Kidwell/Mark Weltz follow:
c       frrrok = 1.8467 * wcf(iplane)
c       frrres(iplane) = 113.73 * rescov(iplane) ** 3.0
c
        frrrok = 1.8467 * frokr(iplane) * rokcov(iplane)
        frrres(iplane) = 113.73 * (fresr(iplane)*rescov(iplane)) ** 3.0
        frrbas = 125.91 * (fbasr(iplane)*bascov(iplane)
     1           + fcryr(iplane)*crycov(iplane)) ** 0.8
        frrbas = frrbas/2.0

c       frccov(iplane) = frrres(iplane) + frrrok + frrraf
        frccov(iplane) = frrrok + frrres(iplane) + frrbas + frrraf
c
c-------total rill friction
c       frctrl(iplane) = (frcsol(iplane)+frrraf+frrrok+frrres(iplane)-
c    1      15.868)
c
        frctrl(iplane) = (frcsol(iplane)+frrraf+frrrok+frrres(iplane)+
     1                    frrbas - 15.868)
        if (frctrl(iplane).lt.frcsol(iplane)) frctrl(iplane) =
     1      frcsol(iplane)
c
c-------interrill friction calculations
c
c       inrfso=(3.42**clay(1,iplane)/12.42**sand(1,iplane))**0.5
c
c       Use Gilley's value for bare soil inrfso instead of present
c       range equation.
c
        inrfso = 4.07
c
        if (rrough(iplane).gt.0.006) then
          inrfro = 42.76 * (1.0-exp(-77.3*rrough(iplane)))
        else
          inrfro = 42.76 * (1.0-exp(-77.3*0.006))
        end if
c
c       XXX  Changes from Mary Kidwell and Mark Weltz  3/9/95
c       frirok = 1.8467 * wcf(iplane)
c       frires = 113.73 * rescov(iplane) ** 3.0
c       fribas(iplane) = 125.91 * rufcov(iplane) ** 0.8
        frirok = 1.8467 * froki(iplane)*rokcov(iplane)
        frires = 113.73 * (fresi(iplane)*rescov(iplane)) ** 3.0
        fribas(iplane) = 125.91 * (fbasi(iplane)*bascov(iplane) +
     1                   fcryi(iplane)*crycov(iplane)) ** 0.8
        frican(iplane) = 38.95 * cancov(iplane) ** 0.8
        inrfco = fribas(iplane) + frires + frirok + frican(iplane)
        inrfto = (inrfso+inrfro+frirok+frires+fribas(iplane)+
     1      frican(iplane)-15.868) / 0.64924
        if (inrfto.lt.inrfso) inrfto = inrfso
c
c     *** L0 ENDIF ***
      end if
c
c------compute rill area (RILLAR):
      rillar = width(iplane) / rspace(iplane)
c
c
c     Add the following IF/ELSE so that for the special case where
c     the user has specified that the rill width is equal to the
c     rill spacing (broad sheet flow assumption), then the model will
c     use the interrill friction factors only.    dcf   8/22/94
c
      if(rillar .lt. 1.0) then
c       equivalent total fric factor for composite area
        frcteq(iplane) = inrfto + rillar * (frctrl(iplane)-inrfto)
      else
        frctrl(iplane) = inrfto
        frcsol(iplane) = inrfso
        frcteq(iplane) = inrfto
      endif
c
      return
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
      subroutine readfrag
      use derivheader

c  read the frags.out file
c which was created by the fragmentation program
      open(unit=3,file='frags.out',status='old')

      read(3,*)
      read(3,*)nfrag
      read(3,*)
      read(3,*)
      read(3,*)
      allocate(nat0(nfrag))
      do m=1,nfrag
c nfrag = the number of fragments (read in readdata)
       read(3,*)nat0(m)
c nat0 = the number of real atoms in the fragment (readdata)
      enddo
      maxatom=0
      do m=1,nfrag
       if(nat0(m).gt.maxatom)maxatom=nat0(m)
      enddo
      maxatom=maxatom+maxcaps

      allocate(nat(nfrag,maxatom))
      allocate(ncap(nfrag))
      allocate(natcap(nfrag,maxcaps,2))
      allocate(fact(nfrag,maxcaps))

      read(3,*)
      do m=1,nfrag
       read(3,*)(nat(m,i),i=1,nat0(m))
c nat = the whole-molecule atom numbers in each fragment (readdata)
      enddo
      read(3,*)
      do m=1,nfrag
       ncap(m)=0
c ncap = the number of caps in each fragment
      enddo
c  read in the identity of the atoms connected to the caps
444   continue
      read(3,*,end=445)k1,m1,n1,fac
c k1 the fragment number
c m1 and n1 are the atoms which define the cap
c fac is the weight for m1 versus n1
c ncap, natcap and fact were allocated in readdata
      ncap(k1)=ncap(k1)+1
      natcap(k1,ncap(k1),1)=m1
      natcap(k1,ncap(k1),2)=n1
      fact(k1,ncap(k1))=fac
      go to 444
445   continue

c total number of atoms in each fragment
      allocate(natomfrag(nfrag))
      do m=1,nfrag
      natomfrag(m)=nat0(m)+ncap(m)
      enddo
      close(unit=3)

      return
      end

max_stars_repo_name: mickcollins/SMFA

max_issues_repo_name: mickcollins/SMFA

max_forks_repo_name: mickcollins/SMFAPAC


text:
locate type T at "T.fpp"
locate constant a at "a.fpp"

struct S { x: T } default a

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
      SUBROUTINE CFFTI1 (N,WA,IFAC)
      IMPLICIT NONE
      INTEGER :: N, IFAC(*), I, I1, IB, IDO, IDOT, II, IP, IPM, J, K1
      INTEGER :: L1, L2, LD, NF, NL, NQ, NR, NTRY, NTRYH(4)
      REAL :: WA(*), ARG, ARGH, ARGLD, FI, TPI
      DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/3,4,2,5/
      NL = N
      NF = 0
      J = 0
  101 J = J+1
      IF (J.le.4) GO TO 102
      GO TO 103
  102 NTRY = NTRYH(J)
      GO TO 104
  103 NTRY = NTRYH(4)+2*(J-4)
  104 NQ = NL/NTRY
      NR = NL-NTRY*NQ
      IF (NR.eq.0) GO TO 105
      GO TO 101
  105 NF = NF+1
      IFAC(NF+2) = NTRY
      NL = NQ
      IF (NTRY .NE. 2) GO TO 107
      IF (NF .EQ. 1) GO TO 107
      DO 106 I=2,NF
         IB = NF-I+2
         IFAC(IB+2) = IFAC(IB+1)
  106 CONTINUE
      IFAC(3) = 2
  107 IF (NL .NE. 1) GO TO 104
      IFAC(1) = N
      IFAC(2) = NF
      TPI = 6.28318530717959
      ARGH = TPI/FLOAT(N)
      I = 2
      L1 = 1
      DO 110 K1=1,NF
         IP = IFAC(K1+2)
         LD = 0
         L2 = L1*IP
         IDO = N/L2
         IDOT = IDO+IDO+2
         IPM = IP-1
         DO 109 J=1,IPM
            I1 = I
            WA(I-1) = 1.
            WA(I) = 0.
            LD = LD+L1
            FI = 0.
            ARGLD = FLOAT(LD)*ARGH
            DO 108 II=4,IDOT,2
               I = I+2
               FI = FI+1.
               ARG = FI*ARGLD
               WA(I-1) = COS(ARG)
               WA(I) = SIN(ARG)
  108       CONTINUE
            IF (IP .LE. 5) GO TO 109
            WA(I1-1) = WA(I-1)
            WA(I1) = WA(I)
  109    CONTINUE
         L1 = L2
  110 CONTINUE
      RETURN
      END

max_stars_repo_name: magnusja/scipy

max_issues_repo_name: magnusja/scipy

max_forks_repo_name: magnusja/scipy


! { dg-do compile }
! { dg-options "-fcray-pointer" }
!
! Test for PR37039, from an issue on comp.lang.fortran
! http://groups.google.com/group/comp.lang.fortran/msg/8cfa06f222721386

      subroutine test(nnode)
      implicit none
      integer n,nnode
      pointer(ip_tab, tab)
      integer , dimension(1:nnode) :: tab
      do n=1,nnode
         tab(n) = 0
      enddo
      end subroutine test


text:
      program t
      character::c
      data c/1h*/
      print *,c
      end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
!-----------------------------------------------------------------------
!     Module:        fieldlines_init_fline
!     Authors:       S. Lazerson (samuel.lazerson@ipp.mpg.de)
!     Date:          10/30/2020
!     Description:   This subrotine loads a previous fieldlines run and
!                    generates a set of fieldline starting points
!                    from a given field line trace.
!-----------------------------------------------------------------------
      SUBROUTINE fieldlines_init_fline
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stel_kinds, ONLY: rprec
      USE fieldlines_lines, ONLY: nlines
      USE fieldlines_runtime
      USE ez_hdf5
      USE mpi_params
      USE mpi_inc
!-----------------------------------------------------------------------
!     Local Variables
!          ier            Error Flag
!-----------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER :: ier, i
      INTEGER :: npoinc_old, nlines_old
      REAL(rprec) :: dphi_old
      REAL(rprec), ALLOCATABLE, DIMENSION(:,:) :: R_old, Z_old, PHI_old

!      INTEGER, PARAMETER :: line_select = 64
      INTEGER, PARAMETER :: nnew_lines = 2**17
!-----------------------------------------------------------------------
!     Begin Subroutine
!-----------------------------------------------------------------------
      IF (lverb) THEN
         WRITE(6,'(A)')  '----- Generating Starting Points -----'
         WRITE(6,'(A)')  '   FILE: '//TRIM(restart_string)
         WRITE(6,'(A,I6)') '   NEW_LINES: ',nnew_lines*2
         WRITE(6,'(A,I3)') '      LINE #: ',line_select
      END IF
#if defined(MPI_OPT)
      CALL MPI_BARRIER(MPI_COMM_FIELDLINES,ierr_mpi)
#endif
      ! Setup PHI_END before going on
      i = MAXLOC(ABS(PHI_end),DIM=1)
      PHI_end(1:nnew_lines) = PHI_end(i)
      ! Initialize fieldine data
      r_start   = -1
      z_start   = -1
      phi_start =  0
      phi_end(nnew_lines+1:MAXLINES)   =  0
#if defined(LHDF5)
      IF (myworkid == master) THEN
         CALL open_hdf5(TRIM(restart_string),fid,ier,LCREATE=.false.)
         CALL read_scalar_hdf5(fid,'nsteps',ier,INTVAR=npoinc_old)
         CALL read_scalar_hdf5(fid,'nlines',ier,INTVAR=nlines_old)
         IF (lverb) THEN
            WRITE(6,'(A,I3)') '   OLD_LINES: ',nlines_old
            WRITE(6,'(A,I6)') '      NSTEPS: ',npoinc_old
         END IF
         ALLOCATE(R_old(nlines_old,npoinc_old+1), &
                  Z_old(nlines_old,npoinc_old+1), &
                  PHI_old(nlines_old,npoinc_old+1))
         CALL read_var_hdf5(fid,'R_lines',nlines_old,npoinc_old+1,ier,DBLVAR=R_old)
         CALL read_var_hdf5(fid,'Z_lines',nlines_old,npoinc_old+1,ier,DBLVAR=Z_old)
         CALL read_var_hdf5(fid,'PHI_lines',nlines_old,npoinc_old+1,ier,DBLVAR=PHI_old)
         CALL close_hdf5(fid,ier)
         npoinc_old = npoinc_old - 1
         ! Setup new PHI
         dphi_old = PHI_old(line_select,2) - PHI_old(line_select,1)
         CALL RANDOM_NUMBER(PHI_start(1:nnew_lines))
         PHI_start = PHI_start * PHI_old(line_select,npoinc_old)
         CALL spline_it(npoinc_old, &
                        PHI_old(line_select,1:npoinc_old), &
                        R_old(line_select,1:npoinc_old), &
                        nnew_lines, &
                        PHI_start(1:nnew_lines), &
                        R_start(1:nnew_lines), 0)
         CALL spline_it(npoinc_old, &
                        PHI_old(line_select,1:npoinc_old), &
                        Z_old(line_select,1:npoinc_old), &
                        nnew_lines, &
                        PHI_start(1:nnew_lines), &
                        Z_start(1:nnew_lines), 0)
         DEALLOCATE(R_old, PHI_old, Z_old)
         ! Now setup reverse fieldline trace
         i = nnew_lines+1
         PHI_start(i:2*nnew_lines) = PHI_start(1:nnew_lines)
         R_start(i:2*nnew_lines) = R_start(1:nnew_lines)
         Z_start(i:2*nnew_lines) = Z_start(1:nnew_lines)
         PHI_end(i:2*nnew_lines) = -PHI_end(1:nnew_lines)
         nlines = 2*nnew_lines
      END IF
#endif
#if defined(MPI_OPT)
      CALL MPI_BARRIER(MPI_COMM_FIELDLINES,ierr_mpi)
      CALL MPI_BCAST(nlines,           1, MPI_INTEGER, master, MPI_COMM_FIELDLINES, ierr_mpi)
      CALL MPI_BCAST(r_start,   MAXLINES,   MPI_REAL8, master, MPI_COMM_FIELDLINES, ierr_mpi)
      CALL MPI_BCAST(z_start,   MAXLINES,   MPI_REAL8, master, MPI_COMM_FIELDLINES, ierr_mpi)
      CALL MPI_BCAST(phi_start, MAXLINES,   MPI_REAL8, master, MPI_COMM_FIELDLINES, ierr_mpi)
      CALL MPI_BCAST(phi_end,   MAXLINES,   MPI_REAL8, master, MPI_COMM_FIELDLINES, ierr_mpi)
#endif
      RETURN
!-----------------------------------------------------------------------
!     End Subroutine
!-----------------------------------------------------------------------
      END SUBROUTINE fieldlines_init_fline

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
!-----------------------------------------------------------------------
!     Module:        diagno_runtime
!     Authors:       S. Lazerson (lazerson@pppl.gov)
!     Date:          03/02/2012
!     Description:   This module contains various runtime parameters
!                    for the DIAGNO code.  It also contains the
!                    handle_error subroutine.
!
!     Notes:
!       v3.11        Properly handles curtor printout.
!       v3.12        Added BPROBE_TURNS to namelist
!       v3.13        Adjusted ordering of BPROBE output to give accurate
!                    values, adjusted form of rogowski output
!       v3.20        Fixed error in vacuum integrated fields
!       v3.25        Vaccum currents now selectable (luse_extcur)
!                    Output format now ES22.12E3.
!       v3.30        Added -bench benchmarking capability
!                    bprobe_turns increased to 2048 elements for ITER
!       v3.31        Fixed bug in EXTCUR array not reading all elements
!                    in diagno_init_coil.f90
!       v3.32        Fixed bug which caused lfreeb in VMEC to default
!                    to true.
!       v3.40        Modified Rogowski deffinition to handle W7-X style
!                    sensors with multiple unconnected traces in a
!                    given sensor.
!       v3.41        Added LSKIP_FLUX and LSKIP_ROGO to skip
!                    sensor calculations to speed up STELLOPT runs
!                    where not all sensors are used.
!                    Added NFP factor back in to VC
!                    Added EQ_SGNS to handle PHIEDGE sign issues
!       v3.50        Parallelized
!       v3.52        Improved load balancing when skiping calcuations
!                    Increased size of skip variables
!                    Addressed bug regarding short EXTCUR arrays
!       v3.52        Paralle Benchmarks now working
!                    SEGROG_TURNS added
!-----------------------------------------------------------------------
      MODULE diagno_runtime
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stel_kinds, ONLY: rprec
      USE EZspline
!-----------------------------------------------------------------------
!     Module Variables
!          lverb          Logical to control screen output
!          lvmec          Logical to indicate VMEC Equilibria
!          lpies          Logical to indicate PIES Equilibria
!          lspec          Logical to indicate SPEC Equilibria
!          lcoil          Logical to indicate coil file
!          lvac           Logical to indicate vacuum fields only
!          lmut           Logical to control mutual induction calc
!          luse_mut       Logical to control use of mutal induction (not fully implemented)
!          lvc_field      Logical to control use of VIRTUAL Casing
!          nextcur        Number of external currents
!          int_step       Number of sub steps to take for diagnostic integration
!          nu             Number of polodial gridpoints
!          nv             Number of toroidal gridpoints
!          nfp_diagno     Number of field periods
!          flux_turns     Scalar integer multiple for flux loops
!          vc_adapt_tol   Absolute tollerance for adaptive integration
!          vc_adapt_rel   Relative tollerance for adaptive integration
!          units          Scale factor for units
!          phiedge        Total Enclosed toroidal flux
!          id_string      Equilibrium Filename
!          coil_string    Coil Filename
!          flux_diag_file Fluxloop specification file
!          bprobes_file   B-Field diagnostics specification file
!          mirnov_file    Mirnov Array specification file
!          seg_rog_file   Segmented Rogowski Coil specification file
!          int_type       Integration type ('midpoint','simpson','bode')
!          flux_mut_file  Mutual inductance file
!          bfield_points_file  B-Field at a point diagnostic
!          extcur         External currents for MGRID calculation
!----------------------------------------------------------------------
      IMPLICIT NONE

      INTEGER, PARAMETER ::  FILE_OPEN_ERR     = 1
      INTEGER, PARAMETER ::  ALLOC_ERR         = 11
      INTEGER, PARAMETER ::  NAMELIST_READ_ERR = 12
      INTEGER, PARAMETER ::  BAD_INPUT_ERR     = 13
      INTEGER, PARAMETER ::  VMEC_INPUT_ERR    = 2
      INTEGER, PARAMETER ::  VMEC_WOUT_ERR     = 21
      INTEGER, PARAMETER ::  MGRID_ERR         = 22
      INTEGER, PARAMETER ::  EZSPLINE_ERR      = 4
      INTEGER, PARAMETER ::  NETCDF_OPEN_ERR      = 5
      INTEGER, PARAMETER ::  PIES_NETCDF_READ_ERR = 51
      INTEGER, PARAMETER ::  HDF5_ERR       = 6
      INTEGER, PARAMETER ::  D02CJF_ERR = 71
      INTEGER, PARAMETER ::  D05ADF_ERR = 72
      INTEGER, PARAMETER ::  C05AJF_ERR = 73
      INTEGER, PARAMETER :: MPI_ERR = 80
      INTEGER, PARAMETER :: MPI_INIT_ERR = 801
      INTEGER, PARAMETER :: MPI_RANK_ERR = 802
      INTEGER, PARAMETER :: MPI_SIZE_ERR = 803
      INTEGER, PARAMETER :: MPI_BARRIER_ERR = 81
      INTEGER, PARAMETER :: MPI_SEND_ERR = 821
      INTEGER, PARAMETER :: MPI_RECV_ERR = 822
      INTEGER, PARAMETER :: MPI_REDU_ERR = 823
      INTEGER, PARAMETER :: MPI_BCAST_ERR = 83
      INTEGER, PARAMETER :: MPI_FINE_ERR = 89

      INTEGER, PARAMETER ::  MAXLINES   = 256

      LOGICAL         :: lverb, lvmec, lpies, lspec, lcoil, lvac, &
                         lrphiz, lmut, luse_mut, lvc_field
      LOGICAL         :: luse_extcur(512),lskip_flux(2048),lskip_rogo(2048)
      INTEGER         :: nextcur, int_step, nu, nv, nfp_diagno, eq_sgns,&
                         nprocs_diagno, mystart, myend
      REAL(rprec)     :: flux_turns(512), bprobe_turns(2048), segrog_turns(256)
      REAL(rprec)     :: vc_adapt_tol, units, phiedge, vc_adapt_rel
      REAL(rprec), ALLOCATABLE :: extcur(:)
      CHARACTER(256)  :: id_string, flux_diag_file, bprobes_file, &
                         mirnov_file, seg_rog_file, bfield_points_file,&
                         int_type, coil_string, flux_mut_file, rog_mut_file,&
                         mir_mut_file, bprobes_mut_file

      REAL(rprec), PARAMETER :: DIAGNO_VERSION = 4.00
      REAL(rprec), PARAMETER ::      pi = 3.14159265358979312D+00
      REAL(rprec), PARAMETER ::     pi2 = 6.28318530717958623
      REAL(rprec), PARAMETER ::  onerad = 1.74532925199432955E-02
      REAL(rprec), PARAMETER ::     mu0 = 1.25663706143591729E-06
      REAL(rprec), PARAMETER :: sqrtmu0 = 1.12099824327958572E-03
!-----------------------------------------------------------------------
!     Subroutines
!          handle_error  Controls Program Termination
!-----------------------------------------------------------------------
      CONTAINS

      SUBROUTINE handle_error(error_num,string_val,ierr)
      IMPLICIT NONE
      INTEGER,INTENT(in)      :: error_num
      INTEGER,INTENT(inout)      :: ierr
      CHARACTER(*),INTENT(in) :: string_val
      WRITE(6,*) '!!!!! ERROR !!!!!'

      IF (error_num .eq. FILE_OPEN_ERR) THEN
            WRITE(6,*) '  DIAGNO COULD NOT OPEN A FILE.'
            WRITE(6,*) '  FILENAME: ',TRIM(string_val)
            WRITE(6,*) '  IERR:     ',ierr
      ELSEIF (error_num .eq. VMEC_INPUT_ERR) THEN
            WRITE(6,*) '  DIAGNO COULD READ THE VMEC INDATA NAMELIST'
            WRITE(6,*) '  FILENAME: ',TRIM(string_val)
            WRITE(6,*) '  IERR:     ',ierr
      ELSEIF (error_num .eq. VMEC_WOUT_ERR) THEN
            WRITE(6,*) '  DIAGNO COULD READ THE VMEC WOUT FILE'
            WRITE(6,*) '  FILENAME: ',TRIM(string_val)
            WRITE(6,*) '  IERR:     ',ierr
      ELSEIF (error_num .eq. ALLOC_ERR) THEN
            WRITE(6,*) '  DIAGNO ENCOUNTERED AN ALLOCATION ERROR'
            WRITE(6,*) '  VARIABLES: ',TRIM(string_val)
            WRITE(6,*) '  IERR:      ',ierr
      ELSEIF (error_num .eq. EZSPLINE_ERR) THEN
            WRITE(6,*) '  DIAGNO ENCOUNTERED AN EZSPLINE ERROR'
            WRITE(6,*) '  ROUTINE/VAR: ',TRIM(string_val)
            WRITE(6,*) '  IERR:      ',ierr
            CALL EZspline_error(ierr)
      ELSEIF (error_num .eq. MGRID_ERR) THEN
            WRITE(6,*) '  DIAGNO ENCOUNTERED AN ERROR READING THE MGRID FILE'
            WRITE(6,*) '  VARIABLES: ',TRIM(string_val)
            WRITE(6,*) '  IERR:      ',ierr
      ELSEIF (error_num .eq. NETCDF_OPEN_ERR) THEN
            WRITE(6,*) '  DIAGNO ENCOUNTERED AN ERROR OPENING A NETCDF FILE'
            WRITE(6,*) '  VARIABLES: ',TRIM(string_val)
            WRITE(6,*) '  IERR:      ',ierr
      ELSEIF (error_num .eq. NAMELIST_READ_ERR) THEN
            WRITE(6,*) '  DIAGNO ENCOUNTERED AN ERROR READING A NAMELIST'
            WRITE(6,*) '  ',TRIM(string_val)
            WRITE(6,*) '  IERR:      ',ierr
      ELSEIF (error_num .eq. D02CJF_ERR) THEN
            WRITE(6,*) '  DIAGNO ENCOUNTERED A NAG ERROR (D02CJF)'
            WRITE(6,*) '     CALLING FUNCTION',TRIM(string_val)
            WRITE(6,*) '     IERR:      ',ierr
            IF (ierr .eq. 1) THEN
               WRITE(6,*) '     CHECK INPUT PARAMETERS!'
            ELSEIF (ierr .eq. 2) THEN
               WRITE(6,*) '     VALUE OF TOL PREVENTS INTEGRATION!'
            ELSEIF (ierr .eq. 3) THEN
               WRITE(6,*) '     TOL TOO SMALL FOR FIRST STEP!'
            ELSEIF (ierr .eq. 4) THEN
               WRITE(6,*) '     XSOL HAS NOT BEEN RESET OR IS BEHIND X!'
            ELSEIF (ierr .eq. 5) THEN
               WRITE(6,*) '     XSOL HAS NOT BEEN RESET OR IS BEHIND X!'
            ELSEIF (ierr .eq. 6) THEN
               WRITE(6,*) '     ROOT COULD NOT BE FOUND!'
            ELSEIF (ierr .eq. 7) THEN
               WRITE(6,*) '     A SERIOUS NAG ERROR HAS OCCURED!'
            ELSE
               WRITE(6,*) '     UNKNOWN NAG ERROR!'
            END IF
      ELSEIF (error_num .eq. D05ADF_ERR) THEN
            WRITE(6,*) '  DIAGNO ENCOUNTERED A NAG ERROR (D05ADF)'
            WRITE(6,*) '     CALLING FUNCTION',TRIM(string_val)
            WRITE(6,*) '     IERR:      ',ierr
            IF (ierr .eq. 1) THEN
               WRITE(6,*) '     EPS <= 0 or A=B or F(A)xF(B)>0'
            ELSEIF (ierr .eq. 2) THEN
               WRITE(6,*) '     EPS IS TOO SMALL!'
            ELSEIF (ierr .eq. 3) THEN
               WRITE(6,*) '     POLE OF F FOUND!'
            ELSEIF (ierr .eq. 4) THEN
               WRITE(6,*) '     A SERIOUS NAG ERROR HAS OCCURED!'
            ELSE
               WRITE(6,*) '     UNKNOWN NAG ERROR!'
            END IF
      ELSEIF (error_num .eq. C05AJF_ERR) THEN
            WRITE(6,*) '  DIAGNO ENCOUNTERED A NAG ERROR (C05AJF)'
            WRITE(6,*) '     CALLING FUNCTION',TRIM(string_val)
            WRITE(6,*) '     IERR:      ',ierr
            IF (ierr .eq. 1) THEN
               WRITE(6,*) '     EPS <= 0 or NFMAX <= 0'
            ELSEIF (ierr .eq. 2) THEN
               WRITE(6,*) '     BAD SCALE FACTOR!'
            ELSEIF (ierr .eq. 3) THEN
               WRITE(6,*) '     NO ZERO OR ACCURACY TOO HIGH!'
            ELSEIF (ierr .eq. 4) THEN
               WRITE(6,*) '     MORE THAN NFMAX CALLS HAVE BEEN MADE!'
            ELSEIF (ierr .eq. 5) THEN
               WRITE(6,*) '     A SERIOUS NAG ERROR HAS OCCURED!'
            ELSE
               WRITE(6,*) '     UNKNOWN NAG ERROR!'
            END IF
      ELSEIF (error_num .eq. HDF5_ERR) THEN
            WRITE(6,*) '  DIAGNO ENCOUNTERED AN HDF ERROR'
            WRITE(6,*) '  ROUTINE:   ',TRIM(string_val)
            WRITE(6,*) '  IERR:      ',ierr
      ELSE
           WRITE(6,*) '  DIAGNO ENCOUNTERED A GENERIC ERROR'
           WRITE(6,*) '  STRING: ',TRIM(string_val)
           WRITE(6,*) '  ierr:   ',ierr
      END IF
      CALL FLUSH(6)
!DEC$ IF DEFINED (MPI_OPT)
      CALL MPI_FINALIZE(ierr)
!DEC$ ENDIF
      STOP
      END SUBROUTINE handle_error

      END MODULE diagno_runtime


text:
      SUBROUTINE eqsolve(nsval, interp_flag, ier_flag,
     1    lreset, lfirst, lscreen, reset_file_name)
      USE vmec_main
      USE vmec_params, ONLY: ntmax, ns4, bad_init_jac_flag,
     1                       bad_jacobian_flag
      USE realspace
      USE vsvd
      USE xstuff
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER :: nsval, ier_flag
      CHARACTER*(*) :: reset_file_name
      LOGICAL :: interp_flag, lreset, lfirst, lscreen
C-----------------------------------------------
C   L o c a l   P a r a m e t e r s
C-----------------------------------------------
      REAL(rprec), PARAMETER :: p98 = 0.98_dp, p96 = 0.96_dp
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: nsold=0, neqs2_old=0
      REAL(rprec) :: w1, r00s, w0, res0, wdota, r0dot
      REAL(rprec) :: delt0
      LOGICAL :: liter_flag, lreset_internal
C-----------------------------------------------
!
!                INDEX OF LOCAL VARIABLES
!
!        hs      radial mesh size increment
!        iequi   counter used to call -EQFOR- at end of run
!        irzloff offset in xc array between R,Z,L components
!        ijacob  counter for number of times jacobian changes sign
!        irst    counter monitoring sign of jacobian; resets R, Z, and
!                Lambda when jacobian changes sign and decreases time step
!        signgs  sign of Jacobian : must be =1 (right-handed) or =-1 (left-handed)
!        iterj   stores position in main iteration loop (j=1,2)
!        itfsq   counter for storing FSQ into FSQT for plotting
!        ivac    counts number of free-boundary iterations
!        ndamp   number of iterations over which damping is averaged
!        meven   parity selection label for even poloidal modes of R and Z
!        modd    parity selection label for odd poloidal modes of R and
!        gc      stacked array of R, Z, Lambda Spectral force coefficients (see readin for stack order)
!        xc      stacked array of scaled R, Z, Lambda Fourier coefficients

!
!       INITIALIZE MESH-DEPENDENT SCALARS
!
      ns = nsval
      ns1 = ns - 1
      hs = one/ns1
      ohs = one/hs
      mns = ns*mnsize
      irzloff = ntmax*mns
      nrzt = nznt*ns
      neqs = 3*irzloff
      neqs1 = neqs + 1
      neqs2 = neqs1 + 1
      ijacob = 0
      delt0   = delt

      WRITE (nthreed, 10) ns, mnmax, ftolv
      IF (lscreen) PRINT 10, ns, mnmax, ftolv
   10 FORMAT(/'  NS = ',i4,' NO. FOURIER MODES = ',i4,' FTOLV = ',
     1   1p,e10.3)
      IF (imovephi .gt. 0 .and. lscreen) PRINT *,
     1   'Turning on Ip Matching by varying Phi-Edge.'

!
!     ALLOCATE NS-DEPENDENT ARRAYS
!
      CALL allocate_ns(lreset, interp_flag, neqs2_old)

      lreset_internal = (ns .gt. nsold)

      IF (lreset_internal .and. neqs2_old.gt.0) gc(1:neqs2_old) =
     1  scalxc(1:neqs2_old)*xstore(1:neqs2_old)

      lreset_internal = lreset_internal .or. lreset

 1000 CONTINUE

      itfsq = 0
      fsq     = one
      rsfac   = one
      w1      = zero
      r00s    = zero
      gphifac = zero
      grmse   = zero

!
!     COMPUTE INITIAL R, Z AND MAGNETIC FLUX PROFILES
!
   20 CONTINUE
      iter2 = 1
      irst = 1
      CALL profil1d (xc, xcdot, lreset)
      IF (lfirst .and. .not.lreset)               !!Command line with lreset=F
     1    CALL load_xc_from_wout(xc(1), xc(1+irzloff), xc(1+2*irzloff),
     2         ntor, mpol1, ns, reset_file_name, lreset_internal)
      CALL profil3d (xc(1), xc(1+irzloff), lreset_internal)

!
!     INTERPOLATE FROM COARSE TO NEXT FINER RADIAL GRID
!
      IF (interp_flag) CALL interp (xc, gc, scalxc, ns, nsold)
      nsold = ns
      neqs2_old = neqs2

!
!     Store XC,XCDOT for possible restart
!
      CALL restart_iter(delt0)
      iter1 = iter2
      liter_flag = .true.
      ier_flag = 0

!
!     ENTER FORCE ITERATION LOOP
!
      iter_loop: DO WHILE (liter_flag)
!
!     ADVANCE FOURIER AMPLITUDES OF R, Z, AND LAMBDA
!
         CALL evolve (delt0, ier_flag, liter_flag, lscreen, r0dot)
         IF (ijacob.eq.0 .and. ier_flag.eq.bad_init_jac_flag
     1                   .and. ns.ge.3) THEN
            IF (lscreen) PRINT *,
     1         ' TRYING TO IMPROVE INITIAL MAGNETIC AXIS GUESS'
            CALL guess_axis (r1, z1, ru0, zu0)
            lreset_internal = .true.
            ijacob = 1
            GOTO 20
         ELSE IF (ier_flag .ne. 0) THEN
            IF (ier_flag .eq. bad_init_jac_flag)
     1          ier_flag = bad_jacobian_flag
            RETURN
         ENDIF

         w0 = wb + wp/(gamma - one)

!
!     ADDITIONAL STOPPING CRITERION (set liter_flag to FALSE)
!
         IF (ijacob .eq. 25) THEN
            irst = 2
            CALL restart_iter(delt0)
            delt0 = p98*delt
            IF (lscreen) PRINT 120, delt0
            GOTO 1000
         ELSE IF (ijacob .eq. 50) THEN
            irst = 2
            CALL restart_iter(delt0)
            delt0 = p96*delt
            IF (lscreen) PRINT 120, delt0
            GOTO 1000
         ELSE IF (ijacob .ge. 75) THEN
            ier_flag = 2
            liter_flag = .false.
         ELSE IF (iter2.ge.niter .and. liter_flag) THEN
            ier_flag = 4
            liter_flag = .false.
         ENDIF

!
!       TIME STEP CONTROL
!
         IF (iter2 .eq. iter1) res0 = fsq
         res0 = MIN(res0,fsq)
!       Store current state (irst=1)
         IF (fsq .le. res0 .and. iter2-iter1 .gt. 10) THEN
            CALL restart_iter(delt0)
!       Residuals are growing in time, reduce time step
         ELSE IF (fsq .gt. 100.0_dp*res0 .and. iter2 .gt. iter1) THEN
            irst = 2
         ELSE IF (iter2 - iter1 .gt. ns4/2 .and. iter2 .gt. 2*ns4
     1        .and. fsqr+fsqz .gt. c1pm2) THEN
            irst = 3
         ENDIF

         IF (irst .ne. 1) THEN
!       Retrieve previous good state
            CALL restart_iter(delt0)
            iter1 = iter2
         ELSE
!       Increment time step and Printout every nstep iterations
            IF (MOD(iter2,nstep).eq.0 .or. iter2.eq.1 .or.
     1         .not.liter_flag) CALL printout(iter2, delt0, w0, lscreen)
            iter2 = iter2 + 1
         ENDIF

!       Store force residual, wdot for plotting
         wdota = ABS(w0 - w1)/w0
         r0dot = ABS(r00 - r00s)/r00
         r00s = r00
         w1 = w0
         IF (ivac.eq.1 .and. lreset) THEN
            IF (lscreen) PRINT 110, iter2
            WRITE (nthreed, 110) iter2
            ivac = ivac + 1
         ENDIF
!
!       STORE FSQ FOR PLOTTING. EVENTUALLY, STORE FOR EACH RADIAL MESH
!
         IF (MOD(iter2,niter/nstore_seq + 1).eq.0 .and. ns.eq.
     1      ns_array(multi_ns_grid)) THEN
            IF (itfsq .lt. nstore_seq) THEN
              itfsq = itfsq + 1
              fsqt(itfsq) = fsqr + fsqz
              wdot(itfsq) = MAX(wdota,c1pm13)
            END IF
         END IF

      END DO iter_loop

      WRITE (nthreed, 60) wdota, r0dot, w0*twopi**2
      IF (lrecon) WRITE (nthreed, 70) r00*fsqsum0/wb

   60 FORMAT(/,' d(ln W)/dt = ',1p,e14.3,' d(ln R0)/dt = ',e14.3,/,
     1         ' MHD Energy = ',e14.6)
   70 FORMAT(' Average radial force balance: Int[FR(m=0)]',
     1   '/Int(B**2/R) = ',1p,e12.5,' (should tend to zero)'/)
  110 FORMAT(/,2x,'VACUUM PRESSURE TURNED ON AT ',i4,' ITERATIONS'/)
  120 FORMAT(2x,'HAVING A CONVERGENCE PROBLEM: RESETTING DELT TO ',f8.3,
     1  /,2x,'If this does NOT resolve problem, try changing ',
     2       '(decrease OR increase) the value of DELT')

      END SUBROUTINE eqsolve


text:
!+ Data module for organizational variables for the parallel program
!------------------------------------------------------------------------------

MODULE data_parallel

!------------------------------------------------------------------------------
!
! Description:
!  This module contains all variables that are necessary for the organization
!  of the parallel program (number of processors, virtual topology, neighbors,
!  etc.). In a sequential environment some of these variables are also
!  necessary and are set accordingly (e.g. nproc = 1 and my_*_id = 0).
!
! Current Code Owner: DWD, Ulrich Schaettler
!  phone:  +49  69  8062 2739
!  fax:    +49  69  8062 3721
!  email:  ulrich.schaettler@dwd.de
!
! History:
! Version    Date       Name
! ---------- ---------- ----
! 1.1        1998/03/11 Ulrich Schaettler
!  Initial release
! 1.10       1998/09/29 Ulrich Schaettler
!  Changed igroup_diag and icomm_diag to allocatable arrays and
!  eliminated use of parameters for grid point and diagnostic output.
! 1.29       1999/05/11 Ulrich Schaettler
!  Added variables for later use of MPE_IO
! 1.39       2000/05/03 Ulrich Schaettler
!  Added declarations of buffers for distributing NAMELIST Input.
! 2.8        2001/07/06 Ulrich Schaettler
!  Added new variable lreorder for specifying whether the PEs may be reordered
!  while building the cartesian MPI-communicator
! 2.17       2002/05/08 Ulrich Schaettler
!  Added variable imp_grib for specifying the MPI type of variables with
!  KIND parameter irealgrib
! 3.5        2003/09/02 Ulrich Schaettler
!  Added variables ncomm_type, ldatatypes, ltime_barrier for the optimization
!  of the boundary exchange (as Namelist parameteres)
! V4_13        2010/05/11 Michael Gertz
!  Adaptions to SVN
! V4_17        2011/02/24 Ulrich Blahak
!  Eliminated my_peri_neigh
!
! Code Description:
! Language: Fortran 90.
! Software Standards: "European Standards for Writing and
! Documenting Exchangeable Fortran 90 Code".
!==============================================================================
!
! Declarations:
!
! Modules used:
USE data_parameters , ONLY :   &
           ireals,    & ! KIND-type parameters for real variables
           iintegers    ! kind-type parameter for "normal" integer variables

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

IMPLICIT NONE

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

! Global (i.e. public) Declarations:

! 1. Information about the processor grid
! ---------------------------------------

  LOGICAL                           ::           &
    ldatatypes,      & ! if .TRUE.: use MPI-Datatypes for some communications
    ltime_barrier,   & ! if .TRUE.: use additional barriers for determining the
                       ! load-imbalance
    lasync_io          ! if .TRUE.: the model runs with extra PEs for
                       ! asynchronous IO

  INTEGER   (KIND=iintegers)       ::           &
    nprocx,          & ! number of processors in x-direction
    nprocy,          & ! number of processors in y-direction
    nprocio,         & ! number of extra processors for doing asynchronous IO
    nproc,           & ! total number of processors: nprocx * nprocy
    num_compute,     & ! number of compute PEs
    num_io,          & ! number of IO PEs

    nboundlines,     & ! number of boundary lines of the domain for which
                       ! no forecast is computed = overlapping boundary
                       ! lines of the subdomains
    ncomm_type,      & ! type of communication for boundary exchange

    my_world_id,     & ! rank of this subdomain in the global communicator
    my_cart_id,      & ! rank of this subdomain in the cartesian communicator
    my_cart_pos(2),  & ! position of this subdomain in the cartesian grid
                       ! in x- and y-direction
    my_cart_neigh(4)   ! neighbors of this subdomain in the cartesian grid
! UB >>
!!$    my_peri_neigh(4)   ! periodic neighbors of this subdomain in the periodic
!!$                       ! grid if it lies at the boundary of the total domain

  INTEGER   (KIND=iintegers), ALLOCATABLE       ::           &
    isubpos(:,:)       ! positions of the subdomains in the total domain. Given
                       ! are the i- and the j-indices of the lower left and the
                       ! upper right grid point in the order
                       !                  i_ll, j_ll, i_ur, j_ur.
                       ! Only the interior of the domains are considered, not
                       ! the boundary lines.

! 2. Further information for MPI
! ------------------------------

  INTEGER   (KIND=iintegers)       ::           &
    igroup_world,        & ! group that belongs to MPI_COMM_WORLD, i.e. all
                           ! processors
    icomm_world,         & ! communicator that belongs to igroup_world, i.e.
                           ! = MPI_COMM_WORLD
    icomm_compute,       & ! communicator for the group of compute PEs
    igroup_cart,         & ! group of the compute PEs
    icomm_cart,          & ! communicator for the virtual cartesian topology
    icomm_row,           & ! communicator for a east-west row of processors
    iexch_req(4),        & ! stores the sends requests for the neighbor-exchange
                           ! that can be used by MPI_WAIT to identify the send
    imp_reals,           & ! determines the correct REAL type used in the model
                           ! for MPI
    imp_grib,            & ! determines the REAL type used for the GRIB library
    imp_integers,        & ! determines the correct INTEGER type used in the
                           ! model for MPI
    imp_byte,            & ! determines the correct BYTE type used in the model
                           ! for MPI
    imp_character,       & ! determines the correct CHARACTER type used in the
                           ! model for MPI
    imp_logical            ! determines the correct LOGICAL   type used in the
                           ! model for MPI

  LOGICAL                                       &
    lcompute_pe,         & ! indicates whether this is a compute PE or not
    lreorder               ! during the creation of the virtual topology the
                           ! ranking of the processors may be reordered

  INTEGER (KIND=iintegers), ALLOCATABLE  ::           &
    igroup_diag(:),      & ! groups for the diagnostic domains
    icomm_diag (:)         ! communicators for the diagnostic domains

! 3. Static Send buffers
! ----------------------

  REAL (KIND=ireals), ALLOCATABLE  ::           &
    sendbuf(:,:)       ! sending buffer for boundary exchange:
                       ! 1-4 are used for sending, 5-8 are used for receiving

  INTEGER (KIND=iintegers)         ::           &
    isendbuflen        ! length of one column of sendbuf

  ! Buffers for distributing the Namelists
  INTEGER (KIND=iintegers), ALLOCATABLE ::   intbuf  (:)
  REAL    (KIND=ireals)   , ALLOCATABLE ::   realbuf (:)
  LOGICAL                 , ALLOCATABLE ::   logbuf  (:)
  CHARACTER (LEN=100)     , ALLOCATABLE ::   charbuf (:)

!==============================================================================
! input file suffix
  integer :: cosmo_input_suffix

END MODULE data_parallel

max_stars_repo_name: niklaswr/TSMP

max_issues_repo_name: niklaswr/TSMP

max_forks_repo_name: niklaswr/TSMP


text:
!======================================================================!
  subroutine GloSum(phi)
!----------------------------------------------------------------------!
!   Estimates global summ among all processors.                        !
!----------------------------------------------------------------------!
  implicit none
!------------------------------[Include]-------------------------------!
  include 'mpif.h'
!-----------------------------[Arguments]------------------------------!
  real    :: phi
!-------------------------------[Locals]-------------------------------!
  real    :: phi_new
  integer :: error
!======================================================================!

!================================================
      call MPI_ALLREDUCE      &
!-----------------------------------+------------
             (phi,            & ! send buffer
              phi_new,        & ! recv buffer
              1,              & ! length
              MPI_REAL,       & ! datatype
              MPI_SUM,        & ! operation
              MPI_COMM_WORLD, &
              error)
!================================================

  phi = phi_new

  end subroutine

max_stars_repo_name: mhadziabdic/T-FlowS

max_issues_repo_name: mhadziabdic/T-FlowS

max_forks_repo_name: mhadziabdic/T-FlowS


           program demo_isblank
           use M_strings, only : isblank
           implicit none
           integer                    :: i
           character(len=1),parameter :: string(*)=[(char(i),i=0,127)]
              write(*,'(*(g0,1x))')'ISXBLANK: ',&
              & iachar(pack( string, isblank(string) ))
           end program demo_isblank


! { dg-do compile }
!
! PR fortran/34657
!
module test_mod
interface
  subroutine my_sub (a)
    real a
  end subroutine
end interface
end module

subroutine my_sub (a)
  use test_mod, gugu => my_sub
  real a
  print *, a
end subroutine

END


text:
! Copyright (c) 2018, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

module m
implicit none
integer :: i = 0
contains
  subroutine p()
    i = 1
  end subroutine p

  subroutine foo(fun_ptr)
    procedure(p), pointer, intent(out) :: fun_ptr
    fun_ptr => p
  end subroutine
end module m

program test
use m
implicit none
procedure(), pointer :: x
call foo(x)
call x()
if (i == 1) then
  write(*,*) 'PASS'
else
  write(*,*) 'FAIL'
  STOP 1
end if
end program test

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
! RUN: bbc %s -o - | FileCheck %s

! CHECK-LABEL: fir.global @_QMc_interoperability_testEthis_thing : !fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}> {
! CHECK:         %[[VAL_0:.*]] = arith.constant 0 : i64
! CHECK:         %[[VAL_1:.*]] = fir.undefined !fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}>
! CHECK:         %[[VAL_2:.*]] = fir.undefined !fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>
! CHECK:         %[[VAL_3:.*]] = fir.insert_value %[[VAL_2]], %[[VAL_0]], ["__address", !fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>] : (!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>, i64) -> !fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>
! CHECK:         %[[VAL_4:.*]] = fir.insert_value %[[VAL_1]], %[[VAL_3]], ["cptr", !fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}>] : (!fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}>, !fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>) -> !fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}>
! CHECK:         fir.has_value %[[VAL_4]] : !fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}>
! CHECK:       }

! CHECK-LABEL: func @_QMc_interoperability_testPget_a_thing() -> !fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}> {
! CHECK:         %[[VAL_2:.*]] = fir.alloca !fir.box<!fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}>>
! CHECK:         %[[VAL_3:.*]] = fir.address_of(@_QMc_interoperability_testEthis_thing) : !fir.ref<!fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}>>
! CHECK:         %[[VAL_4:.*]] = fir.alloca !fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}> {bindc_name = "get_a_thing", uniq_name = "_QMc_interoperability_testFget_a_thingEget_a_thing"}
! CHECK:         %[[VAL_5:.*]] = fir.embox %[[VAL_4]] : (!fir.ref<!fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}>>) -> !fir.box<!fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}>>
! CHECK:         %[[VAL_10:.*]] = fir.embox %[[VAL_3]] : (!fir.ref<!fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}>>) -> !fir.box<!fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}>>
! CHECK:         return %{{.*}} : !fir.type<_QMc_interoperability_testTthing_with_pointer{cptr:!fir.type<_QM__fortran_builtinsT__builtin_c_ptr{__address:i64}>}>
! CHECK:       }

module c_interoperability_test
  use iso_c_binding, only: c_ptr, c_null_ptr

  type thing_with_pointer
     type(c_ptr) :: cptr = c_null_ptr
  end type thing_with_pointer

  type(thing_with_pointer) :: this_thing

contains
  function get_a_thing()
    type(thing_with_pointer) :: get_a_thing
    get_a_thing = this_thing
  end function get_a_thing
end module c_interoperability_test

max_stars_repo_name: akashhansda/llvm-project

max_issues_repo_name: akashhansda/llvm-project

max_forks_repo_name: akashhansda/llvm-project