text:
!** Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
!** See https://llvm.org/LICENSE.txt for license information.
!** SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception

! check whether inlining with optional arguments works
! if the argument is not type scalar integer
! this test: inline with derived type array arguments
module m
 type dt
  integer :: m1
  real :: m2
  integer, dimension(3) :: m3
 end type
 contains

 subroutine s1( a, b, c )
  type(dt) :: a(3)
  integer :: b
  type(dt),optional :: c(2)

  a(1)%m1 = b+1
  a(1)%m2 = b*2
  a(1)%m3(1) = b+2
  a(1)%m3(2) = b+3
  a(1)%m3(3) = b+4
  a(2)%m1 = b+11
  a(2)%m2 = b*2 + 1
  a(2)%m3(1) = b+12
  a(2)%m3(2) = b+13
  a(2)%m3(3) = b+14
  a(3)%m1 = b+21
  a(3)%m2 = b*2 + 2
  a(3)%m3(1) = b+22
  a(3)%m3(2) = b+23
  a(3)%m3(3) = b+24
  if( present(c) ) then
   c(1)%m1 = b+5
   c(1)%m2 = b*3
   c(1)%m3(1) = b+6
   c(1)%m3(2) = b+7
   c(1)%m3(3) = b+8
   c(2)%m1 = b+15
   c(2)%m2 = b*3 + 8
   c(2)%m3(1) = b+16
   c(2)%m3(2) = b+17
   c(2)%m3(3) = b+18
  endif
 end subroutine
end module

program p
 use m
 type(dt) a(3),c(2)
 integer b
 integer result(40)
 integer expect(40)
 data expect /101,200,102,103,104,111,201,112,113,114,121,202, &
              122,123,124,105,300,106,107,108,115,308,116,117, &
              118,201,400,202,203,204,211,401,212,213,214,221, &
              402,222,223,224/

 b = 100
 call s1( a, b, c )
 result( 1) = a(1)%m1
 result( 2) = a(1)%m2
 result( 3) = a(1)%m3(1)
 result( 4) = a(1)%m3(2)
 result( 5) = a(1)%m3(3)
 result( 6) = a(2)%m1
 result( 7) = a(2)%m2
 result( 8) = a(2)%m3(1)
 result( 9) = a(2)%m3(2)
 result(10) = a(2)%m3(3)
 result(11) = a(3)%m1
 result(12) = a(3)%m2
 result(13) = a(3)%m3(1)
 result(14) = a(3)%m3(2)
 result(15) = a(3)%m3(3)
 result(16) = c(1)%m1
 result(17) = c(1)%m2
 result(18) = c(1)%m3(1)
 result(19) = c(1)%m3(2)
 result(20) = c(1)%m3(3)
 result(21) = c(2)%m1
 result(22) = c(2)%m2
 result(23) = c(2)%m3(1)
 result(24) = c(2)%m3(2)
 result(25) = c(2)%m3(3)
 call s1( a, 200 )
 result(26) = a(1)%m1
 result(27) = a(1)%m2
 result(28) = a(1)%m3(1)
 result(29) = a(1)%m3(2)
 result(30) = a(1)%m3(3)
 result(31) = a(2)%m1
 result(32) = a(2)%m2
 result(33) = a(2)%m3(1)
 result(34) = a(2)%m3(2)
 result(35) = a(2)%m3(3)
 result(36) = a(3)%m1
 result(37) = a(3)%m2
 result(38) = a(3)%m3(1)
 result(39) = a(3)%m3(2)
 result(40) = a(3)%m3(3)
 !print *,result
 call check(result,expect,40)
end program

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
      SUBROUTINE SBORN(P,ANS_SUMMED)
C
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
C     By the MadGraph5_aMC@NLO Development Team
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
C
C
C     Return the sum of the split orders which are required in
C      orders.inc (BORN_ORDERS)
C     Also the values needed for the counterterms are stored in the
C      C_BORN_CNT common block
C
C
C     Process: u u~ > t t~ [ real = QCD QED ] QCD^2<=6 QED^2<=0
C     Process: c c~ > t t~ [ real = QCD QED ] QCD^2<=6 QED^2<=0
C     Process: d d~ > t t~ [ real = QCD QED ] QCD^2<=6 QED^2<=0
C     Process: s s~ > t t~ [ real = QCD QED ] QCD^2<=6 QED^2<=0
C
C
C     CONSTANTS
C
      IMPLICIT NONE
      INCLUDE 'nexternal.inc'
      INTEGER NAMPSO, NSQAMPSO
      PARAMETER (NAMPSO=1, NSQAMPSO=1)
      INTEGER NGRAPHS
      PARAMETER (NGRAPHS=   1)
C
C     ARGUMENTS
C
      REAL*8 P(0:3,NEXTERNAL-1)
C
C     VARIABLES
C
      INTEGER I,J,K
      INCLUDE 'orders.inc'
      DOUBLE PRECISION ANS_SUMMED
      COMPLEX*16 ANS(2,0:NSQAMPSO), ANS_CNT(2, NSPLITORDERS)
      LOGICAL KEEP_ORDER(NSQAMPSO), KEEP_ORDER_CNT(NSPLITORDERS,
     $  NSQAMPSO), FIRSTTIME
      DATA KEEP_ORDER / NSQAMPSO * .TRUE. /
      COMMON /C_KEEP_ORDER_CNT/ KEEP_ORDER_CNT
      COMMON /C_BORN_CNT/ ANS_CNT
      INTEGER ORD_SUBTRACT
      DATA FIRSTTIME / .TRUE. /
      INTEGER AMP_ORDERS(NSPLITORDERS)
      DOUBLE PRECISION TINY
      PARAMETER (TINY = 1D-12)
      DOUBLE PRECISION MAX_VAL
      DOUBLE PRECISION       WGT_ME_BORN,WGT_ME_REAL
      COMMON /C_WGT_ME_TREE/ WGT_ME_BORN,WGT_ME_REAL

      DOUBLE PRECISION AMP2B(1), JAMP2B(0:2,0:NAMPSO)
      COMMON/TO_AMPS_BORN/  AMP2B,       JAMP2B
      DOUBLE PRECISION AMP2(1), JAMP2(0:2)
      COMMON/TO_AMPS/  AMP2,       JAMP2
      LOGICAL SPLIT_TYPE_USED(NSPLITORDERS)
      COMMON/TO_SPLIT_TYPE_USED/SPLIT_TYPE_USED
C
C     FUNCTIONS
C
      INTEGER GETORDPOWFROMINDEX_B
      INTEGER ORDERS_TO_AMP_SPLIT_POS
C
C     BEGIN CODE
C
C     look for orders which match the born order constraint

      IF (FIRSTTIME) THEN
        DO I = 1, NSQAMPSO
C         this is for the orders of the born to integrate
          DO J = 1, NSPLITORDERS
            IF(GETORDPOWFROMINDEX_B(J, I) .GT. BORN_ORDERS(J)) THEN
              KEEP_ORDER(I) = .FALSE.
              EXIT
            ENDIF
          ENDDO

C         this is for the orders of the counterterms
          DO J = 1, NSPLITORDERS
            KEEP_ORDER_CNT(J,I) = .TRUE.
            DO K = 1, NSPLITORDERS
              IF (J.EQ.K) THEN
                ORD_SUBTRACT=2
              ELSE
                ORD_SUBTRACT=0
              ENDIF
              IF(GETORDPOWFROMINDEX_B(K, I) .GT. NLO_ORDERS(K)
     $         -ORD_SUBTRACT) THEN
                KEEP_ORDER_CNT(J,I) = .FALSE.
                EXIT
              ENDIF
            ENDDO
          ENDDO
          IF (KEEP_ORDER(I)) THEN
            WRITE(*,*) 'BORN: keeping split order ', I
          ELSE
            WRITE(*,*) 'BORN: not keeping split order ', I
          ENDIF
        ENDDO

        DO J = 1, NSPLITORDERS
          WRITE(*,*) 'counterterm S.O', J, ORDERNAMES(J)
          DO I = 1, NSQAMPSO
            IF (KEEP_ORDER_CNT(J,I)) THEN
              WRITE(*,*) 'BORN: keeping split order', I
            ELSE
              WRITE(*,*) 'BORN: not keeping split order', I
            ENDIF
          ENDDO
        ENDDO
        FIRSTTIME = .FALSE.
      ENDIF

      CALL SBORN_SPLITORDERS(P,ANS)

C     the born to be integrated
      ANS_SUMMED = 0D0
      MAX_VAL = 0D0

C     reset the amp_split array
      AMP_SPLIT(1:AMP_SPLIT_SIZE) = 0D0
      AMP_SPLIT_CNT(1:AMP_SPLIT_SIZE,1:2,1:NSPLITORDERS) = DCMPLX(0D0
     $ ,0D0)

      DO I = 1, NSQAMPSO
        MAX_VAL = MAX(MAX_VAL, ABS(ANS(1,I)))
      ENDDO
      DO I = 1, NSQAMPSO
        IF (KEEP_ORDER(I)) THEN
          ANS_SUMMED = ANS_SUMMED + ANS(1,I)
C         keep track of the separate pieces correspoinding to
C          different coupling combinations
          DO J = 1, NSPLITORDERS
            AMP_ORDERS(J) = GETORDPOWFROMINDEX_B(J, I)
          ENDDO
          IF(ABS(ANS(1,I)).GT.MAX_VAL*TINY)
     $      AMP_SPLIT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS)) = ANS(1,I)
        ENDIF
      ENDDO
C     this is to avoid fake non-zero contributions
      IF (ABS(ANS_SUMMED).LT.MAX_VAL*TINY) ANS_SUMMED=0D0

      WGT_ME_BORN=ANS_SUMMED

C     fill the amp2 and jamp2 arrays
      AMP2(1:NGRAPHS)=AMP2B(1:NGRAPHS)  ! amp2 just needs to be copyed

      DO I = 0, INT(JAMP2B(0,0))
        JAMP2(I)=0D0
        DO J = 1, NAMPSO
C         here sum all, this may be refined later
          JAMP2(I)=JAMP2(I)+JAMP2B(I,J)
        ENDDO
      ENDDO

C     quantities for the counterterms
      DO J = 1, NSPLITORDERS
        ANS_CNT(1:2,J) = (0D0, 0D0)
        IF (.NOT.SPLIT_TYPE_USED(J)) CYCLE
        DO I = 1, NSQAMPSO
          IF (KEEP_ORDER_CNT(J,I)) THEN
            ANS_CNT(1,J) = ANS_CNT(1,J) + ANS(1,I)
            ANS_CNT(2,J) = ANS_CNT(2,J) + ANS(2,I)
C           keep track of the separate pieces also for counterterms
            DO K = 1, NSPLITORDERS
              AMP_ORDERS(K) = GETORDPOWFROMINDEX_B(K, I)
C             take into account the fact that the counterterm for a
C              given split order
C             will be multiplied by the corresponding squared coupling
              IF (K.EQ.J) AMP_ORDERS(K) = AMP_ORDERS(K) + 2
            ENDDO
C           this is to avoid fake non-zero contributions
            IF (ABS(ANS(1,I)).GT.MAX_VAL*TINY)
     $        AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1,J) =
     $        ANS(1,I)
            IF (ABS(ANS(2,I)).GT.MAX_VAL*TINY)
     $        AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),2,J) =
     $        ANS(2,I)
          ENDIF
        ENDDO
C       this is to avoid fake non-zero contributions
        IF (ABS(ANS_CNT(1,J)).LT.MAX_VAL*TINY) ANS_CNT(1,J)=(0D0,0D0)
        IF (ABS(ANS_CNT(2,J)).LT.MAX_VAL*TINY) ANS_CNT(2,J)=(0D0,0D0)
      ENDDO

 999  RETURN
      END

      SUBROUTINE SBORN_SPLITORDERS(P1,ANS)
C
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
C     By the MadGraph5_aMC@NLO Development Team
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
C
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
C     AND HELICITIES
C     FOR THE POINT IN PHASE SPACE P1(0:3,NEXTERNAL-1)
C
C     Process: u u~ > t t~ [ real = QCD QED ] QCD^2<=6 QED^2<=0
C     Process: c c~ > t t~ [ real = QCD QED ] QCD^2<=6 QED^2<=0
C     Process: d d~ > t t~ [ real = QCD QED ] QCD^2<=6 QED^2<=0
C     Process: s s~ > t t~ [ real = QCD QED ] QCD^2<=6 QED^2<=0
C
      IMPLICIT NONE
C
C     CONSTANTS
C
      INCLUDE 'nexternal.inc'
      INCLUDE 'born_nhel.inc'
      INTEGER     NCOMB
      PARAMETER ( NCOMB=  16 )
      INTEGER NAMPSO, NSQAMPSO
      PARAMETER (NAMPSO=1, NSQAMPSO=1)
      INTEGER    THEL
      PARAMETER (THEL=NCOMB*6)
      INTEGER NGRAPHS
      PARAMETER (NGRAPHS=   1)
C
C     ARGUMENTS
C
      REAL*8 P1(0:3,NEXTERNAL-1)
      COMPLEX*16 ANS(2,0:NSQAMPSO)
C
C     LOCAL VARIABLES
C
      INTEGER IHEL,IDEN,I,J,JJ,GLU_IJ
      REAL*8 BORNS(2,0:NSQAMPSO)
      COMPLEX*16 BORNTILDE
      INTEGER NTRY(6)
      DATA NTRY /6*0/
      COMPLEX*16 T(2,NSQAMPSO)
      INTEGER NHEL(NEXTERNAL-1,NCOMB)
      DATA (NHEL(I,   1),I=1,4) / 1,-1,-1, 1/
      DATA (NHEL(I,   2),I=1,4) / 1,-1,-1,-1/
      DATA (NHEL(I,   3),I=1,4) / 1,-1, 1, 1/
      DATA (NHEL(I,   4),I=1,4) / 1,-1, 1,-1/
      DATA (NHEL(I,   5),I=1,4) / 1, 1,-1, 1/
      DATA (NHEL(I,   6),I=1,4) / 1, 1,-1,-1/
      DATA (NHEL(I,   7),I=1,4) / 1, 1, 1, 1/
      DATA (NHEL(I,   8),I=1,4) / 1, 1, 1,-1/
      DATA (NHEL(I,   9),I=1,4) /-1,-1,-1, 1/
      DATA (NHEL(I,  10),I=1,4) /-1,-1,-1,-1/
      DATA (NHEL(I,  11),I=1,4) /-1,-1, 1, 1/
      DATA (NHEL(I,  12),I=1,4) /-1,-1, 1,-1/
      DATA (NHEL(I,  13),I=1,4) /-1, 1,-1, 1/
      DATA (NHEL(I,  14),I=1,4) /-1, 1,-1,-1/
      DATA (NHEL(I,  15),I=1,4) /-1, 1, 1, 1/
      DATA (NHEL(I,  16),I=1,4) /-1, 1, 1,-1/
      INTEGER IDEN_VALUES(6)
      DATA IDEN_VALUES /36, 36, 36, 36, 36, 36/
      INTEGER IJ_VALUES(6)
      DATA IJ_VALUES /1, 2, 0, 0, 1, 2/
