text:
George B. Hoag ran a general merchandise store that opened in 1899?

His store was also the central station for the Farmers Telephone Company that connected Yolo and Solano counties.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE ML5_0_MP_BORN_AMPS_AND_WFS(P)
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     Computes all the AMP and WFS in quadruple precision for the
C     phase space point P(0:3,NEXTERNAL)
C
C     Process: g g > w- t b~ QCD<=2 QED<=1 [ virt = QCD ]
C
      IMPLICIT NONE
C
C     CONSTANTS
C
      INTEGER NBORNAMPS
      PARAMETER (NBORNAMPS=8)
      INTEGER    NLOOPAMPS, NCTAMPS
      PARAMETER (NLOOPAMPS=396, NCTAMPS=252)
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      INTEGER    NWAVEFUNCS
      PARAMETER (NWAVEFUNCS=28)
      INTEGER    NCOMB
      PARAMETER (NCOMB=48)
      REAL*16     ZERO
      PARAMETER (ZERO=0E0_16)
      COMPLEX*32 IMAG1
      PARAMETER (IMAG1=(0E0_16,1E0_16))

C
C     ARGUMENTS
C
      REAL*16 P(0:3,NEXTERNAL)
C
C     LOCAL VARIABLES
C
      INTEGER I,J,H
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
      DATA IC/NEXTERNAL*1/
C
C     FUNCTIONS
C
      LOGICAL ML5_0_IS_HEL_SELECTED
C
C     GLOBAL VARIABLES
C
      INCLUDE 'mp_coupl_same_name.inc'

      INTEGER NTRY
      LOGICAL CHECKPHASE,HELDOUBLECHECKED
      REAL*8 REF
      COMMON/ML5_0_INIT/NTRY,CHECKPHASE,HELDOUBLECHECKED,REF

      LOGICAL GOODHEL(NCOMB)
      LOGICAL GOODAMP(NLOOPAMPS,NCOMB)
      COMMON/ML5_0_FILTERS/GOODAMP,GOODHEL

      INTEGER HELPICKED
      COMMON/ML5_0_HELCHOICE/HELPICKED

      COMPLEX*32 AMP(NBORNAMPS,NCOMB)
      COMMON/ML5_0_MP_AMPS/AMP
      COMPLEX*16 DPAMP(NBORNAMPS,NCOMB)
      COMMON/ML5_0_AMPS/DPAMP
      COMPLEX*32 W(20,NWAVEFUNCS,NCOMB)
      COMMON/ML5_0_MP_WFS/W

      COMPLEX*32 AMPL(3,NCTAMPS)
      COMMON/ML5_0_MP_AMPL/AMPL

      COMPLEX*16 DPW(20,NWAVEFUNCS,NCOMB)
      COMMON/ML5_0_WFCTS/DPW

      COMPLEX*16 DPAMPL(3,NLOOPAMPS)
      LOGICAL S(NLOOPAMPS)
      COMMON/ML5_0_AMPL/DPAMPL,S

      INTEGER HELC(NEXTERNAL,NCOMB)
      COMMON/ML5_0_HELCONFIGS/HELC

      LOGICAL MP_DONE_ONCE
      COMMON/ML5_0_MP_DONE_ONCE/MP_DONE_ONCE

C     This array specify potential special requirements on the
C      helicities to
C     consider. POLARIZATIONS(0,0) is -1 if there is not such
C      requirement.
      INTEGER POLARIZATIONS(0:NEXTERNAL,0:5)
      COMMON/ML5_0_BEAM_POL/POLARIZATIONS

C     ----------
C     BEGIN CODE
C     ---------

      MP_DONE_ONCE=.TRUE.

C     To be on the safe side, we always update the MP params here.
C     It can be redundant as this routine can be called a couple of
C      times for the same PS point during the stability checks.
C     But it is really not time consuming and I would rather be safe.
      CALL MP_UPDATE_AS_PARAM()

      DO H=1,NCOMB
        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1)
     $   .AND.((CHECKPHASE.OR..NOT.HELDOUBLECHECKED).OR.GOODHEL(H))))
     $    THEN
C         Handle the possible requirement of specific polarizations
          IF ((.NOT.CHECKPHASE)
     $     .AND.HELDOUBLECHECKED.AND.POLARIZATIONS(0,0)
     $     .EQ.0.AND.(.NOT.ML5_0_IS_HEL_SELECTED(H))) THEN
            CYCLE
          ENDIF
          DO I=1,NEXTERNAL
            NHEL(I)=HELC(I,H)
          ENDDO
          CALL MP_VXXXXX(P(0,1),ZERO,NHEL(1),-1*IC(1),W(1,1,H))
          CALL MP_VXXXXX(P(0,2),ZERO,NHEL(2),-1*IC(2),W(1,2,H))
          CALL MP_VXXXXX(P(0,3),MDL_MW,NHEL(3),+1*IC(3),W(1,3,H))
          CALL MP_OXXXXX(P(0,4),MDL_MT,NHEL(4),+1*IC(4),W(1,4,H))
          CALL MP_IXXXXX(P(0,5),MDL_MB,NHEL(5),-1*IC(5),W(1,5,H))
          CALL MP_VVV1P0_1(W(1,1,H),W(1,2,H),GC_4,ZERO,ZERO,W(1,6,H))
          CALL MP_FFV2_1(W(1,4,H),W(1,3,H),GC_11,MDL_MB,ZERO,W(1,7,H))
C         Amplitude(s) for born diagram with ID 1
          CALL MP_FFV1_0(W(1,5,H),W(1,7,H),W(1,6,H),GC_5,AMP(1,H))
          CALL MP_FFV2_2(W(1,5,H),W(1,3,H),GC_11,MDL_MT,MDL_WT,W(1,8,H)
     $     )
C         Amplitude(s) for born diagram with ID 2
          CALL MP_FFV1_0(W(1,8,H),W(1,4,H),W(1,6,H),GC_5,AMP(2,H))
          CALL MP_FFV1_1(W(1,4,H),W(1,1,H),GC_5,MDL_MT,MDL_WT,W(1,9,H))
          CALL MP_FFV1_2(W(1,5,H),W(1,2,H),GC_5,MDL_MB,ZERO,W(1,10,H))
C         Amplitude(s) for born diagram with ID 3
          CALL MP_FFV2_0(W(1,10,H),W(1,9,H),W(1,3,H),GC_11,AMP(3,H))
C         Amplitude(s) for born diagram with ID 4
          CALL MP_FFV1_0(W(1,8,H),W(1,9,H),W(1,2,H),GC_5,AMP(4,H))
          CALL MP_FFV1_2(W(1,5,H),W(1,1,H),GC_5,MDL_MB,ZERO,W(1,11,H))
          CALL MP_FFV1_1(W(1,4,H),W(1,2,H),GC_5,MDL_MT,MDL_WT,W(1,12,H)
     $     )
C         Amplitude(s) for born diagram with ID 5
          CALL MP_FFV2_0(W(1,11,H),W(1,12,H),W(1,3,H),GC_11,AMP(5,H))
C         Amplitude(s) for born diagram with ID 6
          CALL MP_FFV1_0(W(1,11,H),W(1,7,H),W(1,2,H),GC_5,AMP(6,H))
C         Amplitude(s) for born diagram with ID 7
          CALL MP_FFV1_0(W(1,8,H),W(1,12,H),W(1,1,H),GC_5,AMP(7,H))
C         Amplitude(s) for born diagram with ID 8
          CALL MP_FFV1_0(W(1,10,H),W(1,7,H),W(1,1,H),GC_5,AMP(8,H))
          CALL MP_FFV1_2(W(1,5,H),W(1,6,H),GC_5,MDL_MB,ZERO,W(1,13,H))
C         Counter-term amplitude(s) for loop diagram number 9
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,13,H),W(1,7,H),R2_QQQ,R2_QQB
     $     ,AMPL(1,1))
          CALL MP_R2_QQ_2_0(W(1,13,H),W(1,7,H),UV_BMASS_1EPS,AMPL(2,2))
          CALL MP_R2_QQ_2_0(W(1,13,H),W(1,7,H),UV_BMASS,AMPL(1,3))
          CALL MP_FFV1P0_3(W(1,5,H),W(1,7,H),GC_5,ZERO,ZERO,W(1,14,H))
C         Counter-term amplitude(s) for loop diagram number 10
          CALL MP_R2_GG_1_R2_GG_2_0(W(1,6,H),W(1,14,H),R2_GGG_1
     $     ,R2_GGG_2,AMPL(1,4))
C         Counter-term amplitude(s) for loop diagram number 11
          CALL MP_FFV1_0(W(1,5,H),W(1,7,H),W(1,6,H),R2_GQQ,AMPL(1,5))
          CALL MP_FFV1_0(W(1,5,H),W(1,7,H),W(1,6,H),UV_GQQB_1EPS
     $     ,AMPL(2,6))
          CALL MP_FFV1_0(W(1,5,H),W(1,7,H),W(1,6,H),UV_GQQB_1EPS
     $     ,AMPL(2,7))
          CALL MP_FFV1_0(W(1,5,H),W(1,7,H),W(1,6,H),UV_GQQB_1EPS
     $     ,AMPL(2,8))
          CALL MP_FFV1_0(W(1,5,H),W(1,7,H),W(1,6,H),UV_GQQB_1EPS
     $     ,AMPL(2,9))
          CALL MP_FFV1_0(W(1,5,H),W(1,7,H),W(1,6,H),UV_GQQB_1EPS
     $     ,AMPL(2,10))
          CALL MP_FFV1_0(W(1,5,H),W(1,7,H),W(1,6,H),UV_GQQB_1EPS
     $     ,AMPL(2,11))
          CALL MP_FFV1_0(W(1,5,H),W(1,7,H),W(1,6,H),UV_GQQG_1EPS
     $     ,AMPL(2,12))
          CALL MP_FFV1_0(W(1,5,H),W(1,7,H),W(1,6,H),UV_GQQB,AMPL(1,13))
          CALL MP_FFV1_0(W(1,5,H),W(1,7,H),W(1,6,H),UV_GQQT,AMPL(1,14))
          CALL MP_FFV1_1(W(1,4,H),W(1,6,H),GC_5,MDL_MT,MDL_WT,W(1,15,H)
     $     )
C         Counter-term amplitude(s) for loop diagram number 13
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,8,H),W(1,15,H),R2_QQQ,R2_QQT
     $     ,AMPL(1,15))
          CALL MP_R2_QQ_2_0(W(1,8,H),W(1,15,H),UV_TMASS_1EPS,AMPL(2,16)
     $     )
          CALL MP_R2_QQ_2_0(W(1,8,H),W(1,15,H),UV_TMASS,AMPL(1,17))
          CALL MP_FFV1P0_3(W(1,8,H),W(1,4,H),GC_5,ZERO,ZERO,W(1,16,H))
C         Counter-term amplitude(s) for loop diagram number 14
          CALL MP_R2_GG_1_R2_GG_2_0(W(1,6,H),W(1,16,H),R2_GGG_1
     $     ,R2_GGG_2,AMPL(1,18))
C         Counter-term amplitude(s) for loop diagram number 15
          CALL MP_FFV1_0(W(1,8,H),W(1,4,H),W(1,6,H),R2_GQQ,AMPL(1,19))
          CALL MP_FFV1_0(W(1,8,H),W(1,4,H),W(1,6,H),UV_GQQB_1EPS
     $     ,AMPL(2,20))
          CALL MP_FFV1_0(W(1,8,H),W(1,4,H),W(1,6,H),UV_GQQB_1EPS
     $     ,AMPL(2,21))
          CALL MP_FFV1_0(W(1,8,H),W(1,4,H),W(1,6,H),UV_GQQB_1EPS
     $     ,AMPL(2,22))
          CALL MP_FFV1_0(W(1,8,H),W(1,4,H),W(1,6,H),UV_GQQB_1EPS
     $     ,AMPL(2,23))
          CALL MP_FFV1_0(W(1,8,H),W(1,4,H),W(1,6,H),UV_GQQB_1EPS
     $     ,AMPL(2,24))
          CALL MP_FFV1_0(W(1,8,H),W(1,4,H),W(1,6,H),UV_GQQB_1EPS
     $     ,AMPL(2,25))
          CALL MP_FFV1_0(W(1,8,H),W(1,4,H),W(1,6,H),UV_GQQG_1EPS
     $     ,AMPL(2,26))
          CALL MP_FFV1_0(W(1,8,H),W(1,4,H),W(1,6,H),UV_GQQB,AMPL(1,27))
          CALL MP_FFV1_0(W(1,8,H),W(1,4,H),W(1,6,H),UV_GQQT,AMPL(1,28))
C         Counter-term amplitude(s) for loop diagram number 17
          CALL MP_FFV2_0(W(1,13,H),W(1,4,H),W(1,3,H),R2_BXTW,AMPL(1,29)
     $     )
C         Counter-term amplitude(s) for loop diagram number 21
          CALL MP_FFV2_0(W(1,5,H),W(1,15,H),W(1,3,H),R2_BXTW,AMPL(1,30)
     $     )
          CALL MP_FFV2_1(W(1,9,H),W(1,3,H),GC_11,MDL_MB,ZERO,W(1,17,H))
C         Counter-term amplitude(s) for loop diagram number 22
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,10,H),W(1,17,H),R2_QQQ,R2_QQB
     $     ,AMPL(1,31))
          CALL MP_R2_QQ_2_0(W(1,10,H),W(1,17,H),UV_BMASS_1EPS,AMPL(2
     $     ,32))
          CALL MP_R2_QQ_2_0(W(1,10,H),W(1,17,H),UV_BMASS,AMPL(1,33))
          CALL MP_FFV2_2(W(1,10,H),W(1,3,H),GC_11,MDL_MT,MDL_WT,W(1,18
     $     ,H))
C         Counter-term amplitude(s) for loop diagram number 23
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,18,H),W(1,9,H),R2_QQQ,R2_QQT
     $     ,AMPL(1,34))
          CALL MP_R2_QQ_2_0(W(1,18,H),W(1,9,H),UV_TMASS_1EPS,AMPL(2,35)
     $     )
          CALL MP_R2_QQ_2_0(W(1,18,H),W(1,9,H),UV_TMASS,AMPL(1,36))
C         Counter-term amplitude(s) for loop diagram number 24
          CALL MP_FFV2_0(W(1,10,H),W(1,9,H),W(1,3,H),R2_BXTW,AMPL(1,37)
     $     )