C
C     GLOBAL VARIABLES
C
      DOUBLE PRECISION AMP2(NGRAPHS), JAMP2(0:2,0:NAMPSO)
      COMMON/TO_AMPS_BORN/  AMP2,       JAMP2
      DATA JAMP2(0,0) /   2/
      LOGICAL GOODHEL(NCOMB,6)
      COMMON /C_GOODHEL/GOODHEL
      DOUBLE COMPLEX SAVEAMP(NGRAPHS,MAX_BHEL)
      COMMON/TO_SAVEAMP/SAVEAMP
      DOUBLE PRECISION SAVEMOM(NEXTERNAL-1,2)
      COMMON/TO_SAVEMOM/SAVEMOM
      DOUBLE PRECISION HEL_FAC
      INTEGER GET_HEL,SKIP(6)
      COMMON/CBORN/HEL_FAC,GET_HEL,SKIP
      LOGICAL CALCULATEDBORN
      COMMON/CCALCULATEDBORN/CALCULATEDBORN
      INTEGER NFKSPROCESS
      COMMON/C_NFKSPROCESS/NFKSPROCESS
      LOGICAL COND_IJ
C     ----------
C     BEGIN CODE
C     ----------
      IDEN=IDEN_VALUES(NFKSPROCESS)
      GLU_IJ = IJ_VALUES(NFKSPROCESS)
      NTRY(NFKSPROCESS)=NTRY(NFKSPROCESS)+1
      IF (NTRY(NFKSPROCESS).LT.2) THEN
        IF (GLU_IJ.EQ.0) THEN
          SKIP(NFKSPROCESS)=0
        ELSE
          SKIP(NFKSPROCESS)=1
          DO WHILE(NHEL(GLU_IJ ,SKIP(NFKSPROCESS)).NE.-NHEL(GLU_IJ ,1))
            SKIP(NFKSPROCESS)=SKIP(NFKSPROCESS)+1
          ENDDO
          SKIP(NFKSPROCESS)=SKIP(NFKSPROCESS)-1
        ENDIF
      ENDIF
      DO JJ=1,NGRAPHS
        AMP2(JJ)=0D0
      ENDDO
      DO I=0,NAMPSO
        DO JJ=1,INT(JAMP2(0,0))
          JAMP2(JJ,I)=0D0
        ENDDO
      ENDDO
      IF (CALCULATEDBORN) THEN
        DO J=1,NEXTERNAL-1
          IF (SAVEMOM(J,1).NE.P1(0,J) .OR. SAVEMOM(J,2).NE.P1(3,J))
     $      THEN
            CALCULATEDBORN=.FALSE.
            WRITE (*,*) 'momenta not the same in Born'
            STOP
          ENDIF
        ENDDO
      ENDIF
      IF (.NOT.CALCULATEDBORN) THEN
        DO J=1,NEXTERNAL-1
          SAVEMOM(J,1)=P1(0,J)
          SAVEMOM(J,2)=P1(3,J)
        ENDDO
        DO J=1,MAX_BHEL
          DO JJ=1,NGRAPHS
            SAVEAMP(JJ,J)=(0D0,0D0)
          ENDDO
        ENDDO
      ENDIF
      DO I=0,NSQAMPSO
        ANS(1,I) = 0D0
        ANS(2,I) = 0D0
      ENDDO
      HEL_FAC=1D0
      DO IHEL=1,NCOMB
          ! the following lines are to avoid segfaults when glu_ij=0
        COND_IJ=SKIP(NFKSPROCESS).EQ.0
        IF (.NOT.COND_IJ) COND_IJ=COND_IJ.OR.NHEL(GLU_IJ,IHEL)
     $   .EQ.NHEL(GLU_IJ,1)
          !if (nhel(glu_ij,ihel).EQ.NHEL(GLU_IJ,1).or.skip(nfksprocess).eq.0) then
        IF (COND_IJ) THEN
          IF ((GOODHEL(IHEL,NFKSPROCESS) .OR. GOODHEL(IHEL
     $     +SKIP(NFKSPROCESS),NFKSPROCESS) .OR. NTRY(NFKSPROCESS) .LT.
     $      2) ) THEN

            CALL BORN(P1,NHEL(1,IHEL),IHEL,T,BORNS)
            DO I=1,NSQAMPSO
              ANS(1,I)=ANS(1,I)+T(1,I)
              ANS(2,I)=ANS(2,I)+T(2,I)
            ENDDO
            IF ( BORNS(1,0).NE.0D0 .AND. .NOT. GOODHEL(IHEL
     $       ,NFKSPROCESS) ) THEN
              GOODHEL(IHEL,NFKSPROCESS)=.TRUE.
            ENDIF
            IF ( BORNS(2,0).NE.0D0 .AND. .NOT. GOODHEL(IHEL
     $       +SKIP(NFKSPROCESS),NFKSPROCESS) ) THEN
              GOODHEL(IHEL+SKIP(NFKSPROCESS),NFKSPROCESS)=.TRUE.
            ENDIF
          ENDIF
        ENDIF
      ENDDO
      DO I=1,NSQAMPSO
        ANS(1,I)=ANS(1,I)/DBLE(IDEN)
        ANS(2,I)=ANS(2,I)/DBLE(IDEN)
        ANS(1,0)=ANS(1,0)+ANS(1,I)
        ANS(2,0)=ANS(2,0)+ANS(2,I)
      ENDDO
      CALCULATEDBORN=.TRUE.
      END

      SUBROUTINE BORN(P,NHEL,HELL,ANS,BORNS)
C
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
C     By the MadGraph5_aMC@NLO Development Team
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
C     FOR THE POINT WITH EXTERNAL LINES W(0:6,NEXTERNAL-1)

C     Process: u u~ > t t~ [ real = QCD QED ] QCD^2<=6 QED^2<=0
C     Process: c c~ > t t~ [ real = QCD QED ] QCD^2<=6 QED^2<=0
C     Process: d d~ > t t~ [ real = QCD QED ] QCD^2<=6 QED^2<=0
C     Process: s s~ > t t~ [ real = QCD QED ] QCD^2<=6 QED^2<=0
C
      IMPLICIT NONE
C
C     CONSTANTS
C
      INTEGER NAMPSO, NSQAMPSO
      PARAMETER (NAMPSO=1, NSQAMPSO=1)
      INTEGER    NGRAPHS,    NEIGEN
      PARAMETER (NGRAPHS=   1,NEIGEN=  1)
      INTEGER    NWAVEFUNCS, NCOLOR
      PARAMETER (NWAVEFUNCS=5, NCOLOR=2)
      REAL*8     ZERO
      PARAMETER (ZERO=0D0)
      COMPLEX*16 IMAG1
      PARAMETER (IMAG1 = (0D0,1D0))
      INCLUDE 'nexternal.inc'
      INCLUDE 'born_nhel.inc'
      INCLUDE 'coupl.inc'
C
C     ARGUMENTS
C
      REAL*8 P(0:3,NEXTERNAL-1),BORNS(2,0:NSQAMPSO)
      INTEGER NHEL(NEXTERNAL-1), HELL
      COMPLEX*16 ANS(2,NSQAMPSO)
C
C     LOCAL VARIABLES
C
      INTEGER I,J,M,N,IHEL,BACK_HEL,GLU_IJ
      INTEGER IC(NEXTERNAL-1),NMO
      PARAMETER (NMO=NEXTERNAL-1)
      DATA IC /NMO*1/
      REAL*8 CF(NCOLOR,NCOLOR)
      COMPLEX*16 ZTEMP, AMP(NGRAPHS), JAMP(NCOLOR,NAMPSO), W(8
     $ ,NWAVEFUNCS), JAMPH(2, NCOLOR,NAMPSO)
      COMPLEX*16 TMP_JAMP(0)
C
C     GLOBAL VARIABLES
C
      DOUBLE PRECISION AMP2(1), JAMP2(0:2,0:NAMPSO)
      COMMON/TO_AMPS_BORN/  AMP2,       JAMP2
      DOUBLE COMPLEX SAVEAMP(NGRAPHS,MAX_BHEL)
      COMMON/TO_SAVEAMP/SAVEAMP
      DOUBLE PRECISION HEL_FAC
      INTEGER GET_HEL,SKIP(6)
      COMMON/CBORN/HEL_FAC,GET_HEL,SKIP
      LOGICAL CALCULATEDBORN
      COMMON/CCALCULATEDBORN/CALCULATEDBORN
      INTEGER NFKSPROCESS
      COMMON/C_NFKSPROCESS/NFKSPROCESS
      INTEGER STEP_HEL
      LOGICAL COND_IJ

C
C     FUNCTION
C
      INTEGER SQSOINDEXB

      INTEGER IJ_VALUES(6)
      DATA IJ_VALUES /1, 2, 0, 0, 1, 2/
C
C     COLOR DATA
C
      DATA (CF(I,  1),I=  1,  2) /9.000000000000000D+00
     $ ,3.000000000000000D+00/
C     1 T(2,1) T(3,4)
      DATA (CF(I,  2),I=  1,  2) /3.000000000000000D+00
     $ ,9.000000000000000D+00/
C     1 T(2,4) T(3,1)
C     ----------
C     BEGIN CODE
C     ----------
      GLU_IJ = IJ_VALUES(NFKSPROCESS)
      DO I = 1, NSQAMPSO
        ANS(1,I)=0D0
        ANS(2,I)=0D0
        BORNS(1,I)=0D0
        BORNS(2,I)=0D0
      ENDDO
      BORNS(1,0)=0D0
      BORNS(2,0)=0D0
      IF (GLU_IJ.NE.0) THEN
        BACK_HEL = NHEL(GLU_IJ)
        IF (BACK_HEL.NE.0) THEN
          STEP_HEL=-2*BACK_HEL
        ELSE
          STEP_HEL=1
        ENDIF
      ELSE
        BACK_HEL=0
        STEP_HEL=1
      ENDIF
      DO IHEL=BACK_HEL,-BACK_HEL,STEP_HEL
        IF (GLU_IJ.NE.0) THEN
          COND_IJ=IHEL.EQ.BACK_HEL.OR.NHEL(GLU_IJ).NE.0
        ELSE
          COND_IJ=IHEL.EQ.BACK_HEL
        ENDIF
        IF (COND_IJ) THEN
          IF (GLU_IJ.NE.0) THEN
            IF (NHEL(GLU_IJ).NE.0) NHEL(GLU_IJ) = IHEL
          ENDIF
          IF (.NOT. CALCULATEDBORN) THEN
            CALL IXXXXX(P(0,1),ZERO,NHEL(1),+1*IC(1),W(1,1))
            CALL OXXXXX(P(0,2),ZERO,NHEL(2),-1*IC(2),W(1,2))
            CALL OXXXXX(P(0,3),MDL_MT,NHEL(3),+1*IC(3),W(1,3))
            CALL IXXXXX(P(0,4),MDL_MT,NHEL(4),-1*IC(4),W(1,4))
            CALL FFV1P0_3(W(1,1),W(1,2),GC_11,ZERO,ZERO,W(1,5))
C           Amplitude(s) for diagram number 1
            CALL FFV1_0(W(1,4),W(1,3),W(1,5),GC_11,AMP(1))
            DO I=1,NGRAPHS
              IF(IHEL.EQ.BACK_HEL)THEN
                SAVEAMP(I,HELL)=AMP(I)
              ELSEIF(IHEL.EQ.-BACK_HEL)THEN
                SAVEAMP(I,HELL+SKIP(NFKSPROCESS))=AMP(I)
              ELSE
                WRITE(*,*) 'ERROR #1 in born.f'
                STOP
              ENDIF
            ENDDO
          ELSEIF (CALCULATEDBORN) THEN
            DO I=1,NGRAPHS
              IF(IHEL.EQ.BACK_HEL)THEN
                AMP(I)=SAVEAMP(I,HELL)
              ELSEIF(IHEL.EQ.-BACK_HEL)THEN
                AMP(I)=SAVEAMP(I,HELL+SKIP(NFKSPROCESS))
              ELSE
                WRITE(*,*) 'ERROR #1 in born.f'
                STOP
              ENDIF
            ENDDO
          ENDIF
C         JAMPs contributing to orders QCD=2 QED=0
          JAMP(1,1) = (1.666666666666667D-01)*AMP(1)
          JAMP(2,1) = (-5.000000000000000D-01)*AMP(1)
          DO M = 1, NAMPSO
            DO I = 1, NCOLOR
              ZTEMP = (0.D0,0.D0)
              DO J = 1, NCOLOR
                ZTEMP = ZTEMP + CF(J,I)*JAMP(J,M)
              ENDDO
              DO N = 1, NAMPSO
                BORNS(2-(1+BACK_HEL*IHEL)/2,SQSOINDEXB(M,N))=BORNS(2
     $           -(1+BACK_HEL*IHEL)/2,SQSOINDEXB(M,N))+ZTEMP
     $           *DCONJG(JAMP(I,N))
              ENDDO
            ENDDO
          ENDDO
          DO I = 1, NGRAPHS
            AMP2(I)=AMP2(I)+AMP(I)*DCONJG(AMP(I))
          ENDDO
          DO J = 1,NAMPSO
            DO I = 1, NCOLOR
              JAMP2(I,J)=JAMP2(I,J)+JAMP(I,J)*DCONJG(JAMP(I,J))
              JAMPH(2-(1+BACK_HEL*IHEL)/2,I,J)=JAMP(I,J)
            ENDDO
          ENDDO
        ENDIF
      ENDDO
      DO I = 1, NSQAMPSO
        BORNS(1,0)=BORNS(1,0)+BORNS(1,I)
        BORNS(2,0)=BORNS(2,0)+BORNS(2,I)
        ANS(1,I) = BORNS(1,I) + BORNS(2,I)
      ENDDO
      DO M = 1, NAMPSO
        DO I = 1, NCOLOR
          ZTEMP = (0.D0,0.D0)
          DO J = 1, NCOLOR
            ZTEMP = ZTEMP + CF(J,I)*JAMPH(2,J,M)
          ENDDO
          DO N = 1, NAMPSO
            ANS(2,SQSOINDEXB(M,N))= ANS(2,SQSOINDEXB(M,N)) + ZTEMP
     $       *DCONJG(JAMPH(1,I,N))
          ENDDO
        ENDDO
      ENDDO
      IF (GLU_IJ.NE.0) NHEL(GLU_IJ) = BACK_HEL
      END

      BLOCK DATA GOODHELS
      INTEGER     NCOMB
      PARAMETER ( NCOMB=  16 )
      INTEGER    THEL
      PARAMETER (THEL=NCOMB*6)
      LOGICAL GOODHEL(NCOMB,6)
      COMMON /C_GOODHEL/GOODHEL
      DATA GOODHEL/THEL*.FALSE./
      END