C         Counter-term amplitude(s) for loop diagram number 25
          CALL MP_FFV1_0(W(1,5,H),W(1,17,H),W(1,2,H),R2_GQQ,AMPL(1,38))
          CALL MP_FFV1_0(W(1,5,H),W(1,17,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,39))
          CALL MP_FFV1_0(W(1,5,H),W(1,17,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,40))
          CALL MP_FFV1_0(W(1,5,H),W(1,17,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,41))
          CALL MP_FFV1_0(W(1,5,H),W(1,17,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,42))
          CALL MP_FFV1_0(W(1,5,H),W(1,17,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,43))
          CALL MP_FFV1_0(W(1,5,H),W(1,17,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,44))
          CALL MP_FFV1_0(W(1,5,H),W(1,17,H),W(1,2,H),UV_GQQG_1EPS
     $     ,AMPL(2,45))
          CALL MP_FFV1_0(W(1,5,H),W(1,17,H),W(1,2,H),UV_GQQB,AMPL(1,46)
     $     )
          CALL MP_FFV1_0(W(1,5,H),W(1,17,H),W(1,2,H),UV_GQQT,AMPL(1,47)
     $     )
C         Counter-term amplitude(s) for loop diagram number 27
          CALL MP_FFV1_0(W(1,8,H),W(1,9,H),W(1,2,H),R2_GQQ,AMPL(1,48))
          CALL MP_FFV1_0(W(1,8,H),W(1,9,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,49))
          CALL MP_FFV1_0(W(1,8,H),W(1,9,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,50))
          CALL MP_FFV1_0(W(1,8,H),W(1,9,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,51))
          CALL MP_FFV1_0(W(1,8,H),W(1,9,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,52))
          CALL MP_FFV1_0(W(1,8,H),W(1,9,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,53))
          CALL MP_FFV1_0(W(1,8,H),W(1,9,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,54))
          CALL MP_FFV1_0(W(1,8,H),W(1,9,H),W(1,2,H),UV_GQQG_1EPS
     $     ,AMPL(2,55))
          CALL MP_FFV1_0(W(1,8,H),W(1,9,H),W(1,2,H),UV_GQQB,AMPL(1,56))
          CALL MP_FFV1_0(W(1,8,H),W(1,9,H),W(1,2,H),UV_GQQT,AMPL(1,57))
          CALL MP_FFV1_1(W(1,9,H),W(1,2,H),GC_5,MDL_MT,MDL_WT,W(1,19,H)
     $     )
C         Counter-term amplitude(s) for loop diagram number 28
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,8,H),W(1,19,H),R2_QQQ,R2_QQT
     $     ,AMPL(1,58))
          CALL MP_R2_QQ_2_0(W(1,8,H),W(1,19,H),UV_TMASS_1EPS,AMPL(2,59)
     $     )
          CALL MP_R2_QQ_2_0(W(1,8,H),W(1,19,H),UV_TMASS,AMPL(1,60))
          CALL MP_FFV1_2(W(1,8,H),W(1,2,H),GC_5,MDL_MT,MDL_WT,W(1,20,H)
     $     )
C         Counter-term amplitude(s) for loop diagram number 29
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,20,H),W(1,9,H),R2_QQQ,R2_QQT
     $     ,AMPL(1,61))
          CALL MP_R2_QQ_2_0(W(1,20,H),W(1,9,H),UV_TMASS_1EPS,AMPL(2,62)
     $     )
          CALL MP_R2_QQ_2_0(W(1,20,H),W(1,9,H),UV_TMASS,AMPL(1,63))
C         Counter-term amplitude(s) for loop diagram number 31
          CALL MP_FFV2_0(W(1,5,H),W(1,19,H),W(1,3,H),R2_BXTW,AMPL(1,64)
     $     )
          CALL MP_FFV2_2(W(1,11,H),W(1,3,H),GC_11,MDL_MT,MDL_WT,W(1,21
     $     ,H))
C         Counter-term amplitude(s) for loop diagram number 35
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,21,H),W(1,12,H),R2_QQQ,R2_QQT
     $     ,AMPL(1,65))
          CALL MP_R2_QQ_2_0(W(1,21,H),W(1,12,H),UV_TMASS_1EPS,AMPL(2
     $     ,66))
          CALL MP_R2_QQ_2_0(W(1,21,H),W(1,12,H),UV_TMASS,AMPL(1,67))
          CALL MP_FFV2_1(W(1,12,H),W(1,3,H),GC_11,MDL_MB,ZERO,W(1,22,H)
     $     )
C         Counter-term amplitude(s) for loop diagram number 36
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,11,H),W(1,22,H),R2_QQQ,R2_QQB
     $     ,AMPL(1,68))
          CALL MP_R2_QQ_2_0(W(1,11,H),W(1,22,H),UV_BMASS_1EPS,AMPL(2
     $     ,69))
          CALL MP_R2_QQ_2_0(W(1,11,H),W(1,22,H),UV_BMASS,AMPL(1,70))
C         Counter-term amplitude(s) for loop diagram number 37
          CALL MP_FFV2_0(W(1,11,H),W(1,12,H),W(1,3,H),R2_BXTW,AMPL(1
     $     ,71))
C         Counter-term amplitude(s) for loop diagram number 38
          CALL MP_FFV1_0(W(1,21,H),W(1,4,H),W(1,2,H),R2_GQQ,AMPL(1,72))
          CALL MP_FFV1_0(W(1,21,H),W(1,4,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,73))
          CALL MP_FFV1_0(W(1,21,H),W(1,4,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,74))
          CALL MP_FFV1_0(W(1,21,H),W(1,4,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,75))
          CALL MP_FFV1_0(W(1,21,H),W(1,4,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,76))
          CALL MP_FFV1_0(W(1,21,H),W(1,4,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,77))
          CALL MP_FFV1_0(W(1,21,H),W(1,4,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,78))
          CALL MP_FFV1_0(W(1,21,H),W(1,4,H),W(1,2,H),UV_GQQG_1EPS
     $     ,AMPL(2,79))
          CALL MP_FFV1_0(W(1,21,H),W(1,4,H),W(1,2,H),UV_GQQB,AMPL(1,80)
     $     )
          CALL MP_FFV1_0(W(1,21,H),W(1,4,H),W(1,2,H),UV_GQQT,AMPL(1,81)
     $     )
C         Counter-term amplitude(s) for loop diagram number 40
          CALL MP_FFV1_0(W(1,11,H),W(1,7,H),W(1,2,H),R2_GQQ,AMPL(1,82))
          CALL MP_FFV1_0(W(1,11,H),W(1,7,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,83))
          CALL MP_FFV1_0(W(1,11,H),W(1,7,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,84))
          CALL MP_FFV1_0(W(1,11,H),W(1,7,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,85))
          CALL MP_FFV1_0(W(1,11,H),W(1,7,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,86))
          CALL MP_FFV1_0(W(1,11,H),W(1,7,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,87))
          CALL MP_FFV1_0(W(1,11,H),W(1,7,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,88))
          CALL MP_FFV1_0(W(1,11,H),W(1,7,H),W(1,2,H),UV_GQQG_1EPS
     $     ,AMPL(2,89))
          CALL MP_FFV1_0(W(1,11,H),W(1,7,H),W(1,2,H),UV_GQQB,AMPL(1,90)
     $     )
          CALL MP_FFV1_0(W(1,11,H),W(1,7,H),W(1,2,H),UV_GQQT,AMPL(1,91)
     $     )
          CALL MP_FFV1_2(W(1,11,H),W(1,2,H),GC_5,MDL_MB,ZERO,W(1,23,H))
C         Counter-term amplitude(s) for loop diagram number 41
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,23,H),W(1,7,H),R2_QQQ,R2_QQB
     $     ,AMPL(1,92))
          CALL MP_R2_QQ_2_0(W(1,23,H),W(1,7,H),UV_BMASS_1EPS,AMPL(2,93)
     $     )
          CALL MP_R2_QQ_2_0(W(1,23,H),W(1,7,H),UV_BMASS,AMPL(1,94))
          CALL MP_FFV1_1(W(1,7,H),W(1,2,H),GC_5,MDL_MB,ZERO,W(1,24,H))
C         Counter-term amplitude(s) for loop diagram number 42
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,11,H),W(1,24,H),R2_QQQ,R2_QQB
     $     ,AMPL(1,95))
          CALL MP_R2_QQ_2_0(W(1,11,H),W(1,24,H),UV_BMASS_1EPS,AMPL(2
     $     ,96))
          CALL MP_R2_QQ_2_0(W(1,11,H),W(1,24,H),UV_BMASS,AMPL(1,97))
C         Counter-term amplitude(s) for loop diagram number 44
          CALL MP_FFV2_0(W(1,23,H),W(1,4,H),W(1,3,H),R2_BXTW,AMPL(1,98)
     $     )
C         Counter-term amplitude(s) for loop diagram number 48
          CALL MP_FFV1_0(W(1,5,H),W(1,22,H),W(1,1,H),R2_GQQ,AMPL(1,99))
          CALL MP_FFV1_0(W(1,5,H),W(1,22,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,100))
          CALL MP_FFV1_0(W(1,5,H),W(1,22,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,101))
          CALL MP_FFV1_0(W(1,5,H),W(1,22,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,102))
          CALL MP_FFV1_0(W(1,5,H),W(1,22,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,103))
          CALL MP_FFV1_0(W(1,5,H),W(1,22,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,104))
          CALL MP_FFV1_0(W(1,5,H),W(1,22,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,105))
          CALL MP_FFV1_0(W(1,5,H),W(1,22,H),W(1,1,H),UV_GQQG_1EPS
     $     ,AMPL(2,106))
          CALL MP_FFV1_0(W(1,5,H),W(1,22,H),W(1,1,H),UV_GQQB,AMPL(1
     $     ,107))
          CALL MP_FFV1_0(W(1,5,H),W(1,22,H),W(1,1,H),UV_GQQT,AMPL(1
     $     ,108))
C         Counter-term amplitude(s) for loop diagram number 50
          CALL MP_FFV1_0(W(1,8,H),W(1,12,H),W(1,1,H),R2_GQQ,AMPL(1,109)
     $     )
          CALL MP_FFV1_0(W(1,8,H),W(1,12,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,110))
          CALL MP_FFV1_0(W(1,8,H),W(1,12,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,111))
          CALL MP_FFV1_0(W(1,8,H),W(1,12,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,112))
          CALL MP_FFV1_0(W(1,8,H),W(1,12,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,113))
          CALL MP_FFV1_0(W(1,8,H),W(1,12,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,114))
          CALL MP_FFV1_0(W(1,8,H),W(1,12,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,115))
          CALL MP_FFV1_0(W(1,8,H),W(1,12,H),W(1,1,H),UV_GQQG_1EPS
     $     ,AMPL(2,116))
          CALL MP_FFV1_0(W(1,8,H),W(1,12,H),W(1,1,H),UV_GQQB,AMPL(1
     $     ,117))
          CALL MP_FFV1_0(W(1,8,H),W(1,12,H),W(1,1,H),UV_GQQT,AMPL(1
     $     ,118))
C         Counter-term amplitude(s) for loop diagram number 51
          CALL MP_FFV1_0(W(1,18,H),W(1,4,H),W(1,1,H),R2_GQQ,AMPL(1,119)
     $     )
          CALL MP_FFV1_0(W(1,18,H),W(1,4,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,120))
          CALL MP_FFV1_0(W(1,18,H),W(1,4,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,121))
          CALL MP_FFV1_0(W(1,18,H),W(1,4,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,122))
          CALL MP_FFV1_0(W(1,18,H),W(1,4,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,123))
          CALL MP_FFV1_0(W(1,18,H),W(1,4,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,124))
          CALL MP_FFV1_0(W(1,18,H),W(1,4,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,125))
          CALL MP_FFV1_0(W(1,18,H),W(1,4,H),W(1,1,H),UV_GQQG_1EPS
     $     ,AMPL(2,126))
          CALL MP_FFV1_0(W(1,18,H),W(1,4,H),W(1,1,H),UV_GQQB,AMPL(1
     $     ,127))
          CALL MP_FFV1_0(W(1,18,H),W(1,4,H),W(1,1,H),UV_GQQT,AMPL(1
     $     ,128))
C         Counter-term amplitude(s) for loop diagram number 53
          CALL MP_FFV1_0(W(1,10,H),W(1,7,H),W(1,1,H),R2_GQQ,AMPL(1,129)
     $     )
          CALL MP_FFV1_0(W(1,10,H),W(1,7,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,130))
          CALL MP_FFV1_0(W(1,10,H),W(1,7,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,131))
          CALL MP_FFV1_0(W(1,10,H),W(1,7,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,132))
          CALL MP_FFV1_0(W(1,10,H),W(1,7,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,133))
          CALL MP_FFV1_0(W(1,10,H),W(1,7,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,134))
          CALL MP_FFV1_0(W(1,10,H),W(1,7,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,135))
          CALL MP_FFV1_0(W(1,10,H),W(1,7,H),W(1,1,H),UV_GQQG_1EPS
     $     ,AMPL(2,136))
          CALL MP_FFV1_0(W(1,10,H),W(1,7,H),W(1,1,H),UV_GQQB,AMPL(1
     $     ,137))
          CALL MP_FFV1_0(W(1,10,H),W(1,7,H),W(1,1,H),UV_GQQT,AMPL(1
     $     ,138))
C         Counter-term amplitude(s) for loop diagram number 56
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),R2_3GG,AMPL(1,139)
     $     )
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),UV_3GB_1EPS
     $     ,AMPL(2,140))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),UV_3GB_1EPS
     $     ,AMPL(2,141))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),UV_3GB_1EPS
     $     ,AMPL(2,142))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),UV_3GB_1EPS
     $     ,AMPL(2,143))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),UV_3GB_1EPS
     $     ,AMPL(2,144))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),UV_3GB_1EPS
     $     ,AMPL(2,145))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),UV_3GG_1EPS
     $     ,AMPL(2,146))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),UV_3GB,AMPL(1,147)
     $     )
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),UV_3GT,AMPL(1,148)
     $     )
C         Counter-term amplitude(s) for loop diagram number 59
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),R2_3GG,AMPL(1,149)
     $     )
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),UV_3GB_1EPS
     $     ,AMPL(2,150))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),UV_3GB_1EPS
     $     ,AMPL(2,151))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),UV_3GB_1EPS
     $     ,AMPL(2,152))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),UV_3GB_1EPS
     $     ,AMPL(2,153))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),UV_3GB_1EPS
     $     ,AMPL(2,154))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),UV_3GB_1EPS
     $     ,AMPL(2,155))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),UV_3GG_1EPS
     $     ,AMPL(2,156))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),UV_3GB,AMPL(1,157)
     $     )
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),UV_3GT,AMPL(1,158)
     $     )
C         Counter-term amplitude(s) for loop diagram number 62
          CALL MP_FFV1_0(W(1,5,H),W(1,24,H),W(1,1,H),R2_GQQ,AMPL(1,159)
     $     )
          CALL MP_FFV1_0(W(1,5,H),W(1,24,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,160))
          CALL MP_FFV1_0(W(1,5,H),W(1,24,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,161))
          CALL MP_FFV1_0(W(1,5,H),W(1,24,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,162))
          CALL MP_FFV1_0(W(1,5,H),W(1,24,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,163))
          CALL MP_FFV1_0(W(1,5,H),W(1,24,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,164))
          CALL MP_FFV1_0(W(1,5,H),W(1,24,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,165))
          CALL MP_FFV1_0(W(1,5,H),W(1,24,H),W(1,1,H),UV_GQQG_1EPS
     $     ,AMPL(2,166))
          CALL MP_FFV1_0(W(1,5,H),W(1,24,H),W(1,1,H),UV_GQQB,AMPL(1
     $     ,167))
          CALL MP_FFV1_0(W(1,5,H),W(1,24,H),W(1,1,H),UV_GQQT,AMPL(1
     $     ,168))
C         Counter-term amplitude(s) for loop diagram number 65
          CALL MP_FFV1_0(W(1,20,H),W(1,4,H),W(1,1,H),R2_GQQ,AMPL(1,169)
     $     )
          CALL MP_FFV1_0(W(1,20,H),W(1,4,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,170))
          CALL MP_FFV1_0(W(1,20,H),W(1,4,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,171))
          CALL MP_FFV1_0(W(1,20,H),W(1,4,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,172))
          CALL MP_FFV1_0(W(1,20,H),W(1,4,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,173))
          CALL MP_FFV1_0(W(1,20,H),W(1,4,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,174))
          CALL MP_FFV1_0(W(1,20,H),W(1,4,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,175))
          CALL MP_FFV1_0(W(1,20,H),W(1,4,H),W(1,1,H),UV_GQQG_1EPS
     $     ,AMPL(2,176))
          CALL MP_FFV1_0(W(1,20,H),W(1,4,H),W(1,1,H),UV_GQQB,AMPL(1
     $     ,177))
          CALL MP_FFV1_0(W(1,20,H),W(1,4,H),W(1,1,H),UV_GQQT,AMPL(1
     $     ,178))
          CALL MP_FFV1_1(W(1,12,H),W(1,1,H),GC_5,MDL_MT,MDL_WT,W(1,25
     $     ,H))
C         Counter-term amplitude(s) for loop diagram number 70
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,8,H),W(1,25,H),R2_QQQ,R2_QQT
     $     ,AMPL(1,179))
          CALL MP_R2_QQ_2_0(W(1,8,H),W(1,25,H),UV_TMASS_1EPS,AMPL(2
     $     ,180))
          CALL MP_R2_QQ_2_0(W(1,8,H),W(1,25,H),UV_TMASS,AMPL(1,181))
          CALL MP_FFV1_2(W(1,8,H),W(1,1,H),GC_5,MDL_MT,MDL_WT,W(1,26,H)
     $     )
C         Counter-term amplitude(s) for loop diagram number 71
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,26,H),W(1,12,H),R2_QQQ,R2_QQT
     $     ,AMPL(1,182))
          CALL MP_R2_QQ_2_0(W(1,26,H),W(1,12,H),UV_TMASS_1EPS,AMPL(2
     $     ,183))
          CALL MP_R2_QQ_2_0(W(1,26,H),W(1,12,H),UV_TMASS,AMPL(1,184))
C         Counter-term amplitude(s) for loop diagram number 73
          CALL MP_FFV2_0(W(1,5,H),W(1,25,H),W(1,3,H),R2_BXTW,AMPL(1
     $     ,185))
          CALL MP_FFV1_2(W(1,10,H),W(1,1,H),GC_5,MDL_MB,ZERO,W(1,27,H))
C         Counter-term amplitude(s) for loop diagram number 77
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,27,H),W(1,7,H),R2_QQQ,R2_QQB
     $     ,AMPL(1,186))
          CALL MP_R2_QQ_2_0(W(1,27,H),W(1,7,H),UV_BMASS_1EPS,AMPL(2
     $     ,187))
          CALL MP_R2_QQ_2_0(W(1,27,H),W(1,7,H),UV_BMASS,AMPL(1,188))
          CALL MP_FFV1_1(W(1,7,H),W(1,1,H),GC_5,MDL_MB,ZERO,W(1,28,H))
C         Counter-term amplitude(s) for loop diagram number 78
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,10,H),W(1,28,H),R2_QQQ,R2_QQB
     $     ,AMPL(1,189))
          CALL MP_R2_QQ_2_0(W(1,10,H),W(1,28,H),UV_BMASS_1EPS,AMPL(2
     $     ,190))
          CALL MP_R2_QQ_2_0(W(1,10,H),W(1,28,H),UV_BMASS,AMPL(1,191))
C         Counter-term amplitude(s) for loop diagram number 80
          CALL MP_FFV2_0(W(1,27,H),W(1,4,H),W(1,3,H),R2_BXTW,AMPL(1
     $     ,192))
C         Counter-term amplitude(s) for loop diagram number 84
          CALL MP_FFV1_0(W(1,5,H),W(1,28,H),W(1,2,H),R2_GQQ,AMPL(1,193)
     $     )
          CALL MP_FFV1_0(W(1,5,H),W(1,28,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,194))
          CALL MP_FFV1_0(W(1,5,H),W(1,28,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,195))
          CALL MP_FFV1_0(W(1,5,H),W(1,28,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,196))
          CALL MP_FFV1_0(W(1,5,H),W(1,28,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,197))
          CALL MP_FFV1_0(W(1,5,H),W(1,28,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,198))
          CALL MP_FFV1_0(W(1,5,H),W(1,28,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,199))
          CALL MP_FFV1_0(W(1,5,H),W(1,28,H),W(1,2,H),UV_GQQG_1EPS
     $     ,AMPL(2,200))
          CALL MP_FFV1_0(W(1,5,H),W(1,28,H),W(1,2,H),UV_GQQB,AMPL(1
     $     ,201))
          CALL MP_FFV1_0(W(1,5,H),W(1,28,H),W(1,2,H),UV_GQQT,AMPL(1
     $     ,202))
C         Counter-term amplitude(s) for loop diagram number 87
          CALL MP_FFV1_0(W(1,26,H),W(1,4,H),W(1,2,H),R2_GQQ,AMPL(1,203)
     $     )
          CALL MP_FFV1_0(W(1,26,H),W(1,4,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,204))
          CALL MP_FFV1_0(W(1,26,H),W(1,4,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,205))
          CALL MP_FFV1_0(W(1,26,H),W(1,4,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,206))
          CALL MP_FFV1_0(W(1,26,H),W(1,4,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,207))
          CALL MP_FFV1_0(W(1,26,H),W(1,4,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,208))
          CALL MP_FFV1_0(W(1,26,H),W(1,4,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,209))
          CALL MP_FFV1_0(W(1,26,H),W(1,4,H),W(1,2,H),UV_GQQG_1EPS
     $     ,AMPL(2,210))
          CALL MP_FFV1_0(W(1,26,H),W(1,4,H),W(1,2,H),UV_GQQB,AMPL(1
     $     ,211))
          CALL MP_FFV1_0(W(1,26,H),W(1,4,H),W(1,2,H),UV_GQQT,AMPL(1
     $     ,212))
C         Counter-term amplitude(s) for loop diagram number 117
          CALL MP_R2_GG_1_0(W(1,6,H),W(1,14,H),R2_GGQ,AMPL(1,213))
          CALL MP_R2_GG_1_0(W(1,6,H),W(1,14,H),R2_GGQ,AMPL(1,214))
          CALL MP_R2_GG_1_0(W(1,6,H),W(1,14,H),R2_GGQ,AMPL(1,215))
          CALL MP_R2_GG_1_0(W(1,6,H),W(1,14,H),R2_GGQ,AMPL(1,216))
C         Counter-term amplitude(s) for loop diagram number 118
          CALL MP_R2_GG_1_0(W(1,6,H),W(1,16,H),R2_GGQ,AMPL(1,217))
          CALL MP_R2_GG_1_0(W(1,6,H),W(1,16,H),R2_GGQ,AMPL(1,218))
          CALL MP_R2_GG_1_0(W(1,6,H),W(1,16,H),R2_GGQ,AMPL(1,219))
          CALL MP_R2_GG_1_0(W(1,6,H),W(1,16,H),R2_GGQ,AMPL(1,220))
C         Counter-term amplitude(s) for loop diagram number 119
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),R2_3GQ,AMPL(1,221)
     $     )
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),R2_3GQ,AMPL(1,222)
     $     )
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),R2_3GQ,AMPL(1,223)
     $     )
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),R2_3GQ,AMPL(1,224)
     $     )
C         Counter-term amplitude(s) for loop diagram number 120
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),R2_3GQ,AMPL(1,225)
     $     )
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),R2_3GQ,AMPL(1,226)
     $     )
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),R2_3GQ,AMPL(1,227)
     $     )
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),R2_3GQ,AMPL(1,228)
     $     )
C         Counter-term amplitude(s) for loop diagram number 123
          CALL MP_R2_GG_1_R2_GG_3_0(W(1,6,H),W(1,14,H),R2_GGQ,R2_GGB
     $     ,AMPL(1,229))
C         Counter-term amplitude(s) for loop diagram number 124
          CALL MP_R2_GG_1_R2_GG_3_0(W(1,6,H),W(1,16,H),R2_GGQ,R2_GGB
     $     ,AMPL(1,230))
C         Counter-term amplitude(s) for loop diagram number 125
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),R2_3GQ,AMPL(1,231)
     $     )
C         Counter-term amplitude(s) for loop diagram number 126
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),R2_3GQ,AMPL(1,232)
     $     )
C         Counter-term amplitude(s) for loop diagram number 129
          CALL MP_R2_GG_1_R2_GG_3_0(W(1,6,H),W(1,14,H),R2_GGQ,R2_GGT
     $     ,AMPL(1,233))
C         Counter-term amplitude(s) for loop diagram number 130
          CALL MP_R2_GG_1_R2_GG_3_0(W(1,6,H),W(1,16,H),R2_GGQ,R2_GGT
     $     ,AMPL(1,234))
C         Counter-term amplitude(s) for loop diagram number 131
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,14,H),R2_3GQ,AMPL(1,235)
     $     )
C         Counter-term amplitude(s) for loop diagram number 132
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,16,H),R2_3GQ,AMPL(1,236)
     $     )
C         Amplitude(s) for UVCT diagram with ID 135
          CALL MP_FFV1_0(W(1,5,H),W(1,7,H),W(1,6,H),GC_5,AMPL(2,237))
          AMPL(2,237)=AMPL(2,237)*(4.0D0*UVWFCT_G_1_1EPS+2.0D0
     $     *UVWFCT_B_0_1EPS)
C         Amplitude(s) for UVCT diagram with ID 136
          CALL MP_FFV1_0(W(1,5,H),W(1,7,H),W(1,6,H),GC_5,AMPL(1,238))
          AMPL(1,238)=AMPL(1,238)*(2.0D0*UVWFCT_G_2+1.0D0*UVWFCT_T_0
     $     +1.0D0*UVWFCT_B_0+2.0D0*UVWFCT_G_1)
C         Amplitude(s) for UVCT diagram with ID 137
          CALL MP_FFV1_0(W(1,8,H),W(1,4,H),W(1,6,H),GC_5,AMPL(2,239))
          AMPL(2,239)=AMPL(2,239)*(4.0D0*UVWFCT_G_1_1EPS+2.0D0
     $     *UVWFCT_B_0_1EPS)
C         Amplitude(s) for UVCT diagram with ID 138
          CALL MP_FFV1_0(W(1,8,H),W(1,4,H),W(1,6,H),GC_5,AMPL(1,240))
          AMPL(1,240)=AMPL(1,240)*(2.0D0*UVWFCT_G_2+1.0D0*UVWFCT_T_0
     $     +1.0D0*UVWFCT_B_0+2.0D0*UVWFCT_G_1)
C         Amplitude(s) for UVCT diagram with ID 139
          CALL MP_FFV2_0(W(1,10,H),W(1,9,H),W(1,3,H),GC_11,AMPL(2,241))
          AMPL(2,241)=AMPL(2,241)*(4.0D0*UVWFCT_G_1_1EPS+2.0D0
     $     *UVWFCT_B_0_1EPS)
C         Amplitude(s) for UVCT diagram with ID 140
          CALL MP_FFV2_0(W(1,10,H),W(1,9,H),W(1,3,H),GC_11,AMPL(1,242))
          AMPL(1,242)=AMPL(1,242)*(2.0D0*UVWFCT_G_2+1.0D0*UVWFCT_T_0
     $     +1.0D0*UVWFCT_B_0+2.0D0*UVWFCT_G_1)
C         Amplitude(s) for UVCT diagram with ID 141
          CALL MP_FFV1_0(W(1,8,H),W(1,9,H),W(1,2,H),GC_5,AMPL(2,243))
          AMPL(2,243)=AMPL(2,243)*(4.0D0*UVWFCT_G_1_1EPS+2.0D0
     $     *UVWFCT_B_0_1EPS)
C         Amplitude(s) for UVCT diagram with ID 142
          CALL MP_FFV1_0(W(1,8,H),W(1,9,H),W(1,2,H),GC_5,AMPL(1,244))
          AMPL(1,244)=AMPL(1,244)*(2.0D0*UVWFCT_G_2+1.0D0*UVWFCT_T_0
     $     +1.0D0*UVWFCT_B_0+2.0D0*UVWFCT_G_1)
C         Amplitude(s) for UVCT diagram with ID 143
          CALL MP_FFV2_0(W(1,11,H),W(1,12,H),W(1,3,H),GC_11,AMPL(2,245)
     $     )
          AMPL(2,245)=AMPL(2,245)*(4.0D0*UVWFCT_G_1_1EPS+2.0D0
     $     *UVWFCT_B_0_1EPS)
C         Amplitude(s) for UVCT diagram with ID 144
          CALL MP_FFV2_0(W(1,11,H),W(1,12,H),W(1,3,H),GC_11,AMPL(1,246)
     $     )
          AMPL(1,246)=AMPL(1,246)*(2.0D0*UVWFCT_G_2+1.0D0*UVWFCT_T_0
     $     +1.0D0*UVWFCT_B_0+2.0D0*UVWFCT_G_1)