C
C     Helper functions to deal with the split orders.
C

      INTEGER FUNCTION SQSOINDEXB(AMPORDERA,AMPORDERB)
C
C     This functions plays the role of the interference matrix. It can
C      be hardcoded or
C     made more elegant using hashtables if its execution speed ever
C      becomes a relevant
C     factor. From two split order indices of the jamps, it return the
C      corresponding
C     index in the squared order canonical ordering.
C
C     CONSTANTS
C
      IMPLICIT NONE
      INTEGER NAMPSO, NSQAMPSO
      PARAMETER (NAMPSO=1, NSQAMPSO=1)
      INTEGER NSPLITORDERS
      PARAMETER (NSPLITORDERS=2)
C
C     ARGUMENTS
C
      INTEGER AMPORDERA, AMPORDERB
C
C     LOCAL VARIABLES
C
      INTEGER I, SQORDERS(NSPLITORDERS)
      INTEGER AMPSPLITORDERS(NAMPSO,NSPLITORDERS)
      DATA (AMPSPLITORDERS(  1,I),I=  1,  2) /    2,    0/
C
C     FUNCTION
C
      INTEGER SQSOINDEXB_FROM_ORDERS
C
C     BEGIN CODE
C
      DO I=1,NSPLITORDERS
        SQORDERS(I)=AMPSPLITORDERS(AMPORDERA,I)
     $   +AMPSPLITORDERS(AMPORDERB,I)
      ENDDO
      SQSOINDEXB=SQSOINDEXB_FROM_ORDERS(SQORDERS)
      END

      INTEGER FUNCTION SQSOINDEXB_FROM_ORDERS(ORDERS)
C
C     From a list of values for the split orders, this function
C      returns the
C     corresponding index in the squared orders canonical ordering.
C
      IMPLICIT NONE
      INTEGER NSQAMPSO
      PARAMETER (NSQAMPSO=1)
      INTEGER NSPLITORDERS
      PARAMETER (NSPLITORDERS=2)
C
C     ARGUMENTS
C
      INTEGER ORDERS(NSPLITORDERS)
C
C     LOCAL VARIABLES
C
      INTEGER I,J
      INTEGER SQSPLITORDERS(NSQAMPSO,NSPLITORDERS)
C     the values listed below are for QCD, QED
      DATA (SQSPLITORDERS(  1,I),I=  1,  2) /    4,    0/
C
C     BEGIN CODE
C
      DO I=1,NSQAMPSO
        DO J=1,NSPLITORDERS
          IF (ORDERS(J).NE.SQSPLITORDERS(I,J)) GOTO 1009
        ENDDO
        SQSOINDEXB_FROM_ORDERS = I
        RETURN
 1009   CONTINUE
      ENDDO

      WRITE(*,*) 'ERROR:: Stopping function sqsoindex_from_orders'
      WRITE(*,*) 'Could not find squared orders ',(ORDERS(I),I=1
     $ ,NSPLITORDERS)
      STOP

      END

      INTEGER FUNCTION GETORDPOWFROMINDEX_B(IORDER, INDX)
C
C     Return the power of the IORDER-th order appearing at position
C      INDX
C     in the split-orders output
C
      IMPLICIT NONE
      INTEGER NSQAMPSO
      PARAMETER (NSQAMPSO=1)
      INTEGER NSPLITORDERS
      PARAMETER (NSPLITORDERS=2)
C
C     ARGUMENTS
C
      INTEGER IORDER, INDX
C
C     LOCAL VARIABLES
C
      INTEGER I
      INTEGER SQSPLITORDERS(NSQAMPSO,NSPLITORDERS)
C     the values listed below are for QCD, QED
      DATA (SQSPLITORDERS(  1,I),I=  1,  2) /    4,    0/
C
C     BEGIN CODE
C
      IF (IORDER.GT.NSPLITORDERS.OR.IORDER.LT.1) THEN
        WRITE(*,*) 'INVALID IORDER B', IORDER
        WRITE(*,*) 'SHOULD BE BETWEEN 1 AND ', NSPLITORDERS
        STOP
      ENDIF

      IF (INDX.GT.NSQAMPSO.OR.INDX.LT.1) THEN
        WRITE(*,*) 'INVALID INDX B', INDX
        WRITE(*,*) 'SHOULD BE BETWEEN 1 AND ', NSQAMPSO
        STOP
      ENDIF

      GETORDPOWFROMINDEX_B=SQSPLITORDERS(INDX, IORDER)
      END

      SUBROUTINE GET_NSQSO_B(NSQSO)
C
C     Simple subroutine returning the number of squared split order
C     contributions returned in ANS when calling SMATRIX_SPLITORDERS
C
      IMPLICIT NONE
      INTEGER NSQAMPSO
      PARAMETER (NSQAMPSO=1)
      INTEGER NSQSO

      NSQSO=NSQAMPSO

      END

ext: f


text:
c-----------------------------------------------------------------------
      subroutine local_solves_fdm(u,v)
c
c     Given an input vector v, this returns the additive Schwarz solution
c
c                       -1
c                  u = M   v
c
c                      T     -1
c     where M = sum ( R_i A_i    R_i )
c                i
c
c     The R_i's are simply index_set restriction operators.
c
c     The local solves are performed with the fast diagonalization method.
c
      include 'SIZE'
      include 'INPUT'
      include 'DOMAIN'
      include 'PARALLEL'
c
      include 'TSTEP'
      include 'CTIMER'
c
      real u(lx2,ly2,lz2,lelv),v(lx2,ly2,lz2,lelv)
      common /scrpre/ v1(lx1,ly1,lz1,lelv)
     $               ,w1(lx1,ly1,lz1),w2(lx1,ly1,lz1)
      common /scrover/ ar(lelv)

      parameter(lxx=lx1*lx1, levb=lelv+lbelv)
      common /fastd/  df(lx1*ly1*lz1,levb)
     $             ,  sr(lxx*2,levb),ss(lxx*2,levb),st(lxx*2,levb)
      integer e,eb,eoff

c
      if (icalld.eq.0) tsolv=0.0
      icalld=icalld+1
      nsolv=icalld
c
      ntot1 = lx1*ly1*lz1*nelv
      ntot2 = lx2*ly2*lz2*nelv
c
c     Fill interiors
      iz1 = 0
      if (if3d) iz1=1
      call rzero(v1,ntot1)
      do e=1,nelv
         do iz=1,lz2
         do iy=1,ly2
         do ix=1,lx2
            v1(ix+1,iy+1,iz+iz1,e) = v(ix,iy,iz,e)
         enddo
         enddo
         enddo
      enddo
      call dface_ext    (v1)
      call dssum        (v1,lx1,ly1,lz1)
      call dface_add1si (v1,-1.)
c
c     Now solve each subdomain problem:
c
      etime1=dnekclock()

      eoff  = 0
      if (ifield.gt.1) eoff  = nelv

      do e = 1,nelv
         eb = e + eoff
         call fastdm1(v1(1,1,1,e),df(1,eb)
     $                           ,sr(1,eb),ss(1,eb),st(1,eb),w1,w2)
      enddo
      tsolv=tsolv+dnekclock()-etime1
c
c     Exchange/add elemental solutions
c
      call s_face_to_int (v1,-1.)
      call dssum         (v1,lx1,ly1,lz1)
      call s_face_to_int (v1, 1.)
      if(param(42).eq.0) call do_weight_op(v1)
c
c     Map back to pressure grid (extract interior values)
c
      do e=1,nelv
         do iz=1,lz2
         do iy=1,ly2
         do ix=1,lx2
            u(ix,iy,iz,e) = v1(ix+1,iy+1,iz+iz1,e)
         enddo
         enddo
         enddo
      enddo
c
      return
      end
c-----------------------------------------------------------------------
      subroutine fastdm1(r,df,sr,ss,st,w1,w2)
c
c     Fast diagonalization solver for FEM on mesh 1
c
      include 'SIZE'
      parameter (lxx=lx1*lx1,lxyz=lx1*ly1*lz1)

      real r(1),df(1),sr(lxx,2),ss(lxx,2),st(lxx,2),w1(1),w2(1)
c
c
c      T
c     S  r
      call tensr3 (w1,lx1,r ,lx1,sr(1,2),ss(1,1),st(1,1),w2)
c
c
c      -1 T
c     D  S  r
c
      call col2   (w1,df,lxyz)
c
c
c        -1 T
c     S D  S  r
c
      call tensr3 (r ,lx1,w1,lx1,sr(1,1),ss(1,2),st(1,2),w2)
c
      return
      end
c-----------------------------------------------------------------------
      subroutine tensr3(v,nv,u,nu,A,Bt,Ct,w)
C
C     -  Tensor product application of v = (C x B x A) u
C        NOTE -- the transpose of B & C must be input, rather than B & C.
C
C     -  scratch arrays: w(nu*nu*nv)
C
C
      include 'SIZE'
      include 'INPUT'
      real v(nv,nv,nv),u(nu,nu,nu)
      real A(1),Bt(1),Ct(1)
      real w(1)

      if (nu.gt.nv) then
         write(6,*) nid,nu,nv,' ERROR in tensr3. Contact P.Fischer.'
         write(6,*) nid,nu,nv,' Memory problem.'
         call exitt
      endif

      if (if3d) then
         nuv = nu*nv
         nvv = nv*nv
         call mxm(A,nv,u,nu,v,nu*nu)
         k=1
         l=1
         do iz=1,nu
            call mxm(v(k,1,1),nv,Bt,nu,w(l),nv)
            k=k+nuv
            l=l+nvv
         enddo
         call mxm(w,nvv,Ct,nu,v,nv)
      else
         call mxm(A,nv,u,nu,w,nu)
         call mxm(w,nv,Bt,nu,v,nv)
      endif
      return
      end
c-----------------------------------------------------------------------
      subroutine s_face_to_int(x,c)
c
c     Scale face and add to interior of element
c
      include 'SIZE'
      include 'INPUT'
      real x(lx1,ly1,lz1,1)
c
      do ie=1,nelv
c
        if (if3d) then
c
          do iz=2,lz1-1
          do ix=2,lx1-1
            x(ix,2    ,iz,ie) = c*x(ix,1  ,iz,ie) + x(ix,2    ,iz,ie)
            x(ix,ly1-1,iz,ie) = c*x(ix,ly1,iz,ie) + x(ix,ly1-1,iz,ie)
          enddo
          enddo
c
          do iz=2,lz1-1
          do iy=2,ly1-1
            x(2    ,iy,iz,ie) = c*x(1  ,iy,iz,ie) + x(2    ,iy,iz,ie)
            x(lx1-1,iy,iz,ie) = c*x(lx1,iy,iz,ie) + x(lx1-1,iy,iz,ie)
          enddo
          enddo
c
          do iy=2,ly1-1
          do ix=2,lx1-1
            x(ix,iy,2    ,ie) = c*x(ix,iy,1  ,ie) + x(ix,iy,2    ,ie)
            x(ix,iy,lz1-1,ie) = c*x(ix,iy,lz1,ie) + x(ix,iy,lz1-1,ie)
          enddo
          enddo
c
        else
c         2D
          do ix=2,lx1-1
            x(ix,2    ,1,ie) = c*x(ix,1  ,1,ie) + x(ix,2    ,1,ie)
            x(ix,ly1-1,1,ie) = c*x(ix,ly1,1,ie) + x(ix,ly1-1,1,ie)
          enddo
          do iy=2,ly1-1
            x(2    ,iy,1,ie) = c*x(1  ,iy,1,ie) + x(2    ,iy,1,ie)
            x(lx1-1,iy,1,ie) = c*x(lx1,iy,1,ie) + x(lx1-1,iy,1,ie)
          enddo
        endif
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine dface_ext(x)
c     Extend interior to face of element
c
      include 'SIZE'
      include 'INPUT'
      real x(lx1,ly1,lz1,1)
c
      do ie=1,nelv
c
         if (if3d) then
c
          do iz=2,lz1-1
          do ix=2,lx1-1
            x(ix,1  ,iz,ie) = x(ix,2    ,iz,ie)
            x(ix,ly1,iz,ie) = x(ix,ly1-1,iz,ie)
          enddo
          enddo
c
          do iz=2,lz1-1
          do iy=2,ly1-1
            x(1  ,iy,iz,ie) = x(2    ,iy,iz,ie)
            x(lx1,iy,iz,ie) = x(lx1-1,iy,iz,ie)
          enddo
          enddo
c
          do iy=2,ly1-1
          do ix=2,lx1-1
            x(ix,iy,1  ,ie) = x(ix,iy,2    ,ie)
            x(ix,iy,lz1,ie) = x(ix,iy,lz1-1,ie)
          enddo
          enddo
c
        else
c
          do ix=2,lx1-1
            x(ix,1  ,1,ie) = x(ix,2    ,1,ie)
            x(ix,ly1,1,ie) = x(ix,ly1-1,1,ie)
          enddo
          do iy=2,ly1-1
            x(1  ,iy,1,ie) = x(2    ,iy,1,ie)
            x(lx1,iy,1,ie) = x(lx1-1,iy,1,ie)
          enddo
c
        endif
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine dface_add1si(x,c)
c     Scale interior and add to face of element
c
      include 'SIZE'
      include 'INPUT'
      real x(lx1,ly1,lz1,1)
c
      do ie=1,nelv
c
        if (if3d) then
c
          do iz=2,lz1-1
          do ix=2,lx1-1
            x(ix,1  ,iz,ie) = x(ix,1  ,iz,ie) + c*x(ix,2    ,iz,ie)
            x(ix,ly1,iz,ie) = x(ix,ly1,iz,ie) + c*x(ix,ly1-1,iz,ie)
          enddo
          enddo