C         Amplitude(s) for UVCT diagram with ID 145
          CALL MP_FFV1_0(W(1,11,H),W(1,7,H),W(1,2,H),GC_5,AMPL(2,247))
          AMPL(2,247)=AMPL(2,247)*(4.0D0*UVWFCT_G_1_1EPS+2.0D0
     $     *UVWFCT_B_0_1EPS)
C         Amplitude(s) for UVCT diagram with ID 146
          CALL MP_FFV1_0(W(1,11,H),W(1,7,H),W(1,2,H),GC_5,AMPL(1,248))
          AMPL(1,248)=AMPL(1,248)*(2.0D0*UVWFCT_G_2+1.0D0*UVWFCT_T_0
     $     +1.0D0*UVWFCT_B_0+2.0D0*UVWFCT_G_1)
C         Amplitude(s) for UVCT diagram with ID 147
          CALL MP_FFV1_0(W(1,8,H),W(1,12,H),W(1,1,H),GC_5,AMPL(2,249))
          AMPL(2,249)=AMPL(2,249)*(4.0D0*UVWFCT_G_1_1EPS+2.0D0
     $     *UVWFCT_B_0_1EPS)
C         Amplitude(s) for UVCT diagram with ID 148
          CALL MP_FFV1_0(W(1,8,H),W(1,12,H),W(1,1,H),GC_5,AMPL(1,250))
          AMPL(1,250)=AMPL(1,250)*(2.0D0*UVWFCT_G_2+1.0D0*UVWFCT_T_0
     $     +1.0D0*UVWFCT_B_0+2.0D0*UVWFCT_G_1)
C         Amplitude(s) for UVCT diagram with ID 149
          CALL MP_FFV1_0(W(1,10,H),W(1,7,H),W(1,1,H),GC_5,AMPL(2,251))
          AMPL(2,251)=AMPL(2,251)*(4.0D0*UVWFCT_G_1_1EPS+2.0D0
     $     *UVWFCT_B_0_1EPS)
C         Amplitude(s) for UVCT diagram with ID 150
          CALL MP_FFV1_0(W(1,10,H),W(1,7,H),W(1,1,H),GC_5,AMPL(1,252))
          AMPL(1,252)=AMPL(1,252)*(2.0D0*UVWFCT_G_2+1.0D0*UVWFCT_T_0
     $     +1.0D0*UVWFCT_B_0+2.0D0*UVWFCT_G_1)
C         Copy the qp wfs to the dp ones as they are used to setup the
C          CT calls.
          DO I=1,NWAVEFUNCS
            DO J=1,20
              DPW(J,I,H)=W(J,I,H)
            ENDDO
          ENDDO
C         Same for the counterterms amplitudes
          DO I=1,NCTAMPS
            DO J=1,3
              DPAMPL(J,I)=AMPL(J,I)
              S(I)=.TRUE.
            ENDDO
          ENDDO
          DO I=1,NBORNAMPS
            DPAMP(I,H)=AMP(I,H)
          ENDDO
        ENDIF
      ENDDO

      END


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      SUBROUTINE TABDIM ( LUREP , NSVICH, NSVOCH, NFLCH , NSTCH ,
     +                    LASWI , ICONF ,         NO_ACT, ACTLST,
     +                    ACTUSE, NBPR  , NO_INS, NO_INE, NO_OUS,
     +                    NO_OUE, NO_FLU, NO_STO, NO_DIS, NO_VEL,
     +                    NOINFO, NOWARN, IERROR)
!
!     Deltares     SECTOR WATERRESOURCES AND ENVIRONMENT
!
!     CREATED: nov -1992 by Jan van Beek
!
!     FUNCTION            : Reads the binary proces definition file
!
!     LOGICAL UNITNUMBERS : LUREP   - report file
!
!     SUBROUTINES CALLED  : SRSTOP, stops execution
!
!     PARAMETERS          :
!
!     NAME    KIND     LENGTH     FUNCT.  DESCRIPTION
!     ----    -----    ------     ------- -----------
!     LUREP   INTEGER       1     INPUT   Report file
!     NSVICH  INTEGER       1     INPUT   No of extra input variables for charon
!     NSVOCH  INTEGER       1     INPUT   No of extra output variables for charon
!     NFLCH   INTEGER       1     INPUT   No of fluxes for charon
!     NSTCH   INTEGER       1     INPUT   No of stochiometric terms for charon
!     LASWI   LOGICAL       1     INPUT   Active processes only switch
!     ICONF   INTEGER       1     INPUT   Requested configuration index
!     NO_ACT  INTEGER       1     INPUT   Number of activated processes
!     ACTLST  CHARACTER     *     INPUT   List of activated processes
!     ACTUSE  INTEGER       *     INPUT   activated processes status indicator
!     NBPR    INTEGER       1     OUTPUT  number of active processes
!     NO_INS  INTEGER       1     OUTPUT  Number of input var on segment
!     NO_INE  INTEGER       1     OUTPUT  Number of input var on exchange
!     NO_OUS  INTEGER       1     OUTPUT  Number of output var on segment
!     NO_OUE  INTEGER       1     OUTPUT  Number of output var on exchange
!     NO_FLU  INTEGER       1     OUTPUT  Number of fluxes
!     NO_STO  INTEGER       1     OUTPUT  Number of stochis
!     NO_DIS  INTEGER       1     OUTPUT  Number of dispersions
!     NO_VEL  INTEGER       1     OUTPUT  Number of velocities
!     NOINFO  INTEGER       1     IN/OUT  Cummulative information count
!     NOWARN  INTEGER       1     IN/OUT  Cummulative warning count
!     IERROR  INTEGER       1     OUTPUT  Error indicatior
!
!     IMPLICIT NONE for extra compiler checks
!
      use timers       !   performance timers

      IMPLICIT NONE
!
!     Declaration of arguments
!
      INTEGER        LUREP           , NSVICH          ,
     +               NSVOCH          , NFLCH           ,
     +               NSTCH           , ICONF           ,
     +               NO_ACT          , NBPR            ,
     +               NO_INS          , NO_INE          ,
     +               NO_OUS          , NO_OUE          ,
     +               NO_FLU          , NO_STO          ,
     +               NO_DIS          , NO_VEL          ,
     +               NOINFO          , NOWARN          ,
     +               IERROR
      INTEGER        ACTUSE(*)
      LOGICAL        LASWI
      CHARACTER*(*)          :: ACTLST(*)
!
!     Common declarations
!
      INCLUDE 'data.inc'
!
!     Local declarations
!
      INTEGER        INS   , INE   , OUS   , OUE   , FLU   ,
     +               STO   , DIS   , VEL   , IPROC , IOFFSE,
     +               NAANTA, IINPU , IITEM , JNDEX , IOUTP ,
     +               IOFFS2, NAANT2, IDISP , IVELO , IOUTF ,
     +               ISTOC , IFLUX , IPRCNF, IGET  , IACT  ,
     +               ILIC  , IC
      LOGICAL        SWIT2D
      integer(4) :: ithndl = 0
      if (timon) call timstrt( "tabdim", ithndl )
!
!     Some init
!
      SWIT2D = .FALSE.

!
!     first loop, determine maximal dimensions
!
      NO_INS = 0
      NO_INE = 0
      NO_OUS = 0
      NO_OUE = 0
      NO_FLU = 0
      NO_STO = 0
      NO_DIS = 0
      NO_VEL = 0
      NBPR = 0
      do iproc=1,nproc
!
!        Check if process is requested and licensed
!
         IF ( .NOT. LASWI ) THEN
            IPRCNF = (IPROC-1)*NCONF + ICONF
            IGET   = ICNPRO(IPRCNF)
         ELSE
            CALL ZOEK ( PROCID(IPROC), NO_ACT, ACTLST, 10    , IACT  )
            IF ( IACT .GT. 0 ) THEN
!              ILIC = 0
!              DO IC = 1 , NCONF
!                 IPRCNF = (IPROC-1)*NCONF + IC
!                 IF ( ICNPRO(IPRCNF).GT.0 .AND. LICONF(IC).GT.0 ) THEN
!                    ILIC = 1
!                 ENDIF
!              ENDDO
!              IF ( ILIC .EQ. 0 ) THEN
!                 WRITE(LUREP,*)
!    +               ' ERROR: no valid license for activated process:',
!    +               PROCID(IPROC)
!                 CALL SRSTOP(1)
!              ENDIF
               IGET = 1
               ACTUSE(IACT) = 1
            ELSE
               IGET = 0
            ENDIF
         ENDIF
         IF ( IGET .GT. 0 ) THEN
            ins = 0
            ine = 0
            ous = 0
            oue = 0
            flu = 0
            sto = 0
            dis = 0
            vel = 0
            NBPR = NBPR+1

!           write (LUREP,'(''Process: '',a10)') procid(iproc)
!           write (*,'(''+Process: '',a10)') procid(iproc)

!------------------------------------------------------------------------c
!           CONSTRUCT PROCESS
!------------------------------------------------------------------------c

!           Fill PDF structure, only dimensions

!           INPUT ITEMS ON SEGMENT LEVEL/EXCHANGE LEVEL

!           scan input items table for FIRST occurence of proces

            do iinpu = 1 , ninpu
               call zoek ( procid(iproc), 1, inpupr(iinpu), 10, jndex)
               if ( jndex .gt.0 ) then
                  if ( inpusx(iinpu) .eq. 1 ) then
!                     ins = ins + 1
                      ins = max(ins,inpunm(iinpu)) ! this is incorrect, because ins and ine are mixed, resulting in overestimation
                  else
!                     ine = ine + 1
                      ine = max(ine,inpunm(iinpu)) ! this is incorrect, because ins and ine are mixed, resulting in overestimation
                  endif
               endif
            enddo

!           OUTPUT ITEMS ON SEGMENT LEVEL/EXCHANGE LEVEL

            do ioutp = 1 , noutp
               call zoek ( procid(iproc), 1, outppr(ioutp),10 , jndex)
               if ( jndex .gt.0 ) then
                  call zoek ( outpit(ioutp), nitem, itemid, 10, iitem)
                  if ( iitem .le. 0 ) stop 'unknown ITEM'

!                 Find item properties and store in PDF structure

                  if ( outpsx(ioutp) .eq. 1 ) then
!                     ous = ous + 1
                      ous = max(ous,outpnm(ioutp)) ! this is incorrect, because ous and oue are mixed, resulting in overestimation of almost a factor 2
                  else
!                     oue = oue + 1
                      oue = max(oue,outpnm(ioutp)) ! this is incorrect, because ous and oue are mixed, resulting in overestimation of almost a factor 2

!                     SCAN VELO and DISP TABLES FOR LINES ASSOCIATED WITH CURRENT OUTPUT ITEM ON EXCHANGE LEVEL

                      if ( .not. swit2d ) then
                         do idisp = 1 , ndisp
                            call zoek ( itemid(iitem)  , 1, dispit(idisp), 10, jndex)
                            if ( jndex .gt. 0 ) then
                               dis = dis + 1
                            endif
                         enddo
                      endif

!                     scan velocity lines table

                      if ( .not. swit2d ) then
                         do ivelo = 1 , nvelo
                            call zoek ( itemid(iitem)  , 1, veloit(ivelo), 10, jndex)
                            if ( jndex .gt. 0 ) then
                               vel = vel + 1
                            endif
                         enddo
                      endif

!                     END of processing output item on exchange level!

                  endif

               endif

            enddo

!           FLUXES

            do ioutf = 1 , noutf
               call zoek ( procid(iproc), 1, outfpr(ioutf), 10, jndex)
               if ( jndex .gt.0 ) then

!                 flu = flu + 1
                  flu = max(flu,outfnm(ioutf))

!                 Process current row

!                 Lookup flux in item  table

                  call zoek ( outffl(ioutf), nitem, itemid, 10, iflux)
                  if ( iflux .le. 0 ) THEN
                     WRITE(LUREP,*) 'ERROR: unknown flux:',OUTFFL(IOUTF)
                     IERROR = 3
                     GOTO 900
                  endif

                  ! scan stochi table for lines associated with present flux

                  do istoc = 1 , nstoc
                     call zoek ( itemid(iflux), 1, stocfl(istoc), 10, jndex)
                     if ( jndex .gt.0 ) then
                        sto = sto + 1
                     endif
                  enddo

               endif
            enddo

!           reserve for extra input/output if process is charon

            if ( procfo(iproc) .eq. 'D40CHA' ) then
               ins = ins + nsvich
               ous = ous + nsvoch
               flu = flu + nflch
               sto = sto + nstch
            endif

            NO_INS = NO_INS + INS
            NO_INE = NO_INE + INE
            NO_OUS = NO_OUS + OUS
            NO_OUE = NO_OUE + OUE
            NO_FLU = NO_FLU + FLU
            NO_STO = NO_STO + STO
            NO_DIS = NO_DIS + DIS
            NO_VEL = NO_VEL + VEL

         ENDIF
!
      enddo
!
  900 CONTINUE
!
!
      if (timon) call timstop( ithndl )
      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


! { dg-do run }
!
! PR fortran/51055
! PR fortran/49110
! PR fortran/60334

subroutine test()
  implicit none
  integer :: i = 5
  character(len=:), allocatable :: s1
  character(len=:), pointer :: s2
  character(len=5), target :: fifeC = 'FIVEC'
  call sub(s1, i)
  if (len(s1) /= 5) call abort()
  if (s1 /= "ZZZZZ") call abort()
  s2 => subfunc()
  if (len(s2) /= 5) call abort()
  if (s2 /= "FIVEC") call abort()
  s1 = addPrefix(subfunc())
  if (len(s1) /= 7) call abort()
  if (s1 /= "..FIVEC") call abort()
contains
  subroutine sub(str,j)
    character(len=:), allocatable :: str
    integer :: j
    str = REPEAT("Z",j)
    if (len(str) /= 5) call abort()
    if (str /= "ZZZZZ") call abort()
  end subroutine sub
  function subfunc() result(res)
    character(len=:), pointer :: res
    res => fifec
    if (len(res) /= 5) call abort()
    if (res /= "FIVEC") call abort()
  end function subfunc
  function addPrefix(str) result(res)
    character(len=:), pointer :: str
    character(len=:), allocatable :: res
    res = ".." // str
  end function addPrefix
end subroutine test

program a
 character(len=:),allocatable :: s
 integer :: j=2
 s = repeat ('x', j)
 if (len(repeat(' ',j)) /= 2) call abort()
 if (repeat('y',j) /= "yy") call abort()
 if (len(s) /= 2) call abort()
 if (s /= "xx") call abort()
 call test()
end program a


text:
program P2410
    implicit none
    integer :: n, i, gap = 1
    integer, dimension(0:1000000) :: save
    read *, n
    save(0) = 1
    do while (gap <= n)
        do i = gap, n
            save(i) = mod(save(i) + save(i - gap), 1000000000)
        end do
        gap = lshift(gap, 1)
    end do
    write (*, '(I0)') save(n)
end program P2410

max_stars_repo_name: daily-boj/kiwiyou

max_issues_repo_name: daily-boj/kiwiyou

max_forks_repo_name: daily-boj/kiwiyou


text:
c	gasd calcula las presiones parciales

      subroutine gasd(t,pe,pg,pp)
      parameter (ncontr=28)
      real*4 alfai(ncontr),chi1(ncontr),chi2(ncontr),
     *u0(ncontr),u1(ncontr),u2(ncontr)
      real*4 du0,du1,du2
      real*8 pp(*),p(28)
      real*8 t,pe,pg,g1,g2,g3,g4,g5,thetad,g3i
      real*8 a,b,c,d,e,c1,c2,c3,f1,f2,f3,f4,f5,fe,phtot
      real*4 theta,cmol(91),dcmol(91)

      theta=5040./t
      thetad=5040.d0/t
      call molecb(theta,cmol,dcmol)
      g4=pe*10.**(cmol(1))
      g5=pe*10.**(cmol(2))

c	ahora calculo las funciones de particion u0,u1,u2 y sus derivadas
      do 5 i=1,ncontr
      		iii=i
5     		call neldatb(iii,0.,weight,alfai(i),chi1(i),chi2(i))
6     do 4 i=1,ncontr
      	  iii=i
	  t0=t
     	  call nelfctb(iii,t0,u0(iii),u1(iii),u2(iii),du0,du1,du2)
4     end do

      g2=sahadouble(thetad,chi1(1),u0(1),u1(1),pe)   ! p(h+)/p(h)
      g3=1./sahadouble(thetad,0.754,1.,u0(1),pe)     ! p(h-)/p(h)
      pp(10)=g3
      g1=0.
      do 1 i=2,ncontr
        a=sahadouble(thetad,chi1(i),u0(i),u1(i),pe) !a=n(FeII)/n(FeI)
        b=sahadouble(thetad,chi2(i),u1(i),u2(i),pe) !b=n(FeIII)/n(FeII)

	c=1.+a*(1.+b)                               !c=n(Fe)/n(FeI)
        p(i)=alfai(i)/c   ! p/ph' for neutral he,li, ... o sea n(FeI)/n(Htot)

1	g1=g1+p(i)*a*(1.+2.*b)   !numero total de e que no vienen de ionizar H

        pp(2)=p(2)
        pp(3)=p(6)
        pp(4)=p(11)
        pp(5)=p(12)

        a=1.+g2+g3
        b=2.*(1.+g2/g5*g4)
        c=g5
        d=g2-g3
        e=g2/g5*g4
        c1=c*b**2+a*d*b-e*a**2
        c2=2.*a*e-d*b+a*b*g1
        c3=-(e+b*g1)
        f1=0.5*c2/c1
        f1=-f1+sign(1.d0,c1)*sqrt(f1**2-c3/c1)
        f5=(1.-a*f1)/b
        f4=e*f5
        f3=g3*f1
        f2=g2*f1
        fe=f2-f3+f4+g1
	phtot=pe/fe

        if(f5.gt.1.e-4) goto 2
          const6=g5/pe*f1**2
          const7=f2-f3+g1

	  do 3 i=1,5
      	       f5=phtot*const6
	       f4=e*f5
	       fe=const7+f4
	       phtot=pe/fe
3	  continue