c
          do iz=2,lz1-1
          do iy=2,ly1-1
            x(1  ,iy,iz,ie) = x(1  ,iy,iz,ie) + c*x(2    ,iy,iz,ie)
            x(lx1,iy,iz,ie) = x(lx1,iy,iz,ie) + c*x(lx1-1,iy,iz,ie)
          enddo
          enddo
c
          do iy=2,ly1-1
          do ix=2,lx1-1
            x(ix,iy,1  ,ie) = x(ix,iy,1  ,ie) + c*x(ix,iy,2    ,ie)
            x(ix,iy,lz1,ie) = x(ix,iy,lz1,ie) + c*x(ix,iy,lz1-1,ie)
          enddo
          enddo
c
        else
c
c         2D
c
          do ix=2,lx1-1
            x(ix,1  ,1,ie) = x(ix,1  ,1,ie) + c*x(ix,2    ,1,ie)
            x(ix,ly1,1,ie) = x(ix,ly1,1,ie) + c*x(ix,ly1-1,1,ie)
          enddo
          do iy=2,ly1-1
            x(1  ,iy,1,ie) = x(1  ,iy,1,ie) + c*x(2    ,iy,1,ie)
            x(lx1,iy,1,ie) = x(lx1,iy,1,ie) + c*x(lx1-1,iy,1,ie)
          enddo
c
        endif
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine init_weight_op

      include 'SIZE'
      include 'INPUT'
      include 'TSTEP'
      parameter(levb=lelv+lbelv)
      common /swaplengths/ l(lx1,ly1,lz1,lelv)
      common /weightop/ w(lx2,lz2,2,3,levb)
      real l
      real w
      integer e,e0,eb

      e0  = 0
      if (ifield.gt.1) e0 = nelv

      n=lx2+1
      if (if3d) then
         do e=1,nelv
            call rzero(l(1,1,1,e),lx1*ly1*lz1)
            do k=2,n
               do j=2,n
                  l(2,j,k,e)=1
                  l(n,j,k,e)=1
               enddo
            enddo
            do k=2,n
               do i=2,n
                  l(i,2,k,e)=1
                  l(i,n,k,e)=1
               enddo
            enddo
            do j=2,n
               do i=2,n
                  l(i,j,2,e)=1
                  l(i,j,n,e)=1
               enddo
            enddo
         enddo
      else
         do e=1,nelv
            call rzero(l(1,1,1,e),lx1*ly1*lz1)
            do j=2,n
               l(2,j,1,e)=1
               l(n,j,1,e)=1
            enddo
            do i=2,n
               l(i,2,1,e)=1
               l(i,n,1,e)=1
            enddo
         enddo
      endif

      call dface_ext(l)
      call dssum(l,lx1,ly1,lz1)
      call dface_add1si(l,-1.)
      call s_face_to_int(l,1.)
c     l now holds the count matrix C on the outer pressure nodes
      if (if3d) then
         do e=1,nelv
            eb = e0+e
            do k=1,lz2
               do j=1,ly2
                  w(j,k,1,1,eb)=1.0/l(2,j+1,k+1,e)
                  w(j,k,2,1,eb)=1.0/l(n,j+1,k+1,e)
               enddo
            enddo
            do k=1,lz2
               do i=1,lx2
                  w(i,k,1,2,eb)=1.0/l(i+1,2,k+1,e)
                  w(i,k,2,2,eb)=1.0/l(i+1,n,k+1,e)
               enddo
            enddo
            do j=1,ly2
               do i=1,lx2
                  w(i,j,1,3,eb)=1.0/l(i+1,j+1,2,e)
                  w(i,j,2,3,eb)=1.0/l(i+1,j+1,n,e)
               enddo
            enddo
         enddo
      else
         do e=1,nelv
            eb = e0+e
            do j=1,ly2
               w(j,1,1,1,eb)=1.0/l(2,j+1,1,e)
               w(j,1,2,1,eb)=1.0/l(n,j+1,1,e)
            enddo
            do i=1,lx2
               w(i,1,1,2,eb)=1.0/l(i+1,2,1,e)
               w(i,1,2,2,eb)=1.0/l(i+1,n,1,e)
            enddo
         enddo
      endif
      end
c-----------------------------------------------------------------------
      subroutine do_weight_op(x)
      include 'SIZE'
      include 'INPUT'
      include 'TSTEP'
      parameter(levb=lelv+lbelv)
      common /weightop/ w(lx2,lz2,2,3,levb)
      real w

      real x(0:lx1-1,0:ly1-1,0:lz1-1,1)
      integer e,e0,eb

      e0  = 0
      if (ifield.gt.1) e0 = nelv

      if (if3d) then
         do e=1,nelv
            eb = e0 + e
            do k=1,lz2
               do j=1,ly2
                  x(  1,j,k,e)=w(j,k,1,1,eb)*x(  1,j,k,e)
                  x(lx2,j,k,e)=w(j,k,2,1,eb)*x(lx2,j,k,e)
               enddo
            enddo
            do k=1,lz2
               do i=2,lx2-1
                  x(i,  1,k,e)=w(i,k,1,2,eb)*x(i,  1,k,e)
                  x(i,ly2,k,e)=w(i,k,2,2,eb)*x(i,ly2,k,e)
               enddo
            enddo
            do j=2,ly2-1
               do i=2,lx2-1
                  x(i,j,  1,e)=w(i,j,1,3,eb)*x(i,j,  1,e)
                  x(i,j,lz2,e)=w(i,j,2,3,eb)*x(i,j,lz2,e)
               enddo
            enddo
         enddo
      else
         do e=1,nelv
            eb = e0 + e
            do j=1,ly2
               x(  1,j,0,e)=w(j,1,1,1,eb)*x(  1,j,0,e)
               x(lx2,j,0,e)=w(j,1,2,1,eb)*x(lx2,j,0,e)
            enddo
            do i=2,lx2-1
               x(i,  1,0,e)=w(i,1,1,2,eb)*x(i,  1,0,e)
               x(i,ly2,0,e)=w(i,1,2,2,eb)*x(i,ly2,0,e)
            enddo
         enddo
      endif
      end

max_stars_repo_name: RonRahaman/nekRS

max_issues_repo_name: neams-th-coe/nekRS

max_forks_repo_name: neams-th-coe/nekRS


text:
!
!    Copyright 2011 Sebastian Heimann
!
!    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 differentiation

    implicit none

    private

    public :: differentiate

    contains

    subroutine differentiate( dt, f, df, method )

        real, intent(in)                     :: dt
        real, dimension(:), intent(in)       :: f
        real, dimension(:), intent(out)      :: df
        integer, optional, intent(in)        :: method

    ! this subroutine differentiates a sampled function f given at
    ! equitistantly spaced points with the spacing dt and
    ! puts the result to df.

    ! the optional argument method is currently ignored.

    ! the differentiation is currently done by doing central differences:
    !        f'(t) = (f(t+dt)-f(t-dt)) / (2 dt)
    ! boundary points are done using forward or backward differences

        integer :: i,j,il,iu,jl,ju

        if (present(method)) return

        if (size(f) /= size(df)) then
            stop "sizes of input and output array for differentiation do not match."
        end if

        if (size(f) < 2) then
            stop "sizes of arrays for differentiation are too short."
        end if

        jl = lbound(df,1)
        ju = ubound(df,1)
        il = lbound(f,1)
        iu = ubound(f,1)

        j=jl+1
        do i=il+1,iu-1
            df(j) = ( f(i+1) - f(i-1) ) / (dt*2)
            j = j+1
        end do

        df(jl) = ( f(il+1)-f(il) ) / dt
        df(ju) = ( f(iu)-f(iu-1) ) / dt

    end subroutine differentiate

end module differentiation

max_stars_repo_name: emolch/kiwi

max_issues_repo_name: alirezaniki/kiwi

max_forks_repo_name: emolch/kiwi


text:
      PROGRAM driver

      INCLUDE 'KPP_ROOT_Parameters.h'
      INCLUDE 'KPP_ROOT_Global.h'

      KPP_REAL DVAL(NSPEC)
      INTEGER i

C ---- TIME VARIABLES ------------------

      TSTART = 0
      TEND = TSTART + 600
      DT = 60.
      TEMP = 298

      STEPMIN = 0.01
      STEPMAX = 900

      RTOLS = 1e-8
      do i=1,NVAR
        RTOL(i) = RTOLS
        ATOL(i) = 1e-3
      end do

      CALL Initialize()

C ********** TIME LOOP *************************

      CALL InitSaveData()

      write(6,990) (SPC_NAMES[MONITOR(i)], i=1,NMONITOR),
     *             (SMASS(i), i=1,NMASS )
990   FORMAT('done[%] Time[h] ',20(4X,A6))

      TIME = TSTART
      do while (TIME .lt. TEND)

        CALL GetMass( C, DVAL )
        write(6,991) (TIME-TSTART)/(TEND-TSTART)*100, TIME/3600.,
     *               (C(MONITOR(i))/CFACTOR, i=1,NMONITOR),
     *               (DVAL(i)/CFACTOR, i=1,NMASS)
991   FORMAT(F6.1,'% ',F7.2,3X,20(E10.4,2X))

        CALL SaveData()

        CALL Update_SUN()
        CALL Update_RCONST()

        CALL INTEGRATE( TIME, TIME+DT )

      end do

      CALL GetMass( C, DVAL )
      write(6,991) (TIME-TSTART)/(TEND-TSTART)*100, TIME/3600.,
     *               (C(MONITOR(i))/CFACTOR, i=1,NMONITOR),
     *               (DVAL(i)/CFACTOR, i=1,NMASS)

      CALL SaveData()

C *********** END TIME LOOP ********

      CALL CloseSaveData()
      CALL GenerateMatlab(' ')

      STOP
      END

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


text:
      SUBROUTINE AUX(YP, Y, X, PARAM)                                   ABSH0206
C***********************************************************************ABSH0207
C*****AUX EVALUATES THE DERIVATIVE OF Y(X) AND TRANSFORMS THE SYSTEM OF*ABSH0208
C*****SECOND ORDER ODES INTO A SYSTEM OF FIRST ORDER ODES - AL*YP=BL.  *ABSH0209
C*****REFERENCES:                                                      *ABSH0210
C*****L.L.LAO,S.P.HIRSHMAN,R.M.WIELAND,ORNL/TM-7616 (1981).            *ABSH0211
C*****                                ,PHYS FLUIDS 24 (1981) AUG       *ABSH0212
C*****L.L.LAO,R.M.WIELAND,W.A.HOULBERG,S.P.HIRSHMAN,ORNL/TM-7871 (1981)*ABSH0213
C*****LAST REVISION: 6/81 L.L.LAO, R.M.WIELAND, AND W.A.HOULBERG ORNL. *ABSH0214
C*****CALCULATED PARAMETERS:                                           *ABSH0215
C*****YP(I)-DERIVATIVE OF Y(I) WITH RESPECT TO X.                      *ABSH0216
C*****INPUT PARAMETERS:                                                *ABSH0217
C*****Y(I)-SOLUTION OF FIRST ORDER ODE I AT X.                         *ABSH0218
C*****X-ABSCISSA AT WHICH YP(I) ARE TO BE EVALUATED.                   *ABSH0219
C*****OTHER COMMENTS:                                                  *ABSH0220
C*****AUX IS REQUIRED BY D02AGF ODE SOLVER IN NAG LIBRARY.             *ABSH0221
C***********************************************************************ABSH0222
      PARAMETER MNMX=61                                                 ABSH0223
      PARAMETER MTHETA=11,MM1=10,MNMOMS=3                               ABSH0224
      PARAMETER MNEQ=6                                                  ABSH0225
      COMMON /CFEVAL/ NFL                                               ABSH0226
      COMMON /LOCMEQ/ XL0, DX, R0, A0, BT0, P0, AI0, E0, E1,            ABSH0227
     *     D0, ARC, SLPE, SLPR0, BETAP0, KMHD, IUNDER                   ABSH0228
      COMMON /NIO/ NIN, NOUT, NPLOT, NDEBUG, NTTY                       ABSH0229
      COMMON /SCFUNC/ S1TH(MTHETA), S2TH(MTHETA), C1TH(MTHETA),         ABSH0230
     *     C2TH(MTHETA)                                                 ABSH0231
      COMMON /SWORK2/ BSIFT(MM1), CSIFT(MM1), DSIFT(MM1),               ABSH0232
     *     BELLI(MM1), CELLI(MM1), DELLI(MM1), BTRIA(MM1),              ABSH0233
     *     CTRIA(MM1), DTRIA(MM1), R(MTHETA), RX(MTHETA),               ABSH0234
     *     RT(MTHETA), RXT(MTHETA), RTT(MTHETA), Z(MTHETA),             ABSH0235
     *     ZX(MTHETA), ZT(MTHETA), ZXT(MTHETA), ZTT(MTHETA),            ABSH0236
     *     GTT(MTHETA), TAU(MTHETA), GSQRT(MTHETA),                     ABSH0237
     *     GTTX(MTHETA), ZMOM1(MTHETA), RMOM0(MTHETA),                  ABSH0238
     *     TMOM(MTHETA), RMOM2(MTHETA), RTAU2(MTHETA),                  ABSH0239
     *     RTAU3(MTHETA), GS10(MTHETA), GS21(MTHETA),                   ABSH0240
     *     GS4(MTHETA), GS5(MTHETA), GS6(MTHETA), GS7(MTHETA),          ABSH0241
     *     GS8(MTHETA), GS2(MTHETA), GS12(MTHETA),                      ABSH0242
     *     GS22(MTHETA), CS13N(MTHETA), SS20(MTHETA),                   ABSH0243
     *     CS1N(MTHETA)                                                 ABSH0244
      COMMON /WEIGHT/ WFCN(MTHETA), CNORM                               ABSH0245
      DIMENSION YP(MNEQ), Y(MNEQ), PARAM(MNEQ)                          ABSH0246
      DATA NFL /0/                                                      ABSH0247
      NFL = NFL + 1                                                     ABSH0248
      NTHETA = MTHETA                                                   ABSH0249
      CALL PRSFUN(X, PX, PPX)                                           ABSH0250
      CALL CURFUN(X, AIX, AIPX)                                         ABSH0251
      AIPOI = AIPX/AIX                                                  ABSH0252
      AI = BETAP0/2.0*PPX/AIX**2                                        ABSH0253