2       pg=pe*(1.+(f1+f2+f3+f4+f5+0.1014)/fe)

        pp(1)=f1   ! p(h)/p(h')
        pp(6)=f2   ! p(h+)/p(h')
        pp(7)=f5   ! p(h2)/p(h')
        pp(8)=fe   ! pe/p(h')
        pp(9)=pe/(1.38054e-16*t) ! n(e)=pe/kt

      return
      end

max_stars_repo_name: cdiazbas/SIRcode

max_issues_repo_name: cdiazbas/SIRcode

max_forks_repo_name: cdiazbas/SIRcode


text:
Stan Forbes served on the City Council from 4/20/1996 to 6/20/2000.  He ran for City Council three more times, each time unsuccessfully, in 2000, 2004 and again in 2006.  Prior to that he served as Trustee for the Davis JUSD School Board from 1985 to 1996.  Stan splits his time between almond farming and running the Avid Reader in Sacramento.  His exwife, Alzada Knickerbocker, owns the The Avid Reader Avid Reader bookstore in downtown downtown Davis.

In http://www.californiaaggie.com/media/paper981/news/2004/02/20/FrontPage/Forbes.Suggests.Woodland.For.Student.Housing1316592.shtml?norewrite200603172010&sourcedomainwww.californiaaggie.com February 2004 Forbes suggested students could move to Woodland to find affordable housing.

Stan was a registered Republican until April 26th, 2006 according to James Schwab.  Currently Forbes is registered as decline to state.

Source:

McEwing, Shearon, 13 Oct. 1993, Forbes brings 8 years of experience to third campaign for school board, Davis Enterprise, pg. A1

20070712 13:37:24 nbsp He had my vote Users/StevenDaubert


!+++++++++++++++++++++++++++++++++++++++++++++++++++
! Test programe to make use of Alex' nml.f90 module
!+++++++++++++++++++++++++++++++++++++++++++++++++++
program example
  use ioparams, only: group1_t, control_t, &
    read_nml, write_nml, &
    parse_command_argument, print_help

  implicit none
  type(group1_t) :: par1
  type(control_t) :: ctr
  integer :: i, iostat
  character(len=256) :: arg

  ! read namelist
  open(88, file="namelist.nml")
  call read_nml(88, par1)
  call read_nml(88, ctr)
  close(88)

  ! parse command-line arguments
  i = 1
  do while(i <= command_argument_count())
    call parse_command_argument(par1, i, iostat)
    if (iostat==0) cycle
    call parse_command_argument(ctr, i, iostat, arg=arg)
    if (iostat == -1) stop  ! help
    ! if (arg == "--help".or.arg == "-h") stop
    if (iostat/=0) then
      write(*,'("Parameter not found: ",A20,". Type -h for help")') trim(arg)
      stop
    endif
  enddo

  write(*,*)
  if (ctr%print_nml) then
    write(*,*) "-------------------------------------"
    write(*,*) "Print namelist : "
    write(*,*) "-------------------------------------"
    call write_nml(6, par1)
    write(*,*) "-------------------------------------"
  endif
  write(*,*)
  write(*,*) "Nice ! Call with -h flag to get help on existing parameters."
  write(*,*)

end program


text:
MODULE neo_spline
! Spline arrays
  USE neo_precision
  INTEGER, PARAMETER       ::   mt = 1
  INTEGER, PARAMETER       ::   mp = 1
  INTEGER                  ::   theta_ind, phi_ind
  INTEGER                  ::   ierr

  REAL(kind=dp) ::  theta_d, phi_d

! Spline array for modb
  REAL(kind=dp),    DIMENSION(:,:,:,:), ALLOCATABLE :: b_spl
! Spline array for geodesic curviture
  REAL(kind=dp),    DIMENSION(:,:,:,:), ALLOCATABLE :: k_spl
! Spline array for sqrg11
  REAL(kind=dp),    DIMENSION(:,:,:,:), ALLOCATABLE :: g_spl
! Spline array for parallel derivative
  REAL(kind=dp),    DIMENSION(:,:,:,:), ALLOCATABLE :: p_spl
! Spline array for quasi-toroidal phi component of b
  REAL(kind=dp),    DIMENSION(:,:,:,:), ALLOCATABLE :: q_spl
END MODULE neo_spline

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
! RUN: %python %S/test_symbols.py %s %flang_fc1
! Test host association in internal subroutine of main program.

!DEF: /main MainProgram
program main
 !DEF: /main/x ObjectEntity INTEGER(4)
 integer x
 !DEF: /main/s (Subroutine) Subprogram
 call s
contains
 !REF: /main/s
 subroutine s
  !DEF: /main/s/y (Implicit) ObjectEntity REAL(4)
  !DEF: /main/s/x HostAssoc INTEGER(4)
  y = x
 contains
  !DEF: /main/s/s2 (Subroutine) Subprogram
  subroutine s2
   !DEF: /main/s/s2/z (Implicit) ObjectEntity REAL(4)
   !DEF: /main/s/s2/x HostAssoc INTEGER(4)
   z = x
  end subroutine
 end subroutine
end program

!DEF: /s (Subroutine) Subprogram
subroutine s
 !DEF: /s/x ObjectEntity REAL(4)
 real x(100, 100)
 !DEF: /s/s1 (Subroutine) Subprogram
 call s1
contains
 !REF: /s/s1
  subroutine s1
    !DEF: /s/s1/x HostAssoc REAL(4)
    print *, x(10, 10)
  end subroutine
end subroutine

!DEF: /sb (Subroutine) Subprogram
subroutine sb
 !DEF: /sb/x TARGET ObjectEntity REAL(4)
 real, target :: x
 !DEF: /sb/s1 (Subroutine) Subprogram
 call s1
contains
 !REF: /sb/s1
 subroutine s1
  !DEF: /sb/s1/p POINTER ObjectEntity REAL(4)
  real, pointer :: p
  !REF: /sb/s1/p
  !DEF: /sb/s1/x TARGET HostAssoc REAL(4)
  p => x
 end subroutine
end subroutine

! Test host associated symbols are also created for symbols that are use
! associated in the host.

!DEF: /m1 Module
module m1
 !DEF: /m1/x PUBLIC ObjectEntity REAL(4)
 real x(100,100)
 !DEF: /m1/x_target PUBLIC, TARGET ObjectEntity REAL(4)
 real, target :: x_target
end module

!DEF: /s_use (Subroutine) Subprogram
subroutine s_use
 !REF: /m1
 use :: m1
 !DEF: /s_use/x Use REAL(4)
 print *, x
 !DEF: /s_use/s1 (Subroutine) Subprogram
 call s1
contains
 !REF: /s_use/s1
 subroutine s1
  !DEF: /s_use/s1/x HostAssoc REAL(4)
  print *, x(10,10)
 end subroutine
end subroutine

!DEF: /sb_use (Subroutine) Subprogram
subroutine sb_use
 !REF: /m1
 use :: m1
 !DEF: /sb_use/x_target TARGET Use REAL(4)
 print *, x_target
 !DEF: /sb_use/s1 (Subroutine) Subprogram
 call s1
contains
 !REF: /sb_use/s1
 subroutine s1
  !DEF: /sb_use/s1/p POINTER ObjectEntity REAL(4)
  real, pointer :: p
  !REF: /sb_use/s1/p
  !DEF: /sb_use/s1/x_target TARGET HostAssoc REAL(4)
  p => x_target
 end subroutine
end subroutine

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm


!----------------------------------------------------------------------------
!BOP
!
! !MODULE:  wrf_binary_interface
!
! !DESCRIPTION: This stub provides the default interfaces to the
!               WRF binary capability in GSI.
!
! !REVISION HISTORY:
!
!  08Jun2016 Mahajan - Initial code
!
!EOP
!-------------------------------------------------------------------------
module get_wrf_binary_interface_mod
use abstract_get_wrf_binary_interface_mod
  type, extends(abstract_get_wrf_binary_interface_class) :: get_wrf_binary_interface_class
  contains
    procedure, pass(this) :: convert_binary_mass => convert_binary_mass_dummy
    procedure, pass(this) :: convert_binary_nmm => convert_binary_nmm_dummy
    procedure, pass(this) :: convert_nems_nmmb => convert_nems_nmmb_dummy
  end type get_wrf_binary_interface_class
contains

  subroutine convert_binary_mass_dummy(this)
    implicit none
    class(get_wrf_binary_interface_class), intent(inout) :: this
  end subroutine convert_binary_mass_dummy

  subroutine convert_binary_nmm_dummy(this,update_pint,ctph0,stph0,tlm0)
    use kinds, only: r_kind
    implicit none
    class(get_wrf_binary_interface_class), intent(inout) :: this
    logical     ,intent(inout) :: update_pint
    real(r_kind),intent(  out) :: ctph0,stph0,tlm0
  end subroutine convert_binary_nmm_dummy

  subroutine convert_nems_nmmb_dummy(this,update_pint,ctph0,stph0,tlm0)
    use kinds, only: r_kind
    implicit none
    class(get_wrf_binary_interface_class), intent(inout) :: this
    logical     ,intent(inout) :: update_pint
    real(r_kind),intent(  out) :: ctph0,stph0,tlm0
  end subroutine convert_nems_nmmb_dummy

end module get_wrf_binary_interface_mod


text:
*###[ ffzli2:
	subroutine ffzli2(zdilog,zlog,cx,ier)
***#[*comment:***********************************************************
*									*
*	Computes the dilogarithm (Li2, Sp) for any (complex) cx 	*
*	to a precision precc.  It assumes that cx is already in the	*
*	area |cx|<=1, Re(cx)<=1/2.  As it is available it also returns	*
*	log(1-cx) = zlog.						*
*									*
*	Input:	cx	(complex)					*
*									*
*	Output:	zdilog	(complex) Li2(cx)				*
*		zlog	(complex) log(1-cx) = -Li1(cx)			*
*		ier	(integer) 0=OK,1=num,2=err			*
*									*
*	Calls:	log,zfflo1,(d/a)imag,real/dble				*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments
*
	integer ier
	DOUBLE COMPLEX cx,zlog,zdilog
*
*	local variables
*
	DOUBLE PRECISION xprec,bdn02,bdn05,bdn10,bdn15,
     +		xi,xr,xdilog,xlog,absc,xa,a,ffbnd
	DOUBLE COMPLEX cc,cz,cz2,zfflo1
	external ffbnd,zfflo1
	save xprec,bdn02,bdn05,bdn10,bdn15
*
*	common blocks
*
	include 'ff.h'
*
*	statement function
*
	absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
*  #] declarations:
*  #[ initialisations:
	data xprec /-1D0/
	if ( xprec .ne. precc ) then
	    xprec = precc
	    bdn02 = ffbnd(1,2,bf)
	    bdn05 = ffbnd(1,5,bf)
	    bdn10 = ffbnd(1,10,bf)
	    bdn15 = ffbnd(1,15,bf)
*		we don't have bf(21) ...
	endif
*  #] initialisations:
*  #[ exceptional cases:
	xi  = DIMAG(cx)
	xr  = DBLE(cx)
	if ( xi .eq. 0) then
	    call ffxli2(xdilog,xlog,xr,ier)
	    zdilog = xdilog
	    zlog = xlog
	    return
	endif
	xa = abs(xi) + abs(xr)
	if ( xa .lt. precc ) then
	    zdilog = cx
	    zlog = -cx
	    return
	endif
*  #] exceptional cases:
*  #[ get log,dilog:
	if ( xa .lt. xloss**2 ) then
	    zlog = zfflo1(cx,ier)
	else
	    zlog = log(1-cx)
	endif
	cz = -zlog
	if ( absc(cz) .lt. xclog2 ) then
	    zdilog = cz
	else
	cz2 = cz*cz
	a = xa**2
	if ( a .gt. bdn15 ) then
	    zdilog = cz2*(DBLE(bf(16)) + cz2*(DBLE(bf(17))
     +		   + cz2*(DBLE(bf(18)) + cz2*(DBLE(bf(19))
     +		   + cz2*(DBLE(bf(20)))))))
	else
	    zdilog = 0
	endif
	if ( a .gt. bdn10 ) then
	    zdilog = cz2*(DBLE(bf(11)) + cz2*(DBLE(bf(12))
     +		   + cz2*(DBLE(bf(13)) + cz2*(DBLE(bf(14))
     +		   + cz2*(DBLE(bf(15)) + zdilog)))))
	endif
	if ( a .gt. bdn05 ) then
	    zdilog = cz2*(DBLE(bf(6)) + cz2*(DBLE(bf(7))
     +		   + cz2*(DBLE(bf(8)) + cz2*(DBLE(bf(9))
     +		   + cz2*(DBLE(bf(10)) + zdilog)))))
	endif
	if ( a .gt. bdn02 ) then
	    zdilog = cz2*(DBLE(bf(3)) + cz2*(DBLE(bf(4))
     +		   + cz2*(DBLE(bf(5)) + zdilog)))
	endif
*	watch the powers of z.
	zdilog = cz + cz2*(DBLE(bf(1)) + cz*(DBLE(bf(2)) + zdilog))
	endif
*  #] get log,dilog:
*###] ffzli2:
	end

*###[ ffzzdl:
	subroutine ffzzdl(zdilog,ipi12,zlog,cx,ier)
***#[*comment:***************************************************
*								*
*	Computes the dilogarithm (Li2, Sp) for any (complex) cx *
*	to about 15 significant figures. This can be improved	*
*	by adding more of the bf's. For real cx > 1 an error is	*
*	generated as the imaginary part is undefined then.	*
*	For use in ffcdbd zlog = log(1-cx) is also calculated	*
*								*
*	Input:	cx	(complex)				*
*								*
*	Output:	zdilog	(complex) Li2(cx) mod factors pi^2/12	*
*		ipi12	(integer) these factors			*
*		zlog	(complex) log(1-cx)			*
*								*
*	Calls:	log,zfflo1,(d/a)imag,real/dble			*
*								*
***#]*comment:***************************************************
*  #[ declarations:
	implicit none
*
*	arguments
*
	integer ipi12,ier
	DOUBLE COMPLEX zdilog,zlog,cx
*
*	local variables
*
	integer jsgn
	DOUBLE PRECISION xprec,bdn02,bdn05,bdn10,bdn15,
     +		xi,xr,s1,s2,xa,a,absc,ffbnd
	DOUBLE COMPLEX cfact,cx1,cy,cz,cz2,zfflo1,c
	external ffbnd,zfflo1
	save xprec,bdn02,bdn05,bdn10,bdn15
*
*	common blocks
*
	include 'ff.h'
*
*	statement function
*
	absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
*  #] declarations:
*  #[ initialisations:
	data xprec /-1D0/
	if ( xprec .ne. precc ) then
	    xprec = precc
	    bdn02 = ffbnd(1,2,bf)
	    bdn05 = ffbnd(1,5,bf)
	    bdn10 = ffbnd(1,10,bf)
	    bdn15 = ffbnd(1,15,bf)
	endif
*  #] initialisations:
*  #[ exceptional cases:
	xi  = DIMAG(cx)
	xr  = DBLE(cx)
	if ( xi .eq. 0 ) then
	    if ( xr .gt. 1 ) call fferr(31,ier)
	    call ffzxdl(zdilog,ipi12,zlog,xr,1,ier)
	    return
	endif
	if ( abs(xi) .lt. xalog2 ) then
	    s1 = 0
	else
	    s1 = xi**2
	endif
	if ( abs(xr) .lt. xalog2 ) then
	    s2 = 0
	else
	    s2 = xr**2
	endif
	xa = sqrt(s1 + s2)
	if ( xa .lt. precc ) then
	    zdilog = cx
	    zlog = -cx
	    ipi12 = 0
	    return
	endif
*  #] exceptional cases:
*  #[ transform to |x|<1, Re(x) < 0.5:
	if ( xr .le. .5D0) then
	    if (xa .gt. 1) then
		if ( 1/xa .lt. xalogm ) then
		    cfact = 0
		elseif ( 1/xa .lt. xclogm ) then
		    cx1 = cx*DBLE(1/xa)
		    cfact = log(-cx1) + log(DBLE(xa))
		else
		    cfact = log(-cx)
		endif
		cy = - cfact**2/2
		ipi12 = -2
		if ( xa*xloss**2 .gt. 1) then
		    if ( 1/xa .lt. xclogm ) then
			cx1 = cx*DBLE(1/xa)
			cx1 = 1/cx1
			cx1 = cx1*DBLE(1/xa)
		    else
			cx1 = 1/cx
		    endif
		    cz = -zfflo1(cx1,ier)
		else
		    cz = -log(1-1/cx)
		endif
		zlog = log(1-cx)
		jsgn = -1
	    else
		cy = 0
		ipi12 = 0
		if ( xa .lt. xloss**2 ) then
		    zlog = zfflo1(cx,ier)
		else
		    zlog = log(1-cx)
		endif
		cz = -zlog
		jsgn = 1
	    endif
	else
	    if (xa .le. sqrt(2*xr)) then
		cz = -log(cx)
		if ( abs(xr-1) + abs(xi) .lt. xclogm ) then
		    cy = 0
		else
		    zlog = log(1-cx)
		    cy = cz*zlog
		endif
		ipi12 = 2
		jsgn = -1
	    else
		if ( 1/xa .lt. xalogm ) then
		    cfact = 0
		elseif ( 1/xa .lt. xclogm ) then
		    cx1 = cx*DBLE(1/xa)
		    cfact = log(-cx1) + log(DBLE(xa))
		else
		    cfact = log(-cx)
		endif
		cy = - cfact**2/2
		ipi12 = -2
		if ( xa*xloss .gt. 1) then
		    if ( 1/xa .lt. xclogm ) then
			cx1 = cx*DBLE(1/xa)
			cx1 = 1/cx1
			cx1 = cx1*DBLE(1/xa)
		    else
			cx1 = 1/cx
		    endif
		    cz = -zfflo1(cx1,ier)
		else
		    cz = -log(1-1/cx)
		endif
		zlog = log(1-cx)
		jsgn = -1
	    endif
	endif
*  #] transform to |x|<1, Re(x) < 0.5:
*  #[ get dilog:
	if ( absc(cz) .lt. xclogm ) then
	    zdilog = cz
	else
	cz2 = cz*cz
	a = DBLE(cz)**2 + DIMAG(cz)**2
	if ( a .gt. bdn15 ) then
	    zdilog = cz2*(DBLE(bf(16)) + cz2*(DBLE(bf(17))
     +		   + cz2*(DBLE(bf(18)) + cz2*(DBLE(bf(19))
     +		   + cz2*(DBLE(bf(20)))))))
	else
	    zdilog = 0
	endif
	if ( a .gt. bdn10 ) then
	    zdilog = cz2*(DBLE(bf(11)) + cz2*(DBLE(bf(12))
     +		   + cz2*(DBLE(bf(13)) + cz2*(DBLE(bf(14))
     +		   + cz2*(DBLE(bf(15)) + zdilog)))))
	endif
	if ( a .gt. bdn05 ) then
	    zdilog = cz2*(DBLE(bf(6)) + cz2*(DBLE(bf(7))
     +		   + cz2*(DBLE(bf(8)) + cz2*(DBLE(bf(9))
     +		   + cz2*(DBLE(bf(10)) + zdilog)))))
	endif
	if ( a .gt. bdn02 ) then
	    zdilog = cz2*(DBLE(bf(3)) + cz2*(DBLE(bf(4))
     +		   + cz2*(DBLE(bf(5)) + zdilog)))
	endif
*	watch the powers of z.
	zdilog = cz + cz2*(DBLE(bf(1)) + cz*(DBLE(bf(2)) + zdilog))
	endif
	if(jsgn.eq.1)then
	    zdilog =  zdilog + cy
	else
	    zdilog = -zdilog + cy
	endif
*  #] get dilog:
*###] ffzzdl:
	end

*###[ zfflog:
	DOUBLE COMPLEX function zfflog(cx,ieps,cy,ier)
***#[*comment:***********************************************************
*									*
*	Calculate the complex logarithm of cx.  The following cases	*
*	are treted separately:						*
*		|cx| too small:		give warning and return 0	*
*					(for Absoft, Apollo DN300)	*
*		Im(cx) = 0, Re(cx) < 0:	take sign according to ieps	*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
*
*	arguments
*
	implicit none
	integer ieps,ier
	DOUBLE COMPLEX cx,cy
*
*	local variables
*
	DOUBLE COMPLEX c,ctroep
	DOUBLE PRECISION absc,xa,xlog1p
*
*	common blocks, statement function
*
	include 'ff.h'
	absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
*  #] declarations:
*  #[ calculations:
	xa = absc(cx)
	if ( xa .lt. xalogm ) then
	    if ( cx .ne. 0 ) call fferr(23,ier)
	    zfflog = 0
	elseif ( DBLE(cx) .lt. 0 .and. DIMAG(cx) .eq. 0 ) then
*     +		 abs(DIMAG(cx)) .lt. precc*abs(DBLE(cx)) ) then
	    xlog1p = log(-DBLE(cx))
*	    checked imaginary parts 19-May-1988
	    if ( abs(ieps) .eq. 1 ) then
		if ( ieps*DBLE(cy) .lt. 0 ) then
		    zfflog = DCMPLX(xlog1p,-pi)
		elseif ( ieps*DBLE(cy) .gt. 0 ) then
		    zfflog = DCMPLX(xlog1p,pi)
		else
		    call fferr(51,ier)
		    zfflog = DCMPLX(xlog1p,pi)
		endif
	    elseif ( ieps .ge. 2 .and. ieps .le. 3 ) then
		zfflog = DCMPLX(xlog1p,-pi)
	    elseif ( ieps .le. -2 .and. ieps .ge. -3 ) then
		zfflog = DCMPLX(xlog1p,pi)
	    else
		call fferr(51,ier)
		zfflog = DCMPLX(xlog1p,pi)
	    endif
	elseif ( xa .lt. xclogm .or. 1/xa .lt. xclogm ) then
	    ctroep = cx*DBLE(1/xa)
	    zfflog = log(ctroep) + DBLE(log(xa))
	else
*	    print *,'zfflog: neem log van ',cx
	    zfflog = log(cx)
	endif
*  #] calculations:
*###] zfflog:
	end

*###[ zfflo1:
	DOUBLE COMPLEX function zfflo1(cx,ier)
***#[*comment:***************************************************
*	calculates log(1-x) for |x|<.14 in a faster way to ~15	*
*	significant figures.					*
***#]*comment:***************************************************
*  #[ declarations:
	implicit none
	integer ier
	DOUBLE COMPLEX cx,c,zfflog
	DOUBLE PRECISION xprec,bdn01,bdn05,bdn10,bdn15,bdn19,
     +		absc,xa,ffbnd
	external zfflog,ffbnd
	save xprec,bdn01,bdn05,bdn10,bdn15,bdn19
	include 'ff.h'
	absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
*  #] declarations:
*  #[ initialisations:
	data xprec /-1D0/
	if ( precc .ne. xprec ) then
	    xprec = precc
*	    determine the boundaries for 1,5,10,15 terms
	    bdn01 = ffbnd(1,1,xninv)
	    bdn05 = ffbnd(1,5,xninv)
	    bdn10 = ffbnd(1,10,xninv)
	    bdn15 = ffbnd(1,15,xninv)
	    bdn19 = ffbnd(1,19,xninv)
	endif
*  #] initialisations:
*  #[ calculations:
	xa = absc(cx)
	if ( xa .gt. bdn19 ) then
	    c = cx-1
	    xa = absc(c)
	    zfflo1 = zfflog(1-cx,0,czero,ier)
	    return
	endif
	if ( xa .gt. bdn15 ) then
	    zfflo1 = cx*( DBLE(xninv(16)) + cx*( DBLE(xninv(17))
     +		   + cx*( DBLE(xninv(18)) + cx*( DBLE(xninv(19))
     +		   + cx*( DBLE(xninv(20)) )))))
	else
	    zfflo1 = 0
	endif
	if ( xa .gt. bdn10 ) then
	    zfflo1 = cx*( DBLE(xninv(11)) + cx*( DBLE(xninv(12))
     +		   + cx*( DBLE(xninv(13)) + cx*( DBLE(xninv(14))
     +		   + cx*( DBLE(xninv(15)) + zfflo1 )))))
	endif
	if ( xa .gt. bdn05 ) then
	    zfflo1 = cx*( DBLE(xninv(6)) + cx*( DBLE(xninv(7))
     +		   + cx*( DBLE(xninv(8)) + cx*( DBLE(xninv(9))
     +		   + cx*( DBLE(xninv(10)) + zfflo1 )))))
	endif
	if ( xa .gt. bdn01 ) then
	    zfflo1 = cx*( DBLE(xninv(2)) + cx*( DBLE(xninv(3))
     +		   + cx*( DBLE(xninv(4)) + cx*( DBLE(xninv(5))
     +		   + zfflo1 ))))
	endif
	zfflo1 = - cx*( DBLE(xninv(1)) + zfflo1 )
*  #] calculations:
*###] zfflo1:
	end

*###[ zfflo2:
	DOUBLE COMPLEX function zfflo2(x,ier)
***#[*comment:***************************************************
*	calculates log(1-x)+x for |x|<.14 in a faster way to	*
*	~15 significant figures.				*
***#]*comment:***************************************************
*  #[ declarations:
	implicit none
	integer ier
	DOUBLE COMPLEX x,zfflo1,cc
	DOUBLE PRECISION bdn01,bdn05,bdn10,bdn15,bdn18,xprec,xa,
     +		ffbnd,absc
	external ffbnd,zfflo1
	save xprec,bdn01,bdn05,bdn10,bdn15,bdn18
	include 'ff.h'
	absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
*  #] declarations:
*  #[ initialisation:
	data xprec /-1D0/
	if ( xprec .ne. precc ) then
	    xprec = precx
	    precx = precc
*	    determine the boundaries for 1,5,10,15 terms
	    bdn01 = ffbnd(1,1,xninv(2))
	    bdn05 = ffbnd(1,5,xninv(2))
	    bdn10 = ffbnd(1,10,xninv(2))
	    bdn15 = ffbnd(1,15,xninv(2))
	    bdn18 = ffbnd(1,18,xninv(2))
	    precx = xprec
	    xprec = precc
	endif
*  #] initialisation:
*  #[ calculations:
	xa = absc(x)
	if ( xa .gt. bdn18 ) then
	    zfflo2 = zfflo1(x,ier) + x
	    return
	endif
	if ( xa .gt. bdn15 ) then
	    zfflo2 = x*( DBLE(xninv(17)) + x*( DBLE(xninv(18)) +
     +		x*( DBLE(xninv(19)) + x*( DBLE(xninv(20)) ))))
	else
	    zfflo2 = 0
	endif
	if ( xa .gt. bdn10 ) then
	    zfflo2 = x*( DBLE(xninv(12)) + x*( DBLE(xninv(13)) +
     +		x*( DBLE(xninv(14)) + x*( DBLE(xninv(15)) +
     +		x*( DBLE(xninv(16)) + zfflo2 )))))
	endif
	if ( xa .gt. bdn05 ) then
	    zfflo2 = x*( DBLE(xninv(7)) + x*( DBLE(xninv(8)) +
     +		x*( DBLE(xninv(9)) +x*( DBLE(xninv(10)) +
     +		x*( DBLE(xninv(11)) + zfflo2 )))))
	endif
	if ( xa .gt. bdn01 ) then
	    zfflo2 = x*( DBLE(xninv(3)) + x*( DBLE(xninv(4)) +
     +		x*( DBLE(xninv(5)) + x*( DBLE(xninv(6)) + zfflo2 ))))
	endif
	zfflo2 = - x**2*( DBLE(xninv(2)) + zfflo2 )
*  #] calculations:
*###] zfflo2:
	end

*###[ zfflo3:
	DOUBLE COMPLEX function zfflo3(x,ier)
***#[*comment:***************************************************
*	calculates log(1-x)+x+x^2/2 for |x|<.14 in a faster 	*
*	way to ~15 significant figures.				*
***#]*comment:***************************************************
*  #[ declarations:
	implicit none
	integer ier
	DOUBLE COMPLEX x,zfflo2,cc
	DOUBLE PRECISION bdn01,bdn05,bdn10,bdn15,xprec,xa,ffbnd,
     +		absc
	external zfflo2,ffbnd
	save xprec,bdn01,bdn05,bdn10,bdn15
	include 'ff.h'
	absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
*  #] declarations:
*  #[ initialisation:
	data xprec /-1D0/
	if ( xprec .ne. precx ) then
	    xprec = precx
	    precx = precc
*	    determine the boundaries for 1,5,10,15 terms
	    bdn01 = ffbnd(1,1,xninv(3))
	    bdn05 = ffbnd(1,5,xninv(3))
	    bdn10 = ffbnd(1,10,xninv(3))
	    bdn15 = ffbnd(1,15,xninv(3))
	    precx = xprec
	    xprec = precc
	endif
*  #] initialisation:
*  #[ calculations:
	xa = absc(x)
	if ( xa .gt. bdn15 ) then
	    zfflo3 = zfflo2(x,ier) + x**2/2
	    return
	endif
	if ( xa .gt. bdn10 ) then
	    zfflo3 = x*( DBLE(xninv(13)) + x*( DBLE(xninv(14)) +
     +		x*( DBLE(xninv(15)) + x*( DBLE(xninv(16)) +
     +		x*( DBLE(xninv(17)) )))))
	else
	    zfflo3 = 0
	endif
	if ( xa .gt. bdn05 ) then
	    zfflo3 = x*( DBLE(xninv(8)) + x*( DBLE(xninv(9)) +
     +		x*( DBLE(xninv(10)) + x*( DBLE(xninv(11)) +
     +		x*( DBLE(xninv(12)) + zfflo3 )))))
	endif
	if ( xa .gt. bdn01 ) then
	    zfflo3 = x*( DBLE(xninv(4)) + x*( DBLE(xninv(5)) +
     +		x*( DBLE(xninv(6)) + x*( DBLE(xninv(7)) + zfflo3 ))))
	endif
	zfflo3 = - x**3*( DBLE(xninv(3)) + zfflo3 )
*  #] calculations:
*###] zfflo3:
	end


text:
program intrinsics_12
print*,kind(5) ! 4
print*,kind(5_4) ! 4
print*,kind(5_8) ! 8
print*,kind(0.d0) ! 8
print*,kind(0.0) ! 4
print*,kind(5._4) ! 4
print*,kind(5._8) ! 8
print*,kind(.true.) ! 4
end program

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


text:
!===============================================================================
!> Generic DISPATCH main program for Cartesian meshes
!===============================================================================
PROGRAM dispatch
  USE setup_mod
  USE cartesian_mod
  USE dispatcher_mod
  type (cartesian_t):: cartesian
  !.............................................................................
  call setup%init                               ! Standard setup (MPI, I/O, scaling, ...)
  call dispatcher%init                          ! Set dispatcher method
  call cartesian%init                           ! Initialize the patch list and patches
  call dispatcher%execute (cartesian%task_list) ! Cáll dispatcher
  call setup%end
!===============================================================================
END PROGRAM dispatch

max_stars_repo_name: applejwjcat/dispatch

max_issues_repo_name: applejwjcat/dispatch

max_forks_repo_name: applejwjcat/dispatch


text:
	!Subrotina para calcular o Fu, utilizando o ponto de vista lagrangiano
	!Referencia: Casulli (1990, 1992)

	!!! Implementação 15/04/2014
	! Leonardo Romero Monteiro

	!!! Modificações
	! Leonardo Romero Monteiro - 13/01/2015

	!!! Observações
	! no lagrangiano as velocidades no tempo n = t serão calculadas para o tempo n = t-1

	! são calculadas as Fu, Fv e Fw
	SUBROUTINE convdiff()

	USE velpre
	USE parametros
	USE smag
	USE ls_param
	USE vartempo
	USE mms_m
	USE obst
	IMPLICIT NONE

	!===================================================================================================================
	! auxiliares de velocidades: velocidades lagrangianas
	real(8), dimension(nx1,ny,nz) :: uint
	real(8), dimension(nx,ny1,nz) :: vint
	real(8), dimension(nx,ny,nz1) :: wint

	real(8), dimension(0:nx1,ny,nz) :: ama
	real(8), dimension(nx,0:ny1,nz) :: bmb
	real(8), dimension(nx,ny,0:nz1) :: dmd
	real(8), dimension(nx1,ny1,nz)  :: amb, bma
	real(8), dimension(nx1,ny,nz1)  :: amd, dma
	real(8), dimension(nx,ny1,nz1)  :: bmd, dmb

	!real(8), dimension(0:nx1,ny1,nz)  :: vx
	!real(8), dimension(0:nx1,ny,nz1)  :: wx

	!real(8), dimension(nx1,0:ny1,nz)  :: uy
	!real(8), dimension(nx,0:ny1,nz1)  :: wy

	!real(8), dimension(nx1,ny,0:nz1)  :: uz
	!real(8), dimension(nx,ny1,0:nz1)  :: vz

	real(8), dimension(nx1,ny,nz) :: rhox, dhsdx
	real(8), dimension(nx,ny1,nz) :: rhoy, dhsdy
	real(8), dimension(nx,ny,nz1) :: rhoz, dhsdz, epis_z

	!contadores
	integer :: i, j, k
	integer :: ntal
	real(8)    :: tal

	!auxiliares
	real(8) :: aux1, aux2

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

	!reinicializando a rotina
	prd0 = prd1

	!reinicializando a subrotina
	epis_z = 0.

	!===================================================================================================================
	!RESOLUÇÃO DO PROBLEMA
	!===================================================================================================================

	!%%%!-- Advectivo --!%%%!
	if (adv_type == 0) then
		call lagr(uint,vint,wint)
	elseif (adv_type == 1) then
		call classico(uint,vint,wint)
	elseif (adv_type == 2) then
		call rotacional(uint,vint,wint)
	elseif (adv_type == 3) then
		call antissim(uint,vint,wint)
	else
		write(*,*) "sem termo advectivo"
		STOP
	endif

	!Camada esponja. Ativar somente quando necessário.
	!CALL sponge_layer(epis_z)

	!Cálculo dos Fs

	call interpx_cf(rho,nx,ny,nz,rhox) !(nx1,ny,nz)
	call interpy_cf(rho,nx,ny,nz,rhoy) !(nx,ny1,nz)
	call interpz_cf(rho,nx,ny,nz,rhoz) !(nx,ny,nz1)

	call interpx_cf(hsx*kurv/rho,nx,ny,nz,dhsdx) !(nx1,ny,nz)
	call interpy_cf(hsy*kurv/rho,nx,ny,nz,dhsdy) !(nx,ny1,nz)
	call interpz_cf(hsz*kurv/rho,nx,ny,nz,dhsdz) !(nx,ny,nz1)

!	call interpy_cf(u(1:nx1,1:ny,1:nz),nx1,ny,nz,uy(1:nx1,1:ny1,1:nz)) !(nx1,ny1,nz)
!	call interpz_cf(u(1:nx1,1:ny,1:nz),nx1,ny,nz,uz(1:nx1,1:ny,1:nz1)) !(nx1,ny,nz1)

!	call interpx_cf(v(1:nx,1:ny1,1:nz),nx,ny1,nz,vx(1:nx1,1:ny1,1:nz)) !(nx1,ny1,nz)
!	call interpz_cf(v(1:nx,1:ny1,1:nz),nx,ny1,nz,vz(1:nx,1:ny1,1:nz1)) !(nx,ny1,nz1)

!	call interpx_cf(w(1:nx,1:ny,1:nz1),nx,ny,nz1,wx(1:nx1,1:ny,1:nz1)) !(nx1,ny,nz1)
!	call interpy_cf(w(1:nx,1:ny,1:nz1),nx,ny,nz1,wy(1:nx,1:ny1,1:nz1)) !(nx,ny1,nz1)

!      vx(0,:,:)=vx(1,:,:)
!      wx(0,:,:)=wx(1,:,:)

!      uy(:,0,:)=uy(:,1,:)
!      wy(:,0,:)=wy(:,1,:)

!      uz(:,:,0)=uz(:,:,1)
!      vz(:,:,0)=vz(:,:,1)

	call interpy_cf(xnut,nx1,ny,nz,bma) !(nx1,ny1,nz)
	call interpz_cf(xnut,nx1,ny,nz,dma) !(nx1,ny,nz1)

	call interpx_cf(ynut,nx,ny1,nz,amb) !(nx1,ny1,nz)
	call interpz_cf(ynut,nx,ny1,nz,dmb) !(nx,ny1,nz1)

	call interpx_cf(znut,nx,ny,nz1,amd) !(nx1,ny,nz1)
	call interpy_cf(znut,nx,ny,nz1,bmd) !(nx,ny1,nz1)

      ama(1:nx,1:ny,1:nz)=nut(1:nx,1:ny,1:nz)/dx
      bma=bma/dy
      dma=dma/dz

      amb=amb/dx
      bmb(1:nx,1:ny,1:nz)=nut(1:nx,1:ny,1:nz)/dy
      dmb=dmb/dz

      amd=amd/dx
      bmd=bmd/dy
      dmd(1:nx,1:ny,1:nz)=nut(1:nx,1:ny,1:nz)/dz

      ama(0,:,:)=ama(1,:,:)
      bmb(:,0,:)=bmb(:,1,:)
      dmd(:,:,0)=dmd(:,:,1)

      ama(nx1,:,:)=ama(nx,:,:)
      bmb(:,ny1,:)=bmb(:,ny,:)
      dmd(:,:,nz1)=dmd(:,:,nz)

	do k = 1, nz
	do j = 1, ny
	do i = 1, nx1
		Fu(i,j,k) = -uint(i,j,k) + ((ama(i,j,k)*2.*(u(i+1,j,k)-u(i,j,k))-&
                            ama(i-1,j,k)*2.*(u(i,j,k)-u(i-1,j,k)))/dx + &
                            bma(i,j+1,k)*((u(i,j+1,k)-u(i,j,k))/dy+(v(i,j+1,k)-v(i-1,j+1,k))/dx) - &
                            bma(i,j,k)  *((u(i,j,k)-u(i,j-1,k))/dy+(v(i,j,k)  -v(i-1,j,k))/dx  ) + &
                            dma(i,j,k+1)*((u(i,j,k+1)-u(i,j,k))/dz+(w(i,j,k+1)-w(i-1,j,k+1))/dx) - &
                            dma(i,j,k)*  ((u(i,j,k)-u(i,j,k-1))/dz+(w(i,j,k)  -w(i-1,j,k))/dx  ))/rhox(i,j,k) + gx &
			    -sigma*dhsdx(i,j,k) +tf_u(i,j,k) - gz/(chezy*chezy)*sqrt(ub(i,j,k)*ub(i,j,k))*u(i,j,k)
	enddo
	enddo

	do j = 1, ny1
	do i = 1, nx
		Fv(i,j,k) = -vint(i,j,k) + (amb(i+1,j,k)*((v(i+1,j,k)-v(i,j,k))/dx+(u(i+1,j,k)-u(i+1,j-1,k))/dy) - &
                            		    amb(i,j,k)  *((v(i,j,k)-v(i-1,j,k))/dx+(u(i,j,k)  -u(i,j-1,k))/dy  ) + &
                            (bmb(i,j,k)*2.*(v(i,j+1,k)-v(i,j,k)) - &
                            bmb(i,j-1,k)*2.*(v(i,j,k)-v(i,j-1,k)))/dy + &
                            dmb(i,j,k+1)*((v(i,j,k+1)-v(i,j,k))/dz+(w(i,j,k+1)-w(i,j-1,k+1))/dy ) - &
                            dmb(i,j,k)  *((v(i,j,k)-v(i,j,k-1))/dz+(w(i,j,k)  -w(i,j-1,k))/dy   ))/rhoy(i,j,k) &
			    -sigma*dhsdy(i,j,k) +tf_v(i,j,k) - gz/(chezy*chezy)*sqrt(vb(i,j,k)*vb(i,j,k))*v(i,j,k)
	enddo
	enddo
	enddo

	do k = 1, nz1
	do j = 1, ny
	do i = 1, nx
		Fw(i,j,k) = -wint(i,j,k) + (amd(i+1,j,k)*((w(i+1,j,k)-w(i,j,k))/dx+(u(i+1,j,k)-u(i+1,j,k-1))/dz) - &
                            		    amd(i,j,k)  *((w(i,j,k)-w(i-1,j,k))/dx+(u(i,j,k)  -u(i,j,k-1))/dz  ) + &
                            		    bmd(i,j+1,k)*((w(i,j+1,k)-w(i,j,k))/dy+(v(i,j+1,k)-v(i,j+1,k-1))/dz) - &
                            		    bmd(i,j,k)  *((w(i,j,k)-w(i,j-1,k))/dy+(v(i,j,k)  -v(i,j,k-1))/dz  ) + &
                            (dmd(i,j,k)*2.*(w(i,j,k+1)-w(i,j,k))-&
                            dmd(i,j,k-1)*2.*(w(i,j,k)-w(i,j,k-1)))/dz)/rhoz(i,j,k) &
			    - (w(i,j,k)-0.)* epis_z(i,j,k) - gz &
			    -sigma*dhsdz(i,j,k) +tf_w(i,j,k) - gz/(chezy*chezy)*sqrt(wb(i,j,k)*wb(i,j,k))*w(i,j,k)
	enddo
	enddo
	enddo

!==================================================================================================================
END SUBROUTINE convdiff

max_stars_repo_name: pemiguell/SuLi

max_issues_repo_name: pemiguell/SuLi

max_forks_repo_name: pemiguell/SuLi


text:
Subroutine tablem
  Common /table/xarray(0:1000), earray(0:1000)
  Save
  ptmax = 2.01
  anorm = ptdis(ptmax)
  Do l = 0, 200
     x = 0.01*float(l+1)
     rr = ptdis(x)/anorm
     earray(l) = rr
     xarray(l) = x
  End Do
  Return
End Subroutine tablem

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


! Exercise nested function decomposition, gcc/tree-nested.c.

! { dg-do run }

program collapse3
  integer :: p1, p2, p3, p4, p5, p6, p7, p8, p9
  p1 = 2
  p2 = 6
  p3 = -2
  p4 = 4
  p5 = 13
  p6 = 18
  p7 = 1
  p8 = 1
  p9 = 1
  call test1
  call test2 (p1, p2, p3, p4, p5, p6)
  call test3 (p1, p2, p3, p4, p5, p6, p7, p8, p9)
  call test4
contains
  subroutine test1
    integer :: a(3,3,3), k, kk, kkk, l, ll, lll
    !$acc parallel
    !$acc loop collapse(3)
      do 115 k=1,3
dokk:   do kk=1,3
          do kkk=1,3
            a(k,kk,kkk) = 1
          enddo
        enddo dokk
115   continue
    !$acc end parallel
    if (any(a(1:3,1:3,1:3).ne.1)) call abort
    !$acc parallel
    !$acc loop collapse(3)
dol:  do 120 l=1,3
doll:   do ll=1,3
          do lll=1,3
            a(l,ll,lll) = 2
          enddo
        enddo doll
120   end do dol
    !$acc end parallel
    if (any(a(1:3,1:3,1:3).ne.2)) call abort
    end subroutine test1

  subroutine test2(v1, v2, v3, v4, v5, v6)
    integer :: i, j, k, a(1:7, -3:5, 12:19), b(1:7, -3:5, 12:19)
    integer :: v1, v2, v3, v4, v5, v6
    logical :: l, r
    l = .false.
    r = .false.
    a(:, :, :) = 0
    b(:, :, :) = 0
    !$acc parallel reduction (.or.:l)
    !$acc loop reduction (.or.:l) collapse (3)
      do i = v1, v2
        do j = v3, v4
          do k = v5, v6
            l = l.or.i.lt.2.or.i.gt.6.or.j.lt.-2.or.j.gt.4
            l = l.or.k.lt.13.or.k.gt.18
            if (.not.l) a(i, j, k) = a(i, j, k) + 1
          end do
        end do
      end do
    !$acc end parallel
    do i = v1, v2
      do j = v3, v4
        do k = v5, v6
          r = r.or.i.lt.2.or.i.gt.6.or.j.lt.-2.or.j.gt.4
          r = r.or.k.lt.13.or.k.gt.18
          if (.not.l) b(i, j, k) = b(i, j, k) + 1
        end do
      end do
    end do
    if (l .neqv. r) call abort
    do i = v1, v2
      do j = v3, v4
        do k = v5, v6
           if (a(i, j, k) .ne. b(i, j, k)) call abort
        end do
      end do
    end do
  end subroutine test2

  subroutine test3(v1, v2, v3, v4, v5, v6, v7, v8, v9)
    integer :: i, j, k, a(1:7, -3:5, 12:19), b(1:7, -3:5, 12:19)
    integer :: v1, v2, v3, v4, v5, v6, v7, v8, v9
    logical :: l, r
    l = .false.
    r = .false.
    a(:, :, :) = 0
    b(:, :, :) = 0
    !$acc parallel reduction (.or.:l)
    !$acc loop reduction (.or.:l) collapse (3)
      do i = v1, v2, v7
        do j = v3, v4, v8
          do k = v5, v6, v9
            l = l.or.i.lt.2.or.i.gt.6.or.j.lt.-2.or.j.gt.4
            l = l.or.k.lt.13.or.k.gt.18
            if (.not.l) a(i, j, k) = a(i, j, k) + 1
          end do
        end do
      end do
    !$acc end parallel
    do i = v1, v2, v7
      do j = v3, v4, v8
        do k = v5, v6, v9
          r = r.or.i.lt.2.or.i.gt.6.or.j.lt.-2.or.j.gt.4
          r = r.or.k.lt.13.or.k.gt.18
          if (.not.l) b(i, j, k) = b(i, j, k) + 1
        end do
      end do
    end do
    if (l .neqv. r) call abort
    do i = v1, v2, v7
      do j = v3, v4, v8
        do k = v5, v6, v9
           if (a(i, j, k) .ne. b(i, j, k)) call abort
        end do
      end do
    end do
  end subroutine test3

  subroutine test4
    integer :: i, j, k, a(1:7, -3:5, 12:19), b(1:7, -3:5, 12:19)
    integer :: v1, v2, v3, v4, v5, v6, v7, v8, v9
    logical :: l, r
    l = .false.
    r = .false.
    a(:, :, :) = 0
    b(:, :, :) = 0
    v1 = p1
    v2 = p2
    v3 = p3
    v4 = p4
    v5 = p5
    v6 = p6
    v7 = p7
    v8 = p8
    v9 = p9
    !$acc parallel reduction (.or.:l)
    !$acc loop reduction (.or.:l) collapse (3)
      do i = v1, v2, v7
        do j = v3, v4, v8
          do k = v5, v6, v9
            l = l.or.i.lt.2.or.i.gt.6.or.j.lt.-2.or.j.gt.4
            l = l.or.k.lt.13.or.k.gt.18
            if (.not.l) a(i, j, k) = a(i, j, k) + 1
          end do
        end do
      end do
    !$acc end parallel
    do i = v1, v2, v7
      do j = v3, v4, v8
        do k = v5, v6, v9
          r = r.or.i.lt.2.or.i.gt.6.or.j.lt.-2.or.j.gt.4
          r = r.or.k.lt.13.or.k.gt.18
          if (.not.r) b(i, j, k) = b(i, j, k) + 1
        end do
      end do
    end do
    if (l .neqv. r) call abort
    do i = v1, v2, v7
      do j = v3, v4, v8
         do k = v5, v6, v9
           if (a(i, j, k) .ne. b(i, j, k)) call abort
         end do
      end do
    end do
  end subroutine test4

end program collapse3


text:
subroutine print_help
  WRITE (*, *) "Required Options: --asymmetry path_to_file --ssa path_to_file --solar-zenith angle_in_degrees"
  CALL EXIT(-1)
end subroutine print_help

program disort_albedo
  use f90getopt
  use disort_variables
  implicit none

  ! CLI values
  CHARACTER(len=100) :: ASYM_FILE, SSA_FILE
  INTEGER            :: IO_STATUS

  ! Maximum number of bands to simulate
  INTEGER :: BANDS = 0
  ! Effective radius in micro-meter
  INTEGER :: GRAINSIZE = 148
  ! Solar zenith angle
  REAL :: ANGLE = 0.0

  ! ANGLE conversion
  REAL(kind=4),PARAMETER :: PI = 2.*ASIN(1.0)
  REAL, PARAMETER :: Degree180 = 180.0
  ! Degrees to radians
  REAL            :: D_to_R    = PI/Degree180

  ! File data inputs
  REAL, DIMENSION(:,:), ALLOCATABLE :: ASYMF
  REAL, DIMENSION(:,:), ALLOCATABLE :: SALB
  ! Results
  REAL, DIMENSION(:,:), ALLOCATABLE :: ALBS

  INTEGER :: IU, I, J
  CHARACTER  HEADER*127

  ! Command line options
  type(option_s):: opts(4)
  opts(1) = option_s( "asymmetry",  .TRUE.,  'a' )
  opts(2) = option_s( "ssa",  .TRUE., 's')
  opts(3) = option_s( "solar-zenith",  .TRUE., 'z')
  opts(4) = option_s( "help", .false.,  'h')

  if (command_argument_count() /= 6 ) call print_help

  ! Process command line options one by one
  DO
    SELECT CASE( getopt( "ash:", opts ) )
      CASE( char(0) )
        EXIT
      CASE( 'a' )
        ASYM_FILE = optarg
      CASE( 's' )
        SSA_FILE = optarg
      CASE( 'z' )
        READ(optarg, '(f10.0)') ANGLE
        ANGLE = COS(ANGLE * D_to_R)
        print *, ANGLE
      CASE( 'h' )
        call print_help
    END SELECT
  END DO

  ! Get number of lines, equal to bands, from ASYM_FILE
  OPEN(10, file=ASYM_FILE, iostat=IO_STATUS, status='old')
  IF (IO_STATUS/=0) STOP 'Cannot open assymetry file'

  DO
    READ(10, *, iostat=IO_STATUS)
    IF (IO_STATUS/=0) EXIT
    BANDS = BANDS + 1
  END DO
  CLOSE(10)

  allocate(ASYMF(GRAINSIZE, BANDS))
  allocate(SALB(GRAINSIZE, BANDS))
  allocate(ALBS(GRAINSIZE, BANDS))

  OPEN(unit=8,file=ASYM_FILE)
  READ(8,*) ASYMF
  CLOSE(8)

  OPEN(unit=10,file=SSA_FILE, iostat=IO_STATUS, status='old')
  IF (IO_STATUS/=0) STOP 'Cannot open ssa file'
  READ(10,*) SALB
  CLOSE(10)

  ACCUR = 0.0

  DO I = 1, BANDS
    DO J = 1, GRAINSIZE
      USRTAU    = .FALSE.
      USRANG    = .TRUE.
      LAMBER    = .FALSE.
      PLANK     = .FALSE.
      DO_PSEUDO_SPHERE = .FALSE.
      DELTAMPLUS = .FALSE.

      NSTR = 16
      IF( MOD(NSTR,2) /= 0 ) NSTR = NSTR + 1
      NLYR = 1
      NMOM = NSTR
      NTAU = 4
      IF( .NOT. USRTAU ) NTAU = NLYR + 1
      !NUMU      = 1; IF(.NOT.USRANG) NUMU = NSTR
      NPHI = 1

      IBCND  = 1
      NLYR   = 1
      ONLYFL = .FALSE.
      USRANG = .TRUE.
      NUMU   = 1
      IF( .NOT. USRANG ) NUMU = NSTR

      IF( USRANG .AND. IBCND == 1 ) THEN
        NUMU_O = NUMU
        NUMU = 2*NUMU
      END IF

      call allocate_disort_allocatable_arrays( NLYR, NMOM, NSTR, NUMU, NPHI, NTAU )

      DTAUC( 1 ) = 1000.0
      SSALB( 1 ) = SALB(J,I)
      CALL  GETMOM( 3, ASYMF(J,I), NMOM, PMOM )
      PRNT( 4 )  = .FALSE.
      PRNT( 2 )  = .FALSE.
      PRNT( 1 )  = .FALSE.
      ! Solar zenith angle
      UMU( 1 )   =  ANGLE

      DO IU = 1, NUMU_O
        UMU( IU + NUMU_O ) = UMU( IU )
      ENDDO

      DO IU = 1, NUMU_O
        UMU( IU ) = -UMU( 2*NUMU_O + 1 - IU )
      ENDDO

      ALBEDO = 0.0
      HEADER = ''

      CALL DISORT( NLYR, NMOM, NSTR, NUMU, NPHI, NTAU,           &
               USRANG, USRTAU, IBCND, ONLYFL, PRNT,          &
               PLANK, LAMBER, DELTAMPLUS, DO_PSEUDO_SPHERE,  &
               DTAUC, SSALB, PMOM, TEMPER, WVNMLO, WVNMHI,   &
               UTAU, UMU0, PHI0, UMU, PHI, FBEAM,            &
               FISOT, ALBEDO, BTEMP, TTEMP, TEMIS,           &
               EARTH_RADIUS, H_LYR,                          &
               RHOQ, RHOU, RHO_ACCURATE, BEMST, EMUST,       &
               ACCUR,  HEADER,                               &
               RFLDIR, RFLDN, FLUP, DFDT, UAVG, UU,          &
               ALBMED, TRNMED )

      ALBS(J,I) = ALBMED(1)

      call deallocate_disort_allocatable_arrays()
    ENDDO
  ENDDO

  deallocate(ASYMF)
  deallocate(SALB)

  ! output file with albedo
  OPEN(unit = 1,file = "albedo.csv")
  DO I = 1, BANDS
    WRITE(1,*) (ALBS(J,I), J=1, GRAINSIZE)
  ENDDO
  CLOSE(1)

  deallocate(ALBS)

end program disort_albedo

max_stars_repo_name: UofU-Cryosphere/spectro-dp

max_issues_repo_name: UofU-Cryosphere/spectro-dp

max_forks_repo_name: UofU-Cryosphere/spectro-dp