C*****VECTORIZE THE COEFFICIENTS IN THE MOMENT EQUATIONS.               ABSH0254
      DO 10 I=1,NTHETA                                                  ABSH0255
           R(I) = Y(1) - X*C1TH(I) + Y(5)*C2TH(I)                       ABSH0256
           RX(I) = Y(2) - C1TH(I) + Y(6)*C2TH(I)                        ABSH0257
           RT(I) = X*S1TH(I) - 2.0*Y(5)*S2TH(I)                         ABSH0258
           RXT(I) = S1TH(I) - 2.0*Y(6)*S2TH(I)                          ABSH0259
           RTT(I) = X*C1TH(I) - 4.0*Y(5)*C2TH(I)                        ABSH0260
           Z(I) = Y(3)*(X*S1TH(I)+Y(5)*S2TH(I))                         ABSH0261
           ZX(I) = Y(4)*(X*S1TH(I)+Y(5)*S2TH(I)) +                      ABSH0262
     *          Y(3)*(S1TH(I)+Y(6)*S2TH(I))                             ABSH0263
           ZT(I) = Y(3)*(X*C1TH(I)+2.0*Y(5)*C2TH(I))                    ABSH0264
           ZXT(I) = Y(4)*(X*C1TH(I)+2.0*Y(5)*C2TH(I)) +                 ABSH0265
     *          Y(3)*(C1TH(I)+2.0*Y(6)*C2TH(I))                         ABSH0266
           ZTT(I) = -Y(3)*(X*S1TH(I)+4.0*Y(5)*S2TH(I))                  ABSH0267
           GTT(I) = RT(I)**2 + ZT(I)**2                                 ABSH0268
           TAU(I) = RT(I)*ZX(I) - RX(I)*ZT(I)                           ABSH0269
           GSQRT(I) = R(I)*TAU(I)                                       ABSH0270
           RTAU2(I) = GSQRT(I)*TAU(I)                                   ABSH0271
           RTAU3(I) = RTAU2(I)*TAU(I)                                   ABSH0272
           GS10(I) = ZT(I)*GTT(I)/RTAU3(I)                              ABSH0273
           GS12(I) = GS10(I)*C2TH(I)                                    ABSH0274
           GS2(I) = RT(I)*GTT(I)/RTAU3(I)                               ABSH0275
           GS21(I) = GS2(I)*S1TH(I)                                     ABSH0276
           GS22(I) = GS2(I)*S2TH(I)                                     ABSH0277
           GS4(I) = 1.0/R(I)                                            ABSH0278
           GS5(I) = ZT(I)/R(I)/GSQRT(I)                                 ABSH0279
           CS13N(I) = GTT(I)*(RX(I)*ZXT(I)-ZX(I)*RXT(I))                ABSH0280
           GTTX(I) = RT(I)*RXT(I) + ZT(I)*ZXT(I)                        ABSH0281
           GS6(I) = (CS13N(I)+(RX(I)*RT(I)+ZX(I)*ZT(I))*(RTT(I)*        ABSH0282
     *          ZX(I)+RT(I)*ZXT(I)-RXT(I)*ZT(I)-RX(I)*ZTT(I)))/         ABSH0283
     *          RTAU3(I)                                                ABSH0284
           GS7(I) = (GTTX(I)-RX(I)*RTT(I)-ZX(I)*ZTT(I))/RTAU2(I)        ABSH0285
           GS8(I) = GTT(I)/RTAU2(I)                                     ABSH0286
           SS20(I) = GS5(I) + GS6(I) + GS7(I)                           ABSH0287
           CS1N(I) = (2.0*GTTX(I)-GTT(I)*RX(I)/R(I)+CS13N(I)/           ABSH0288
     *          TAU(I))/GSQRT(I)                                        ABSH0289
           RMOM0(I) = ZT(I)*WFCN(I)                                     ABSH0290
           TMOM(I) = TAU(I)*WFCN(I)                                     ABSH0291
           RMOM2(I) = (Y(3)*RT(I)*S2TH(I)-ZT(I)*C2TH(I))*WFCN(I)        ABSH0292
C          ZMOM1(I) = (-X*S1TH(I)+Y(5)*S2TH(I))*RT(I)*WFCN(I)           ABSH0293
           ZMOM1(I) = (+X*S1TH(I)+Y(5)*S2TH(I))*RT(I)*WFCN(I)           ABSH002A
   10 CONTINUE                                                          ABSH0294
      YP(1) = Y(2)                                                      ABSH0295
      YP(3) = Y(4)                                                      ABSH0296
      YP(5) = Y(6)                                                      ABSH0297
C*****R0 MOMENT TERMS - SHIFT.                                          ABSH0298
      CR010 = SDOT(NTHETA,RMOM0,1,GS10,1)                               ABSH0299
      CR012 = SDOT(NTHETA,RMOM0,1,GS12,1)                               ABSH0300
      CR021 = SDOT(NTHETA,RMOM0,1,GS21,1)                               ABSH0301
      CR022 = SDOT(NTHETA,RMOM0,1,GS22,1)                               ABSH0302
      CR03 = SDOT(NTHETA,RMOM0,1,R,1)                                   ABSH0303
      CR04 = SDOT(NTHETA,RMOM0,1,GS4,1)                                 ABSH0304
      CR08 = SDOT(NTHETA,RMOM0,1,GS8,1)                                 ABSH0305
      CR0SS2 = SDOT(NTHETA,RMOM0,1,SS20,1)                              ABSH0306
      CPSI10 = SDOT(NTHETA,TMOM,1,GS10,1)                               ABSH0307
      CPSI12 = SDOT(NTHETA,TMOM,1,GS12,1)                               ABSH0308
      CPSI21 = SDOT(NTHETA,TMOM,1,GS21,1)                               ABSH0309
      CPSI22 = SDOT(NTHETA,TMOM,1,GS22,1)                               ABSH0310
      CPSI3 = SDOT(NTHETA,TMOM,1,R,1)                                   ABSH0311
      CPSI4 = SDOT(NTHETA,TMOM,1,GS4,1)                                 ABSH0312
      CPSI8 = SDOT(NTHETA,TMOM,1,GS8,1)                                 ABSH0313
      CSISS2 = SDOT(NTHETA,TMOM,1,SS20,1)                               ABSH0314
      CCS1 = SDOT(NTHETA,WFCN,1,CS1N,1)/CPSI8                           ABSH0315
      CR0GS4 = CR04/CPSI4                                               ABSH0316
      CR0GS8 = CR08/CPSI8                                               ABSH0317
      DR010 = CR010 - CR0GS8*CPSI10                                     ABSH0318
      DR012 = CR012 - CR0GS8*CPSI12                                     ABSH0319
      DR021 = CR021 - CR0GS8*CPSI21                                     ABSH0320
      DR022 = CR022 - CR0GS8*CPSI22                                     ABSH0321
      DR03 = CPSI8**2/CNORM**2*(CR03-CR0GS4*CPSI3)                      ABSH0322
      DR040 = CR0SS2 - CR0GS4*CSISS2                                    ABSH0323
      DR041 = CR08 - CR0GS4*CPSI8                                       ABSH0324
      DR05 = CCS1*(CR08-CR0GS4*CPSI8)                                   ABSH0325
      AL11 = DR010                                                      ABSH0326
      AL12 = -(X*DR021+Y(5)*DR022)                                      ABSH0327
      AL13 = DR012 - Y(3)*DR022                                         ABSH0328
      BL1 = -AI*DR03 - DR040 - DR041*AIPOI + DR05 +                     ABSH0329
     *     2.0*Y(4)*(DR021+Y(6)*DR022)                                  ABSH0330
      IF (KMHD.GT.1) GO TO 20                                           ABSH0331
C*****SOLUTION OF MATRIX FOR KMHD=1.                                    ABSH0332
      B3 = 2.0*D0                                                       ABSH0333
      B2 = 2.0*(E0-E1)                                                  ABSH0334
      B1 = (BL1-AL13*BL3-AL12*BL2)/AL11                                 ABSH0335
      GO TO 40                                                          ABSH0336
C*****Z0 MOMENT TERMS - ELONGATION.                                     ABSH0337
   20 CZ110 = SDOT(NTHETA,ZMOM1,1,GS10,1)                               ABSH0338
      CZ112 = SDOT(NTHETA,ZMOM1,1,GS12,1)                               ABSH0339
      CZ121 = SDOT(NTHETA,ZMOM1,1,GS21,1)                               ABSH0340
      CZ122 = SDOT(NTHETA,ZMOM1,1,GS22,1)                               ABSH0341
      CZ13 = SDOT(NTHETA,ZMOM1,1,R,1)                                   ABSH0342
      CZ14 = SDOT(NTHETA,ZMOM1,1,GS4,1)                                 ABSH0343
      CZ18 = SDOT(NTHETA,ZMOM1,1,GS8,1)                                 ABSH0344
      CZ1SS2 = SDOT(NTHETA,ZMOM1,1,SS20,1)                              ABSH0345
      CZ1GS4 = CZ14/CPSI4                                               ABSH0346
      CZ1GS8 = CZ18/CPSI8                                               ABSH0347
      DZ110 = CZ110 - CZ1GS8*CPSI10                                     ABSH0348
      DZ112 = CZ112 - CZ1GS8*CPSI12                                     ABSH0349
      DZ121 = CZ121 - CZ1GS8*CPSI21                                     ABSH0350
      DZ122 = CZ122 - CZ1GS8*CPSI22                                     ABSH0351
      DZ13 = CPSI8**2/CNORM**2*(CZ13-CZ1GS4*CPSI3)                      ABSH0352
      DZ140 = CZ1SS2 - CZ1GS4*CSISS2                                    ABSH0353
      DZ141 = CZ18 - CZ1GS4*CPSI8                                       ABSH0354
      DZ15 = CCS1*(CZ18-CZ1GS4*CPSI8)                                   ABSH0355
      AL21 = DZ110                                                      ABSH0356
      AL22 = -(X*DZ121+Y(5)*DZ122)                                      ABSH0357
      AL23 = DZ112 - Y(3)*DZ122                                         ABSH0358
      BL2 = -AI*DZ13 - DZ140 - DZ141*AIPOI + DZ15 +                     ABSH0359
     *     2.0*Y(4)*(DZ121+Y(6)*DZ122)                                  ABSH0360
      IF (KMHD.EQ.3) GO TO 30                                           ABSH0361
C*****SOLUTION OF MATRIX FOR KMHD=2.                                    ABSH0362
      B3 = 2.0*D0                                                       ABSH0363
      BL1 = BL1 - AL13*BL3                                              ABSH0364
      BL2 = BL2 - AL23*BL3                                              ABSH0365
      DET = AL11*AL22 - AL12*AL21                                       ABSH0366
      IF (DET.EQ.0.0) GO TO 50                                          ABSH0367
      B1 = (BL1*AL22-AL12*BL2)/DET                                      ABSH0368
      B2 = (AL11*BL2-BL1*AL21)/DET                                      ABSH0369
      GO TO 40                                                          ABSH0370
C*****R2 MOMENT TERMS - TRIANGULARITY.                                  ABSH0371
   30 CR210 = SDOT(NTHETA,RMOM2,1,GS10,1)                               ABSH0372
      CR212 = SDOT(NTHETA,RMOM2,1,GS12,1)                               ABSH0373
      CR221 = SDOT(NTHETA,RMOM2,1,GS21,1)                               ABSH0374
      CR222 = SDOT(NTHETA,RMOM2,1,GS22,1)                               ABSH0375
      CR23 = SDOT(NTHETA,RMOM2,1,R,1)                                   ABSH0376
      CR24 = SDOT(NTHETA,RMOM2,1,GS4,1)                                 ABSH0377
      CR28 = SDOT(NTHETA,RMOM2,1,GS8,1)                                 ABSH0378
      CR2SS2 = SDOT(NTHETA,RMOM2,1,SS20,1)                              ABSH0379
      CR2GS4 = CR24/CPSI4                                               ABSH0380
      CR2GS8 = CR28/CPSI8                                               ABSH0381
      DR210 = CR210 - CR2GS8*CPSI10                                     ABSH0382
      DR212 = CR212 - CR2GS8*CPSI12                                     ABSH0383
      DR221 = CR221 - CR2GS8*CPSI21                                     ABSH0384
      DR222 = CR222 - CR2GS8*CPSI22                                     ABSH0385
      DR23 = CPSI8**2/CNORM**2*(CR23-CR2GS4*CPSI3)                      ABSH0386
      DR240 = CR2SS2 - CR2GS4*CSISS2                                    ABSH0387
      DR241 = CR28 - CR2GS4*CPSI8                                       ABSH0388
      DR25 = CCS1*(CR28-CR2GS4*CPSI8)                                   ABSH0389
      AL31 = DR210                                                      ABSH0390
      AL32 = -(X*DR221+Y(5)*DR222)                                      ABSH0391
      AL33 = DR212 - Y(3)*DR222                                         ABSH0392
      BL3 = -AI*DR23 - DR240 - DR241*AIPOI + DR25 +                     ABSH0393
     *     2.0*Y(4)*(DR221+Y(6)*DR222)                                  ABSH0394
C*****SOLUTION OF MATRIX FOR KMHD=3.                                    ABSH0395
      DET1 = AL22*AL33 - AL23*AL32                                      ABSH0396
      DET2 = AL21*AL33 - AL23*AL31                                      ABSH0397
      DET3 = AL21*AL32 - AL22*AL31                                      ABSH0398
      DET = AL11*DET1 - AL12*DET2 + AL13*DET3                           ABSH0399
      IF (DET.EQ.0.0) GO TO 50                                          ABSH0400
      B1 = BL1*DET1 - AL12*(BL2*AL33-AL23*BL3) +                        ABSH0401
     *     AL13*(BL2*AL32-AL22*BL3)                                     ABSH0402
      B1 = B1/DET                                                       ABSH0403
      B2 = AL11*(BL2*AL33-AL23*BL3) - BL1*DET2 +                        ABSH0404
     *     AL13*(AL21*BL3-BL2*AL31)                                     ABSH0405
      B2 = B2/DET                                                       ABSH0406
      B3 = AL11*(AL22*BL3-BL2*AL32) - AL12*(AL21*BL3-BL2*AL31)          ABSH0407
     *     + BL1*DET3                                                   ABSH0408
      B3 = B3/DET                                                       ABSH0409
   40 YP(2) = B1                                                        ABSH0410
      YP(4) = B2                                                        ABSH0411
      YP(6) = B3                                                        ABSH0412
      IUNDER = 0                                                        ABSH0413
      RETURN                                                            ABSH0414
C*****SINGULAR MATRIX--USUALLY UNDERFLOW IN DET.                        ABSH0415
C*****SOMETIMES OCCURS WHEN THE SOLUTION OF THE MOMENT EQUATIONS IS     ABSH0416
C*****ABOUT TO FAIL.                                                    ABSH0417
   50 CONTINUE                                                          ABSH0418
C--      WRITE(NTTY,1010) X,Y                                           ABSH0419
C-- 1010 FORMAT(38H SINGULAR SYSTEM IN "AUX" AT X= ,F5.3,/,4H Y= ,6E10.2ABSH0420
      IUNDER = 1                                                        ABSH0421
      RETURN                                                            ABSH0422
      END                                                               ABSH0423

max_stars_repo_name: jonathanschilling/VMOMS

max_issues_repo_name: jonathanschilling/VMOMS

max_forks_repo_name: jonathanschilling/VMOMS


!> This is a typical NGA driver.
!> @version 1.0
!> @author O. Desjardins
program nga
   use random,     only: random_initialize
   use param,      only: param_init,param_final
   use parallel,   only: parallel_init,parallel_final
   use messager,   only: messager_init,messager_final
   use geometry,   only: geometry_init
   use simulation, only: simulation_init,simulation_run,simulation_final
   implicit none

   ! =======================================
   ! Initialization sequence ===============
   ! =======================================
   ! Initialize parallel environment =======
   call parallel_init
   ! Initialize messaging capabilities ====
   call messager_init
   ! Initialize user interaction ===========
   call param_init
   ! Initialize random number generator ====
   call random_initialize
   ! =======================================

   ! =======================================
   ! Setup the problem grid ================
   ! =======================================
   call geometry_init
   ! =======================================

   ! =======================================
   ! Setup the solver ======================
   ! =======================================
   call simulation_init
   ! =======================================

   ! =======================================
   ! Run the solver ========================
   ! =======================================
   call simulation_run
   ! =======================================

   ! =======================================
   ! Finalize the solver ===================
   ! =======================================
   call simulation_final
   ! =======================================

   ! Code termination ======================
   ! Terminate solvers, data, and geometry
   !call simulation_final
   !call data_final
   !call geometry_final

   ! Terminate I/O, parser, timing, and monitor
   !call io_final

   ! =======================================
   ! Termination sequence ==================
   ! =======================================
   ! Clean up user parameters ==============
   call param_final
   ! Clean up monitoring tools =============
   call messager_final
   ! Clean parallel exit ===================
   call parallel_final
   ! =======================================

end program nga


text:
      subroutine sortpol(l,m,n)

c     Check that polynomes and monomes are printed in an
c     orderly fashion

      real t(100), u(100,100,100)

c     check for linear complexity, expected: l + m + n + ...
      if(n.gt.m) then
         do i = 1, l
            t(i) = 0.
         enddo
         do k = 1, n
            t(k) = 0.
         enddo
         do j = 1, m
            t(j) = 0.
         enddo
      endif

c     check for global correctness for all thrid degree terms, expected:
c     l^3 + l^2*m + l^2*n +
      if(n.gt.m) then
c     check for monome sort, expected: l*m*n
         do i = 1, l
            do k = 1, n
               do j = 1, m
                  u(i,j,k) = 0.
               enddo
            enddo
         enddo
c     three power of 3
         do i = 1, l
            do k = 1, l
               do j = 1, l
                  u(i,j,k) = 0.
               enddo
            enddo
         enddo
         do i = 1, n
            do k = 1, n
               do j = 1, n
                  u(i,j,k) = 0.
               enddo
            enddo
         enddo
         do i = 1, m
            do k = 1, m
               do j = 1, m
                  u(i,j,k) = 0.
               enddo
            enddo
         enddo

c     six third degree terms with one power of 2

         do i = 1, l
            do k = 1, l
               do j = 1, m
                  u(i,j,k) = 0.
               enddo
            enddo
         enddo
         do i = 1, l
            do k = 1, n
               do j = 1, l
                  u(i,j,k) = 0.
               enddo
            enddo
         enddo

         do i = 1, n
            do k = 1, n
               do j = 1, m
                  u(i,j,k) = 0.
               enddo
            enddo
         enddo
         do i = 1, l
            do k = 1, n
               do j = 1, n
                  u(i,j,k) = 0.
               enddo
            enddo
         enddo

         do i = 1, m
            do k = 1, n
               do j = 1, m
                  u(i,j,k) = 0.
               enddo
            enddo
         enddo
         do i = 1, l
            do k = 1, m
               do j = 1, m
                  u(i,j,k) = 0.
               enddo
            enddo
         enddo

      endif

      end


text:
subroutine solve_solid_disp_pa(x,x_curr,kid,ien,node_sbc,xpre1,solid_prevel,solid_preacc,ien_sbc,&
			solid_stress,solid_bcvel,mtype,&
			ne_local,ien_local,nn_local,node_local,send_address,ad_length,&
			global_com,nn_global_com,local_com,nn_local_com)
! subroutine to solve the solid displacement at every time step
! x  -- initial solid configuration
! x_curr --- current solid configuration
! disp --- solid displacement
! kid --- 1st type boundary condition nodes id
! ien --- solid connective matrix
! xpre1 --- last 1 time step solid position
! xpre2 --- last 2 time step solid position
! solid_acc --- solid acceleration at n+1
! solid_preacc --- solid acceleration at n
! solid_vel --- solid velocity at n+1
! solid_prevel --- solid velocity at n
! solid_stress --- solid traction B.C.
! solid_bcvel --- solid velocity B.C.

use mpi_variables
use solid_variables
use run_variables, only: dt
implicit none
include 'mpif.h'

real(8) x(nsd_solid,nn_solid)
real(8) x_curr(nsd_solid,nn_solid)
real(8) disp(nsd_solid,nn_solid)
integer kid(nsd_solid,nn_solid)
integer ien(ne_solid,nen_solid)
integer node_sbc(nn_sbc)
real(8) xpre1(nsd_solid,nn_solid)
real(8) solid_prevel(nsd_solid,nn_solid)
real(8) solid_preacc(nsd_solid,nn_solid)
real(8) solid_acc(nsd_solid,nn_solid)
real(8) solid_vel(nsd_solid,nn_solid)
real(8) solid_bcvel(nsd_solid,nn_solid)
integer mtype(ne_solid)
!real(8) xpre2(nsd_solid,nn_solid)
!real(8) disp1(nsd_solid,nn_solid)
!real(8) disp2(nsd_solid,nn_solid)
integer ien_sbc(ne_sbc,nen_solid+2)
real(8) solid_stress(nsd_solid*2,nn_solid)
!-------------------------------------------------
real(8) sq_solid(0:3,8,8)
integer iq
!real(8) :: xq_solid(nsdpad_solid,nquadpad_solid),wq(nquadpad_solid)
!-------------------------------------------------
integer inner
integer outer
parameter (inner = 300) ! solid equation inner 50 should be sufficient
parameter (outer = 10)  ! solid equation outer 5 should be sufficient
!-------------------------------------------------
real(8) dg(nsd_solid,nn_solid) ! disp correction
real(8) w(nsd_solid,nn_solid)  ! pre-conditioner
real(8) p(nsd_solid,nn_solid)  ! residual vector
real(8) res
real(8) del
integer i
integer node
!----------------------------
        real(8) alpha
        real(8) beta
        real(8) gama
!----------------------------
! MPI mesh partition variables
integer ne_local
integer ien_local(ne_local)
integer nn_local
integer node_local(nn_local)
integer ad_length
integer send_address(ad_length,2)
integer nn_global_com
integer global_com(nn_global_com)  ! global node index for communication
integer nn_local_com
integer local_com(nn_local_com)  ! local index in the communication region on each processor

real(8) solid_bcforce(nsd_solid,nn_solid)
! define the numerical parameters
alpha = -0.05 ! -1/3 < alpha < 0 and alpha == 0 is Newmark method
gama = (1.0 - 2.0 * alpha) * 0.5
beta = ( (1.0 - alpha)**2 ) * 0.25
!beta = (1.0 - alpha**2 ) * 0.25
!----------------------------------

! Get sq for solid mesh to calculate sh in block_solid.f90
if (nsd_solid == 3) then
     do iq=1,nquad_solid
                if(nen_solid.eq.4) then
                  sq_solid(0,1,iq) = xq_solid(1,iq)
                  sq_solid(0,2,iq) = xq_solid(2,iq)
                  sq_solid(0,3,iq) = xq_solid(3,iq)
                  sq_solid(0,4,iq) = 1 - xq_solid(1,iq) - xq_solid(2,iq) - xq_solid(3,iq)
        else
                  sq_solid(0,1,iq) = (1 - xq_solid(1,iq))   &
                           * (1 - xq_solid(2,iq)) * (1 - xq_solid(3,iq)) / 8
                  sq_solid(0,2,iq) = (1 + xq_solid(1,iq))   &
                           * (1 - xq_solid(2,iq)) * (1 - xq_solid(3,iq)) / 8
                  sq_solid(0,3,iq) = (1 + xq_solid(1,iq))   &
                           * (1 + xq_solid(2,iq)) * (1 - xq_solid(3,iq)) / 8
                  sq_solid(0,4,iq) = (1 - xq_solid(1,iq))   &
                           * (1 + xq_solid(2,iq)) * (1 - xq_solid(3,iq)) / 8
                  sq_solid(0,5,iq) = (1 - xq_solid(1,iq))   &
                           * (1 - xq_solid(2,iq)) * (1 + xq_solid(3,iq)) / 8
                  sq_solid(0,6,iq) = (1 + xq_solid(1,iq))   &
                           * (1 - xq_solid(2,iq)) * (1 + xq_solid(3,iq)) / 8
                  sq_solid(0,7,iq) = (1 + xq_solid(1,iq))   &
                           * (1 + xq_solid(2,iq)) * (1 + xq_solid(3,iq)) / 8
                  sq_solid(0,8,iq) = (1 - xq_solid(1,iq))   &
                           * (1 + xq_solid(2,iq)) * (1 + xq_solid(3,iq)) / 8
                  sq_solid(1,1,iq) = - (1 - xq_solid(2,iq)) * (1 - xq_solid(3,iq)) / 8
                  sq_solid(1,2,iq) = + (1 - xq_solid(2,iq)) * (1 - xq_solid(3,iq)) / 8
                  sq_solid(1,3,iq) = + (1 + xq_solid(2,iq)) * (1 - xq_solid(3,iq)) / 8
                  sq_solid(1,4,iq) = - (1 + xq_solid(2,iq)) * (1 - xq_solid(3,iq)) / 8
                  sq_solid(1,5,iq) = - (1 - xq_solid(2,iq)) * (1 + xq_solid(3,iq)) / 8
                  sq_solid(1,6,iq) = + (1 - xq_solid(2,iq)) * (1 + xq_solid(3,iq)) / 8
                  sq_solid(1,7,iq) = + (1 + xq_solid(2,iq)) * (1 + xq_solid(3,iq)) / 8
                  sq_solid(1,8,iq) = - (1 + xq_solid(2,iq)) * (1 + xq_solid(3,iq)) / 8
                  sq_solid(2,1,iq) = - (1 - xq_solid(1,iq)) * (1 - xq_solid(3,iq)) / 8
                  sq_solid(2,2,iq) = - (1 + xq_solid(1,iq)) * (1 - xq_solid(3,iq)) / 8
                  sq_solid(2,3,iq) = + (1 + xq_solid(1,iq)) * (1 - xq_solid(3,iq)) / 8
                  sq_solid(2,4,iq) = + (1 - xq_solid(1,iq)) * (1 - xq_solid(3,iq)) / 8
                  sq_solid(2,5,iq) = - (1 - xq_solid(1,iq)) * (1 + xq_solid(3,iq)) / 8
                  sq_solid(2,6,iq) = - (1 + xq_solid(1,iq)) * (1 + xq_solid(3,iq)) / 8
                  sq_solid(2,7,iq) = + (1 + xq_solid(1,iq)) * (1 + xq_solid(3,iq)) / 8
                  sq_solid(2,8,iq) = + (1 - xq_solid(1,iq)) * (1 + xq_solid(3,iq)) / 8
                  sq_solid(3,1,iq) = - (1 - xq_solid(1,iq)) * (1 - xq_solid(2,iq)) / 8
                  sq_solid(3,2,iq) = - (1 + xq_solid(1,iq)) * (1 - xq_solid(2,iq)) / 8
                  sq_solid(3,3,iq) = - (1 + xq_solid(1,iq)) * (1 + xq_solid(2,iq)) / 8
                  sq_solid(3,4,iq) = - (1 - xq_solid(1,iq)) * (1 + xq_solid(2,iq)) / 8
                  sq_solid(3,5,iq) = + (1 - xq_solid(1,iq)) * (1 - xq_solid(2,iq)) / 8
                  sq_solid(3,6,iq) = + (1 + xq_solid(1,iq)) * (1 - xq_solid(2,iq)) / 8
                  sq_solid(3,7,iq) = + (1 + xq_solid(1,iq)) * (1 + xq_solid(2,iq)) / 8
                  sq_solid(3,8,iq) = + (1 - xq_solid(1,iq)) * (1 + xq_solid(2,iq)) / 8
        endif
          enddo
elseif (nsd_solid == 2) then

      do iq=1,nquad_solid
                if(nen_solid==3) then
                  sq_solid(0,1,iq) = xq_solid(1,iq)
                  sq_solid(0,2,iq) = xq_solid(2,iq)
                  sq_solid(0,3,iq) = 1 - xq_solid(1,iq) - xq_solid(2,iq)
        elseif (nen_solid==4) then
                  sq_solid(0,1,iq) = (1 - xq_solid(1,iq)) * (1 - xq_solid(2,iq)) / 4
                  sq_solid(0,2,iq) = (1 + xq_solid(1,iq)) * (1 - xq_solid(2,iq)) / 4
                  sq_solid(0,3,iq) = (1 + xq_solid(1,iq)) * (1 + xq_solid(2,iq)) / 4
                  sq_solid(0,4,iq) = (1 - xq_solid(1,iq)) * (1 + xq_solid(2,iq)) / 4

                  sq_solid(1,1,iq) = - (1 - xq_solid(2,iq)) / 4
                  sq_solid(1,2,iq) = + (1 - xq_solid(2,iq)) / 4
                  sq_solid(1,3,iq) = + (1 + xq_solid(2,iq)) / 4
                  sq_solid(1,4,iq) = - (1 + xq_solid(2,iq)) / 4

                  sq_solid(2,1,iq) = - (1 - xq_solid(1,iq)) / 4
                  sq_solid(2,2,iq) = - (1 + xq_solid(1,iq)) / 4
                  sq_solid(2,3,iq) = + (1 + xq_solid(1,iq)) / 4
                  sq_solid(2,4,iq) = + (1 - xq_solid(1,iq)) / 4

        endif
          enddo
end if
!------------------------------------------------------------------
solid_acc(:,:) = 0.0d0
solid_vel(:,:) = 0.0d0
solid_acc(:,:) = solid_preacc(:,:)

!do i=1,nn_sbc
!        solid_acc(:,node_sbc(i))=(solid_bcvel(:,node_sbc(i)) - solid_prevel(:,node_sbc(i))) / dt
!end do

w(:,:)=0.0d0
p(:,:)=0.0d0
dg(:,:)=0.0d0
!p_rec(:,:) = 0.0

! Evaluate residual for solid equations in parallel
call block_solid_pa(x,solid_acc,w,p,ien,nsd_solid,nen_solid,ne_solid,&
		nn_solid,nquad_solid,wq_solid,sq_solid,xpre1,&
		solid_prevel,solid_preacc,ien_sbc,ne_sbc,solid_stress,mtype,&
		ne_local,ien_local)

! commute residual for each processor

call communicate_res_ad_subsolid(p,nsd_solid,nn_solid,send_address,ad_length)
call communicate_res_ad_subsolid(w,nsd_solid,nn_solid,send_address,ad_length)

!======================================
! Apply 2nd type boundary
! The MPI version the 2nd type BC has to be applied here
    if (nsd_solid ==  2) then
	call apply_2ndbc_solid2d(xpre1,nsd_solid,nn_solid,ien_sbc,ne_sbc,nen_solid,ien,ne_solid,solid_bcforce,solid_stress)
	do i=1,nn_local
	node=node_local(i)
	p(1:nsd_solid,node) = p(1:nsd_solid,node) + solid_bcforce(1:nsd_solid,node)
	end do
	else
	call apply_2ndbc_solid(xpre1,nsd_solid,nn_solid,ien_sbc,ne_sbc,nen_solid,ien,ne_solid,solid_bcforce,solid_stress)
	p(:,:) = p(:,:) + solid_bcforce(:,:)
	end if

! Set 1st B.C.
call setid_pa(p,nsd_solid,nn_solid,kid,node_local,nn_local)
call getnorm_pa(p,nsd_solid,nn_solid,node_local,nn_local,res)
res=sqrt(res)

if (myid == 0) write(*,*) '===Initial error for solid displacement===', res

!-----------------------
! Take the invese of w as the preconditioner
! Help a lot!!!! However, have not figured out why ...
!w(:,:)=1.0d0/w(:,:)
!----------------------

call gmres_solid_pa(x,w,p,dg,ien,kid,nsd_solid,nn_solid,ne_solid,nen_solid,inner,outer,&
		nquad_solid,wq_solid,sq_solid,xpre1,&
		solid_prevel,solid_preacc,solid_stress,ne_sbc,ien_sbc,mtype,&
		ne_local,ien_local,node_local,nn_local,&
		global_com,nn_global_com,local_com,nn_local_com,send_address,ad_length,nsd_solid)

call getnorm_pa(dg,nsd_solid,nn_solid,node_local,nn_local,del)
 del = sqrt(del)
if (myid == 0) write(*,*) '===solid displacement correction norm===', del

solid_acc(:,:) = solid_acc(:,:) + dg(:,:)
x_curr(:,:) = xpre1(:,:) + dt*solid_prevel(:,:) +&
		 (dt**2)*0.5*( (1.0-2.0*beta)*solid_preacc(:,:) + 2.0*beta*solid_acc(:,:) )

solid_vel(:,:) = solid_prevel(:,:) + dt*( (1-gama)*solid_preacc(:,:) + gama*solid_acc(:,:) )

!================================
! Just for testing - pure explicit 2nd order
!x_curr(:,:) = xpre1(:,:) + dt*solid_vel(:,:)
!================================
!if (myid == 0) then
!write(*,*) "===============", x_curr(1,1), "%%%%%%%%%%%%%%%%%%%%%%%"
!write(*,*) "===============", solid_vel(1,1)*dt, "%%%%%%%%%%%%%%%%%%%%%%%"
!end if

solid_prevel(:,:) = solid_vel(:,:)
solid_preacc(:,:) = solid_acc(:,:)

return
end subroutine solve_solid_disp_pa

max_stars_repo_name: biofluids/IFEM-archive

max_issues_repo_name: biofluids/IFEM-archive

max_forks_repo_name: biofluids/IFEM-archive


text:
program problem16
  implicit none
  integer*16, dimension(16) :: n = (/0_16,0_16,0_16,0_16, 0_16,0_16,0_16,0_16, 0_16,0_16,0_16,0_16, 0_16,0_16,0_16,shiftl(1_16,40)/)
  integer                   :: i,s

  s=0

  ! Repeatedly divide the 1024-bit number by 10 to get the digits of
  ! the decimal representation.  The 1024-bit number is represented as
  ! 16 128-bit numbers (instead of 16 64-bit numbers) because extra
  ! precision is needed to handle the carries during the long
  ! division.
  do while (sum(n) > 0)
     do i=16,1,-1
        if (i > 1) then
           n(i-1)=shiftl(mod(n(i),10),64)+n(i-1) ! carry
        else
           s=s+mod(n(i),10)                      ! tally
        end if
        n(i)=n(i)/10
     end do
  end do

  print *, s

end program problem16

max_stars_repo_name: jamesmcclain/Dioskouroi

max_issues_repo_name: jamesmcclain/Dioskouroi

max_forks_repo_name: jamesmcclain/Dioskouroi


text:
      parameter(Nx=360,Ny=180)
      parameter(Mrey=529,MTP=364)
      parameter(MTPstrt=12)

      real TPtme93(MTP),Reytme93(Mrey),x(2),y(2)
      real sst1(Nx,Ny),sst2(Nx,Ny),sst(Nx,Ny)

      open(20,file='REYv2_sst_93-03_wklypanom',form='unformatted',
     *access='direct',status='old',recl=Nx*Ny*4)
      open(30,file='REY_sstpanom_11-364',form='unformatted',
     *access='direct',status='unknown',recl=Nx*Ny*4)

! set up TP times   (XI)  days
      TPtme93(11)=5.25
      do irep=12,MTP
	TPtme93(irep)=TPtme93(irep-1)+9.9156
      enddo

! set up REYNOLDS times (X)  days
      Reytme93(1)=6.5
      do irec=2,Mrey
	Reytme93(irec)=Reytme93(irec-1)+7
      enddo

      IREC1=2
      DO NR=MTPstrt,MTP

        XI=TPtme93(NR)
	do irec=IREC1,Mrey
	write(6,*)NR,XI,Reytme93(irec-1),Reytme93(irec)
	if(XI.ge.Reytme93(irec-1).and.XI.lt.Reytme93(irec))then
	  x(1)=Reytme93(irec-1)
	  x(2)=Reytme93(irec)
	  IREC1=irec
	  go to 100
	endif
	enddo
	STOP1

  100  read(20,rec=irec-1)sst1
       read(20,rec=irec)sst2
       do j=1,Ny
       do i=1,Nx
	 if(sst1(i,j).eq.-9999..or.sst2(i,j).eq.-9999.)then
	   sst(i,j)=-9999.
	 else
	   y(1)=sst1(i,j)
	   y(2)=sst2(i,j)
	   call xinterp(x,y,XI,2,sst(i,j))
	 endif
       enddo
       enddo
       write(30,rec=NR)sst
	enddo

	write(6,*)' JOB COMPLETE'
	stop
	end
      subroutine xINTERP (X,Y,XI,N,yinterp)
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c  MUST BE MONOTONIC INCREASING!!!!!!!!!!!!!!!!!!!
c  note :  real*4 version
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      REAL X(n),Y(n),xi,yinterp
C
      if(x(2).lt.x(1))then
	 write(6,*)' STOPPiNG XINTERP....not monotonic increasing'
	 stop
      endif
      IF (XI.LT.X(1)) GO TO 800
C
      DO 100 I=2,N
      IF (XI.LE.X(I)) GO TO 400
 100  CONTINUE
      yinterp=Y(N)+((Y(N)-Y(N-1))*(XI-X(N)))/(X(N)-X(N-1))
      write(6,*)' outside range ..XI > x(n)'
      RETURN
C
  400 yinterp=(Y(I-1)*(X(I)-XI)-Y(I)*(X(I-1)-XI))/(X(I)-X(I-1))
      RETURN
C
 800  yinterp=Y(1)+((Y(2)-Y(1))*(XI-X(1)))/(X(2)-X(1))
      write(6,*)' outside range ..XI < x(1)'
      RETURN
      END


program main
  implicit none

  integer :: argc
  argc = command_argument_count ()

  !$acc parallel copyin(argc)
  if (argc .ne. 0) then
     call abort
  end if
  !$acc end parallel

end program main


text: #redirect November 18, 2011 UC Davis Police Response to Occupy UC Davis

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE zndata6 (IFLTAB, IFPOS, JULS, ISECS, CPATH, NPATH,
     *                   JULD, ISECD, IDTYPE, ISTAT)
C
C
C     New Data
C     Return pathnames for new data written to the file since
C     a given date/time.  Pathnames are returned one at a time,
C     for subsequent calls until the end of the file is reached.
C     To start list, set IFPOS to 0.  At the end of the list
C     ISTAT will be set to 1 (with no pathname returned), otherwise
C     ISTAT will be 0 for a normal return with pathname, or <0 for
C     an (somewhat serious) error return.
C
C      Input:
C        IFLTAB - From zopen6
C        IFPOS - File position, set to zero at start
C        ISTIME - Start time, in minutes since Jan 1, 1900,
C                 to check for data later than
C
C     Output:
C        CPATH - Pathname returned
C        NPATH - Length of pathname
C        IPTIME - Time of last write
C        ISTAT -  <0:  Error
C                  0:  Pathname retrieved okay
C                  1:  Reached end of file.  No pathname returned
C
C     Written by Bill Charley at HEC, 2004.
C
C
C
      INCLUDE 'zdsskz.h'
C
      INCLUDE 'zdsslz.h'
C
      INCLUDE 'zdsscz.h'
C
      INCLUDE 'zdssbz.h'
C
      INCLUDE 'zdssnz.h'
C
      INCLUDE 'zdssiz.h'
C
      INCLUDE 'zdssmz.h'
C
C
      INTEGER IFLTAB(*), IFPOS, NPATH, ISTAT
      CHARACTER CPATH*(*)

C
      INTEGER IPBIN(128)
      LOGICAL LMATCH
      CHARACTER CCDATE*10, CCTIME*10
C
C
C
      IF (MLEVEL.GE.11) WRITE (MUNIT,20) IFLTAB(KUNIT)
 20   FORMAT (T6,'-----DSS---Debug:  Enter zndata6;  Unit:',I5)
C
C
C     Have we reached the end of the file
      IF (IFPOS.GE.IFLTAB(KFSIZE)) GO TO 850
C
C
C
      IF (IFPOS.LE.0) THEN
C        Get the first pathname bin address, and
C        the address of the first bin
         CALL zrdprm6(IFLTAB, .FALSE.)
         NBIN = IFLTAB(KBNBLK)
         IADD = IFLTAB(KAFBIN)
         JPBIN = 1
      ELSE
         IADD   = IFPOS
         JPBIN  = IFLTAB(KEY3+1)
         NBIN   = IFLTAB(KEY3+2)
         ISELCA = IFLTAB(KEY3+3)
C
         IF (JPBIN.GT.(IFLTAB(KBNSIZ)-2)) GO TO 200
      ENDIF
C
C
C     Now read pathname bin from the file
 40   CONTINUE
C     Do a quick check of the address
      IF (IADD.LE.0) THEN
         ISTAT = -2
         GO TO 900
      ENDIF
      IF (IADD.GT.IFLTAB(KFSIZE)) THEN
C         Get an update of the file size
          CALL zrdprm6(IFLTAB, .FALSE.)
C         If the address is just a little off, then someone else
C         is writing to the file while we are reading it.
C         If it is way off, then it is a file error
          ISTAT = -3
          IF ((IADD+20000).GT.IFLTAB(KFSIZE)) GO TO 900
		IF (IADD.GT.IFLTAB(KFSIZE)) GO TO 850
      ENDIF
      CALL zgtrec6 (IFLTAB, IPBIN, IFLTAB(KBNSIZ), IADD, .FALSE.)
C
C
C     Loop through bin, looking for pathnames
 100  CONTINUE
C     Any more pathnames left?
      IF (IPBIN(JPBIN).LE.0) GO TO 200
C     Yes - Compute the number of integer words in the pathname
      NBPWPA = ((IPBIN(JPBIN+KBNPAT)-1) / NCPW) + 1
      NBMWPA = ((IPBIN(JPBIN+KBNPAT)-1) / NCMW) + 1
C
C     Record status good? (not deleted or renamed?)
      IF ((IPBIN(JPBIN).EQ.1).OR.(IPBIN(JPBIN).EQ.11)) THEN
C        Valid record.  Get pathname and pathname length
         NPATH = IPBIN(JPBIN+KBNPAT)
         CALL HOL2CH (IPBIN(JPBIN+KBPATH), CKPATH, NBMWPA)
C        Read its information block to compare last written times
         NADD = IPBIN(JPBIN+NBPWPA+KBAINF)
         CALL zgtrec6 (IFLTAB, INFO, NINFO+NBPWPA, NADD, .FALSE.)
C
C        Extract information from this array
         IF (INFO(KIFLAG).NE.NPFLAG) GO TO 900
         IDTYPE = INFO(NBPWPA+KITYPE)
         CCDATE = ' '
         CCTIME = ' '
         CALL HOLCHR (INFO(NBPWPA+KIDATE), 1, NDATEC, CCDATE, 1)
         CALL HOLCHR (INFO(NBPWPA+KITIME), 1, NTIMEC, CCTIME, 1)
         CALL DATJUL(CCDATE(1:NDATEC), JULD, IERR)
         ISECD = IHMS2S(CCTIME(1:NTIMEC))
	   LMATCH = .FALSE.
         IF ((JULD.GE.0).AND.(ISECD.GE.0)) THEN
            IF (JULD.GT.JULS) THEN
               LMATCH = .TRUE.
            ELSE IF (JULD.EQ.JULS) THEN
               IF (ISECD.GE.ISECS) LMATCH = .TRUE.
            ENDIF
         ENDIF
         IF (LMATCH) THEN
            CPATH = CKPATH
C	      Update the bin pointer
            JPBIN = JPBIN + NBPWPA + NLBIN
            GO TO 800
         ENDIF
      ENDIF
C
C     Do an error check on this block
      IF ((IPBIN(JPBIN).LT.-1000).OR.(IPBIN(JPBIN).GT.1000)) THEN
         ISTAT = -4
         GO TO 900
      ENDIF
C
C     Unmatched, deleted or renamed record.
C     Update the bin pointer (to next possible pathname
C     location  within this block).
      JPBIN = JPBIN + NBPWPA + NLBIN
C     Is that pointer too large?
      IF (JPBIN.GT.(IFLTAB(KBNSIZ)-2)) GO TO 200
C     Go back and look for next pathnme within this bin
      GO TO 100
C
C
 200  CONTINUE
C     At this point, there are no more pathnames in the current bin
C     (or it is full and extends into another block).
C     Read the next pathname bin
      NBIN = NBIN - 1
C     Any more bins in this block?
      IF (NBIN.LE.0) THEN
C        No - Get pointer to next bin block (section).
C        Get location of pointers in bin
         I = IFLTAB(KBNSIZ)
C        Any more bins in the DSS file? (Exit to 850 if no more).
         IF (IPBIN(I).EQ.0) GO TO 850
         IADD = IPBIN(I)
         NBIN = IFLTAB(KBNBLK)
      ELSE
C        More bins available within this block - get next one.
         IADD = IADD + IFLTAB(KBNSIZ)
      ENDIF
C
C     Go back up and read next bin
      JPBIN = 1
      GO TO 40
C
C
 800  CONTINUE
      ISTAT = 0
C     Save pointer information
      IFPOS = IADD
      IFLTAB(KEY3+1) = JPBIN
      IFLTAB(KEY3+2) = NBIN
      IFLTAB(KEY3+3) = ISELCA
C
 820  CONTINUE
      IF (MLEVEL.GE.11) WRITE (MUNIT,840)
 840  FORMAT (T6,'-----DSS---Debug:  Exit  zckpat6')
      RETURN
C
C
 850  CONTINUE
C     End of file
      ISTAT = 1
      IFPOS = IFLTAB(KFSIZE)
      CPATH = ' '
      NPATH = 0
      IFLTAB(KEY3+1) = 0
      IFLTAB(KEY3+2) = 0
      IFLTAB(KEY3+3) = 0
      GO TO 820
C
C
 900  CONTINUE
      IF (MLEVEL.GE.1) WRITE(MUNIT,920)
 920  FORMAT(' -----DSS---Error:  File error detected in zndata6')
C      CALL zerror6 (IFLTAB, 100, 'zndata6', 0, IFLTAB(KTABLE), ' ', 0,
C     * ' ', 0)
      IF (ISTAT.EQ.-2) THEN
      IF (MLEVEL.GE.1) WRITE(MUNIT,930) IADD
 930  FORMAT(' zndata6 Error: File address out of range:',I10)
      ELSE IF (ISTAT.EQ.-3) THEN
	IF (MLEVEL.GE.1) WRITE(MUNIT,930) IADD
      ELSE IF (ISTAT.EQ.-4) THEN
      IF (MLEVEL.GE.1) WRITE(MUNIT,940) IPBIN(JPBIN)
 940  FORMAT(' zndata6 Error: Pathname status out of range:',I10)
      ENDIF
      IF (ISTAT.GE.0) ISTAT = -1
      GO TO 820
C
      END

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


text:
     module netcdfai_var
       use netcdf
       use netcdfai_common
       use netcdfai_att
       implicit none

       interface nf90ai_get_var
         module procedure nf90ai_get_var_0D,nf90ai_get_var_1D, &
                          nf90ai_get_var_2D,nf90ai_get_var_3D, &
                          nf90ai_get_var_4D
       end interface nf90ai_get_var

       interface nf90ai_get_var_att
         module procedure nf90ai_get_var_att_text,nf90ai_get_var_att_1D
       end interface nf90ai_get_var_att

       contains

       subroutine nf90ai_new_var(ncid,varid,var)
         integer, intent(in) :: ncid,varid
         type(nf90aiVar), intent(inout) :: var
         integer :: attIdx

         allocate(var%dimids(NF90_MAX_VAR_DIMS))
         call nf90ai_check( &
           nf90_inquire_variable( &
             ncid = ncid, &
             varid = varid, &
             name = var%name, &
             xtype = var%xtype, &
             ndims = var%ndims,  &
             dimids =var%dimids, &
             natts= var%natts &
           ) &
         )

          allocate(var%atts(var%natts))
          do attIdx=1,var%natts
            call nf90ai_new_att(ncid,varid,attIdx,var%atts(attIdx))
          end do

       end subroutine nf90ai_new_var

       function nf90ai_get_var_0D(self,varName,value) result(returnVal)
         ! Dummy Arguments
         type(nf90aiMeta),intent(inout) :: self
         character(len=*),intent(in) :: varName
         real,intent(inout):: value
         integer :: returnVal

         ! Local Variables
         integer :: varIdx
         integer :: i

         if(nf90ai_get_var_index(self,varName,varIdx) == nf90_noerr) then

           if(self%vars(varIdx)%nDims==0) then
             !returnVal =  nf90_open(self%filename, NF90_NOWRITE, self%ncid)
             call nf90ai_check( &
               nf90_open(self%filename, NF90_NOWRITE, self%ncid) &
             )
             !returnVal = nf90_get_var(self%ncid,varIdx,value)
             !print *, 'get_var returnVal: ',returnVal
             call nf90ai_check( &
               nf90_get_var(self%ncid,varIdx,value) &
             )
             !call nf90ai_check( nf90_close(self%ncid) )
             returnVal = nf90_close(self%ncid)
             !print *, 'close returnVal: ',returnVal
             !print *, trim(self%vars(varIdx)%name), &
             !         self%vars(varIdx)%nDims, &
             !         varIdx,values
           else
             returnVal=-999
           endif
         else
           returnVal=-999
         endif

       end function nf90ai_get_var_0D

       function nf90ai_get_var_1D(self,varName,values) result(returnVal)
         ! Dummy Arguments
         type(nf90aiMeta),intent(inout) :: self
         character(len=*),intent(in) :: varName
         real,intent(inout),pointer :: values(:)
         integer :: returnVal

         ! Local Variables
         integer :: varIdx,dimSize
         integer :: i

         nullify(values)

         if(nf90ai_get_var_index(self,varName,varIdx) == nf90_noerr) then

           if(self%vars(varIdx)%nDims==1) then
             !returnVal = nf90_open(self%filename, NF90_NOWRITE, self%ncid)
             !print *, 'get_var_1D: open:',returnVal
             call nf90ai_check( &
               nf90_open(self%filename, NF90_NOWRITE, self%ncid) &
             )
             dimSize=self%dimSizes(self%vars(varIdx)%dimids(1))
             allocate(values(dimSize))
             call nf90ai_check( &
               nf90_get_var(self%ncid,varIdx,values) &
             )
             !returnVal =  nf90_get_var(self%ncid,varIdx,values)
             !print *, 'get_var_1D: get_var:',returnVal
             !call nf90ai_check( nf90_close(self%ncid) )
             returnVal = nf90_close(self%ncid)
             !print *, 'get_var_1D: get_var:',returnVal
             !print *, trim(self%vars(varIdx)%name), &
             !         self%vars(varIdx)%nDims, &
             !         varIdx,dimSize
             !print *, (values(i),i=1,dimSize)
           else
             returnVal=-999
           endif
         else
           returnVal=-999
         endif

       end function nf90ai_get_var_1D

       function nf90ai_get_var_2D(self,varName,values) result(returnVal)
         ! Dummy Arguments
         type(nf90aiMeta),intent(inout) :: self
         character(len=*),intent(in) :: varName
         real,intent(inout),pointer :: values(:,:)
         integer :: returnVal

         ! Local Variables
         integer :: varIdx,dimSize1,dimSize2
         integer :: i,j

         nullify(values)

         if(nf90ai_get_var_index(self,varName,varIdx) == nf90_noerr) then

           if(self%vars(varIdx)%nDims==2) then
             !returnVal = nf90_open(self%filename, NF90_NOWRITE, self%ncid)
             !print *, 'get_var_2D: open: ',returnVal
             call nf90ai_check( &
               nf90_open(self%filename, NF90_NOWRITE, self%ncid) &
             )
             dimSize1=self%dimSizes(self%vars(varIdx)%dimids(1))
             dimSize2=self%dimSizes(self%vars(varIdx)%dimids(2))
             allocate(values(dimSize1,dimSize2))
             !returnVal = nf90_get_var(self%ncid,varIdx,values)
             !print *, 'get_var_2D: get_var: ',returnVal
             call nf90ai_check( &
               nf90_get_var(self%ncid,varIdx,values) &
             )
             !call nf90ai_check( nf90_close(self%ncid) )
             returnVal = nf90_close(self%ncid)
             !print *, 'get_var_2D: close: ',returnVal
             !print *, trim(self%vars(varIdx)%name), &
             !         self%vars(varIdx)%nDims, &
             !         varIdx,dimSize1,dimSize2
             !print *, ((values(i,j),j=1,dimSize2),i=1,dimSize1)
           else
             returnVal=-999
           endif
         else
           returnVal=-999
         endif

       end function nf90ai_get_var_2D

       function nf90ai_get_var_3D(self,varName,values) result(returnVal)
         ! Dummy Arguments
         type(nf90aiMeta),intent(inout) :: self
         character(len=*),intent(in) :: varName
         real,intent(inout),pointer :: values(:,:,:)
         integer :: returnVal

         ! Local Variables
         integer :: varIdx,dimSize1,dimSize2,dimSize3
         integer :: i,j,k

         if(nf90ai_get_var_index(self,varName,varIdx) == nf90_noerr) then

           if(self%vars(varIdx)%nDims==3) then
             call nf90ai_check( &
               nf90_open(self%filename, NF90_NOWRITE, self%ncid) &
             )
             dimSize1=self%dimSizes(self%vars(varIdx)%dimids(1))
             dimSize2=self%dimSizes(self%vars(varIdx)%dimids(2))
             dimSize3=self%dimSizes(self%vars(varIdx)%dimids(3))
             allocate(values(dimSize1,dimSize2,dimSize3))
             call nf90ai_check( &
               nf90_get_var(self%ncid,varIdx,values) &
             )
             call nf90ai_check( nf90_close(self%ncid) )
             !print *, trim(self%vars(varIdx)%name), &
             !         self%vars(varIdx)%nDims, &
             !         varIdx,dimSize1,dimSize2,dimSize3
             !print *, (((values(i,j,k),k=1,dimSize3),j=1,dimSize2),i=1,dimSize1)
           else
             returnVal=-999
           endif
         else
           returnVal=-999
         endif

       end function nf90ai_get_var_3D

       function nf90ai_get_var_4D(self,varName,values) result(returnVal)
         ! Dummy Arguments
         type(nf90aiMeta),intent(inout) :: self
         character(len=*),intent(in) :: varName
         real,intent(inout),pointer :: values(:,:,:,:)
         integer :: returnVal

         ! Local Variables
         integer :: varIdx,dimSize1,dimSize2,dimSize3,dimSize4
         integer :: i,j,k,l

         if(nf90ai_get_var_index(self,varName,varIdx) == nf90_noerr) then

           if(self%vars(varIdx)%nDims==4) then
             call nf90ai_check( &
               nf90_open(self%filename, NF90_NOWRITE, self%ncid) &
             )
             dimSize1=self%dimSizes(self%vars(varIdx)%dimids(1))
             dimSize2=self%dimSizes(self%vars(varIdx)%dimids(2))
             dimSize3=self%dimSizes(self%vars(varIdx)%dimids(3))
             dimSize4=self%dimSizes(self%vars(varIdx)%dimids(4))
             allocate(values(dimSize1,dimSize2,dimSize3,dimSize4))
             call nf90ai_check( &
               nf90_get_var(self%ncid,varIdx,values) &
             )
             call nf90ai_check( nf90_close(self%ncid) )
             !print *, trim(self%vars(varIdx)%name), &
             !         self%vars(varIdx)%nDims, &
             !         varIdx,dimSize1,dimSize2,dimSize3
             !print *, ((((values(i,j,k,l),l=1,dimSize4),k=1,dimSize3), &
             !             j=1,dimSize2),i=1,dimSize1)
           else
             returnVal=-999
           endif
         else
           returnVal=-999
         endif

       end function nf90ai_get_var_4D

       function nf90ai_get_var_index(self,varName,varIdx) result(returnVal)

         type(nf90aiMeta),intent(in) :: self
         character(len=*), intent(in) :: varName
         integer,intent(inout) :: varIdx

         integer :: returnVal
         integer :: i

         returnVal=-999
         varIdx=-1
         if(varName .eq. NF90AI_GLOBAL_VAR_NAME) then
           varIdx=0
           returnVal=nf90_noerr
         else
           do i=1,self%nVars
             if (varName .eq. self%vars(i)%name) then
               varIdx=i
               returnVal=nf90_noerr
               exit
             end if
           end do
         endif

       end function nf90ai_get_var_index

       function nf90ai_get_var_att_text(self,varName,attName,valueString) result(returnVal)
         type(nf90aiMeta),intent(in) :: self
         character(len=*),intent(in) :: varName,attName
         character(len=NF90AI_MAX_CHAR_SIZE),intent(inout) :: valueString

         integer :: varIdx
         integer :: returnVal

         if(nf90ai_get_var_index(self,varName,varIdx) == nf90_noerr) then
           returnVal = nf90ai_get_att(self,varIdx,attName,valueString)
         else
           returnVal = -999
         end if

       end function nf90ai_get_var_att_text

       function nf90ai_get_var_att_1D(self,varName,attName,values) result(returnVal)
         type(nf90aiMeta),intent(in) :: self
         character(len=*),intent(in) :: varName,attName
         real,intent(inout),pointer :: values(:)

         integer :: varIdx
         integer :: returnVal

         nullify(values)

         if(nf90ai_get_var_index(self,varName,varIdx) == nf90_noerr) then
           returnVal = nf90ai_get_att(self,varIdx,attName,values)
         else
           returnVal = -999
         end if

       end function nf90ai_get_var_att_1D

     end module netcdfai_var

max_stars_repo_name: ScottLongmore/netCDF-fortran-API

max_issues_repo_name: ScottLongmore/netCDF-fortran-API

max_forks_repo_name: ScottLongmore/netCDF-fortran-API