text:
c###sigdis.for
      SUBROUTINE SIGDIS
c***********************************************
c          Alex Shovkoplyas made some changes to:
c             curmuf
c             seltxr
c             settxr
c             sigdis
c          Normally if the version (..\database\version.w32)
c             is yy.mmddW then the original version is used.
c                yy.mmddA then Alex's version is used.
c***********************************************
      use cversn
c jw      common /CVERSN/ VERSN
c jw      character VERSN*8
      if(versn(8:8).eq.'A' .or. versn(8:8).eq.'a') then
         call sigdis_alex       !  use Alex's modified version
      else
         call sigdis_orig       !  use original version
      end if
      return
      end
C--------------------------------
      SUBROUTINE SIGDIS_orig
C--------------------------------
C
C  ADJUST SIGNAL DISTRIBUTION TABLES FOR THIS PATH AND SET ABSORPTION
C  LOSS PARAMETERS
C
      use cversn
c jw      common /cversn/ versn
c jw         character versn*8      !  look for "H" to ignore absorption fix
      COMMON / FILES / LUI, LUO, LU2, LU5, LU6, LU15, LU16, LU20, LU25,
     A LU26, LU35
      COMMON/CON/D2R,DCL,GAMA,PI,PI2,PIO2,R2D,RZ,VOFL
      COMMON/FRQ/FREA(13),FREL(29),FREQ,JMODE,ITXRCP(2)
      COMMON/GEOG/GYZ(5),RAT(5),GMDIP(5),CLCK(5),ABIY(5),ARTIC(5),SIGPAT
     A(5), EPSPAT(5)
      COMMON/MUFS/EMUF(24),F1MUF(24),F2MUF(24),ESMUF(24),ALLMUF(24),FOT
     A (24),XLUF(24),HPF(24),ANGMUF(24),MODMUF,SIGL(4),SIGU(4),DELMUF(4)
     B ,HPMUF(4),HTMUF(4),FVMUF(4),AFMUF(4),NHOPMF(4),YFOT(4),YHPF(4)
     C ,YMUF(4)
      COMMON /RON /CLAT(5), CLONG(5), GLAT(5), RD(5), FI(3,5), YI(3,5),
     1HI(3,5), HPRIM(30,5), HTRUE(30,5), FVERT(30,5),KM,KFX, AFAC(30,5),
     2HTR(50,3), FNSQ(50,3)
      COMMON/SIGD/DSL,ASM,DSU,AGLAT,DSLF,ASMF,DSUF,ACAV,FEAV,AFE,BFE,HNU
     A ,HTLOSS,XNUZ,XVE
      COMMON/TON/ADJ,ADS,GNOS,GOT,REL,SL,SLS,SPR,SU,SUS
     A ,XNOISE,ZNOISE,NF
      COMMON/INFORM/INFO,IHSHR,IHLNG
      COMMON /DON /ALATD, AMIN, AMIND, BTR, BTRD, DLONG, DMP, ERTR, GCD,
     1 GCDKM, PMP, PWR, TLAT, TLATD, TLONG, TLONGD, RSN, SIGTR, RLAT,
     2 RLATD,RLONG,RLONGD,BRTD,FLUX,ULAT,ULATD,ULONG,ULONGD,SSN,D90R,
     3 D50R,D10R,D90S,D50S,D10S
      GLAV = 0.0
C.....SET NUMBER OF IONOSPHERES
      KFXX=KFX
ccc      KFXX=KM         !  temp change by Greg for receprocity try

      XKF  = KFXX
      XADJ = 0.0
      XSU  = 0.0
      XSL  = 0.0
      XFMP = 0.0
      XSUP = 0.0
      XSLP = 0.0
      AC = 0.0
      FEAV=0.0
ccc      write(luo,'(''in sigdis, kfx='',2i5)') kfx,km
      DO 115 K = 1,KFXX
C.....ABSORPTION INDEX
      ABIY(K) =  - 0.04 + EXP(-2.937 + 0.8445*FI(1,K) )
c****** Change made 6/2/99 by Greg Hand to fix absorption error
c       that was never fixed in VOACAP.
c******  per change to ICEPAC 3/18/93 by Frank Stewart to fix losses problem
c          this limits ABIY to a lower limit of .1 because that is the
c          lower limit of the curve.
c          IONCAP was wrong when it did not have the lower limit.
c          It was too optimistic with its prediction.
ccc      write(12,'(''abiy('',i1,'')='',f8.4)') k,abiy(k)
c          the "HAM" version does not do the absorption fix
c          to create the "HAM" version, edit the file:
c                  ..\database\version.w32
c             change the "W" to "H".
c          this was changed 2/8/2006 for Dean Straw
c          this was changed 4/18/2006 to "I" from "H". "I"=IONCAP
c          "I" or "a" is IONCAP. "a" is Alex's changes
      if(versn(8:8).eq.'I' .or. versn(8:8).eq.'a') then
      else
         ABIY(K)=MAX(0.1,ABIY(K))
      end if
C.....SYSTEM LOSS TABLE
      DUMMY1 = GLAT(K)
      DUMMY2 = CLCK(K)
      IDP=2
      IF(GCDKM.GT.2500.)IDP=5
      CALL SYSSY(DUMMY1,DUMMY2,IDP,ADJ,SU,SL,FMP,SUP,SLP)
ccc      write(luo,122) k,dummy1,dummy2,idp,adj,su,sl,fmp,sup,slp
ccc     +                ,abiy(k),fi(1,k)
ccc122   format('syssy=',i5,2e15.7,i5,8f10.4)
C.....MEDIAN
      XADJ = XADJ + ADJ
C.....UPPER DECILE
      XSU  = XSU  + SU
C.....LOWER DECILE
      XSL  = XSL  + SL
C.....PREDICTION ERROR, MEDIAN
      XFMP = XFMP + FMP
C.....UPPER DECILE
      XSUP = XSUP + SUP
C.....LOWER DECILE
      XSLP = XSLP + SLP
C.....E CRITICAL
      FEAV = FEAV + FI(1,K)
C.....ABSORBTION INDEX
      AC   = AC   + ABIY(K)
      ARTIC(K) = ADJ
C.....GEOMAGNETIC LATITUDE
      GLAV = GLAV + ABS(GLAT(K))
  115 CONTINUE
C.....BEGIN SECTION TO SAVE THE AVERAGE VALUES
      GLAV= GLAV/XKF
      AGLAT = GLAV
      ACAV= AC/XKF
      FEAV= FEAV/XKF
      ADJ = XADJ/XKF
ccc      write(luo,'(''adj='',4f11.4)') adj,xadj,xkf,acav
      SU  = XSU/XKF
      SL  = XSL/XKF
      ADS = XFMP/XKF
      SUS = XSUP/XKF
      SLS = XSLP/XKF
C.....END SECTION TO SAVE THE AVERAGE VALUES
      K = JMODE
C
C  PARAMETER FOR D-E ABSORTION LOSS.
C  ABIY( )  IS ABSORPTION INDEX. ACAV IS THE AVERAGE FOR THE PATH
C
C.....D - E REGION LOSS ADJUSTMENT FACTOR
      HTLOSS = 88.
      XVE    = XLIN( 90., HTRUE,FVERT,30, K )
      XVE    = XVE/FI(1,K)
      XNUZ = 63.07
      HNU  = 4.39
      IF(FEAV - 2.0) 125,125,120
C.....ADJUSTMENT TO CCIR 252 (HAYDON,LUCAS) LOSS EQUATION FOR E MODES
  120 AFE  = 1.359
      BFE  = 8.617
      GO TO 140
  125 IF(FEAV - 0. 5) 130,130,135
  130 AFE  = 0.0
      BFE  = 0.0
      GO TO 140
  135 AFE  = 1.359 *(FEAV - 0.5)/1.5
      BFE  = 8.617 *(FEAV - 0.5)/1.5
  140 CONTINUE
C
C  SET SIGNAL DISTRIBUTION TABLES. TABLE IS FOR FREQ = FTAB
C
C.....USE FOT, F2 LAYER
      GLAV = GLAV*R2D
      IF( GLAV - 40.)  145,145,150
  145 FTAB = YFOT(3)
      GO TO 165
  150 IF(GLAV  - 50.)  155,155,160
  155 FTAB = YFOT(3)
C.....INTERPOLATE IN BETWEEN
      FTAB = FTAB  - (GLAV - 40.)*(FTAB - 10.)/ 10.
      GO TO 165
C.....SET TO 10 MHZ (NEAR POLES)
  160 FTAB = 10.
  165 CONTINUE
      ESLSM = 0.0
      IF(YMUF(4).gt.0.) then
C   ES CONTRIBUTION TO TABLE( OBSCURTION LOSS).
C
C.....NOW REMOVE ES OBSCURATION AND F2 OVER-THE-MUF AT FTAB FROM TABLES
C.....AND WILL REPLACE FOR EACH MODE AND FREQUENCY AS NECESSARY
         DUMMY = YMUF(4)
         PES = PRBMUF(FTAB,DUMMY,DUMMY,4)
         PES =AMAX1(0.1,PES)
         PES =AMIN1(0.9,PES)
         ESLSM = -10.*ALOG10(1.-PES)
      end if
C  F2 OVER-THE-MUF CONTRIBUTION TO TABLE
C.....PROBABILITY CALCULATION
      DUMMY = YMUF(3)
      PF2 = PRBMUF(FTAB,DUMMY,DUMMY,3)
      PF2 = AMAX1(0.1, PF2)
ccc      PF2 = AMIN1(0.9, PF2)       !  this added 5/12/2005 GRH
      F2LSM = -10.*ALOG10(PF2)
      IF(IAND(INFO,2).GT.0)THEN
        WRITE(99,'(A,A,3F6.2)')' "SIGDIS" OTM(F2 NO SEC) LOSS E',
     1 'XCL FROM TAB 6 FTAB,YMUF,LOS=',FTAB,YMUF(3),F2LSM
      ENDIF
C RESIDUAL (AURORAL) LOSS ADJUSTMENT TO MEDIAN SIGNAL LEVEL.
      ASM = ADJ  - ESLSM - F2LSM
ccc      write(luo,123) asm,adj,eslsm,f2lsm,pf2
ccc123   format('ASM=',5f11.4)
      ASM = AMAX1(ASM, 0.0 )
C  UPPER DECILE
      PES = 0.0
      IF(YFOT(4).gt.0.) then
         DUMMY = YFOT(4)
         PES = PRBMUF(FTAB,DUMMY,DUMMY,4)
         PES = AMAX1( 0.1,PES)
         PES = AMIN1( 0.9,PES)
      end if
      DUMMY = YHPF(3)
      PF2 = PRBMUF(FTAB,DUMMY,DUMMY,3)
      PF2 = AMAX1(0.1, PF2)
ccc      PF2 = AMIN1(0.9, PF2)       !  this added 5/12/2005 GRH
C  UPPER DECILE SIGNAL LEVEL ADJUSTMENT TO MEDIAN
      DSU = 1.28* SL  -( +10.*ALOG10(1.- PES) + ESLSM)
     A                -( +10.*ALOG10(PF2) + F2LSM)
      DSU = AMAX1 ( DSU, 0.5)
C  LOWER DECILE.
      PES = 0.0
      IF(YHPF(4).gt.0.) then
         DUMMY = YHPF(4)
         PES = PRBMUF(FTAB,DUMMY,DUMMY,4)
         PES= AMAX1(0.1,PES)
         PES= AMIN1(0.9, PES)
      end if
      DUMMY = YFOT(3)
      PF2 = PRBMUF(FTAB,DUMMY,DUMMY,3)
      PF2 = AMAX1( 0.1,PF2)
ccc      PF2 = AMIN1(0.9, PF2)       !  this added 5/12/2005 GRH
C LOWER DECILE SIGNAL LEVEL ADJUSTMENT TO MEDIAN
      DSL = 1.28* SU - (-10.* ALOG10(1.-PES) - ESLSM)
     A               - (-10.* ALOG10(PF2) - F2LSM)
      DSL = AMAX1(1.,DSL)
      RETURN
      END
C--------------------------------
c###sigdis.for
      SUBROUTINE SIGDIS_alex
C--------------------------------
C
C  ADJUST SIGNAL DISTRIBUTION TABLES FOR THIS PATH AND SET ABSORPTION
C  LOSS PARAMETERS
C
      use cversn
c jw      common /cversn/ versn
c jw         character versn*8      !  look for "H" to ignore absorption fix
      COMMON / FILES / LUI, LUO, LU2, LU5, LU6, LU15, LU16, LU20, LU25,
     A LU26, LU35
      COMMON/CON/D2R,DCL,GAMA,PI,PI2,PIO2,R2D,RZ,VOFL
      COMMON/FRQ/FREA(13),FREL(29),FREQ,JMODE,ITXRCP(2)
      COMMON/GEOG/GYZ(5),RAT(5),GMDIP(5),CLCK(5),ABIY(5),ARTIC(5),SIGPAT
     A(5), EPSPAT(5)
      COMMON/MUFS/EMUF(24),F1MUF(24),F2MUF(24),ESMUF(24),ALLMUF(24),FOT
     A (24),XLUF(24),HPF(24),ANGMUF(24),MODMUF,SIGL(4),SIGU(4),DELMUF(4)
     B ,HPMUF(4),HTMUF(4),FVMUF(4),AFMUF(4),NHOPMF(4),YFOT(4),YHPF(4)
     C ,YMUF(4)
      COMMON /RON /CLAT(5), CLONG(5), GLAT(5), RD(5), FI(3,5), YI(3,5),
     1HI(3,5), HPRIM(30,5), HTRUE(30,5), FVERT(30,5),KM,KFX, AFAC(30,5),
     2HTR(50,3), FNSQ(50,3)
      COMMON/SIGD/DSL,ASM,DSU,AGLAT,DSLF,ASMF,DSUF,ACAV,FEAV,AFE,BFE,HNU
     A ,HTLOSS,XNUZ,XVE
      COMMON/TON/ADJ,ADS,GNOS,GOT,REL,SL,SLS,SPR,SU,SUS
     A ,XNOISE,ZNOISE,NF
      COMMON/INFORM/INFO,IHSHR,IHLNG
      COMMON /DON /ALATD, AMIN, AMIND, BTR, BTRD, DLONG, DMP, ERTR, GCD,
     1 GCDKM, PMP, PWR, TLAT, TLATD, TLONG, TLONGD, RSN, SIGTR, RLAT,
     2 RLATD,RLONG,RLONGD,BRTD,FLUX,ULAT,ULATD,ULONG,ULONGD,SSN,D90R,
     3 D50R,D10R,D90S,D50S,D10S
c>>>>>A.S.
      DIMENSION LX(3)
      DATA LX/1,3,5/
c<<<<<
      GLAV = 0.0
C.....SET NUMBER OF IONOSPHERES
      KFXX=KFX
ccc      KFXX=KM         !  temp change by Greg for receprocity try

      XKF  = KFXX
      XADJ = 0.0
      XSU  = 0.0
      XSL  = 0.0
      XFMP = 0.0
      XSUP = 0.0
      XSLP = 0.0
      AC = 0.0
      FEAV=0.0
ccc      write(luo,'(''in sigdis, kfx='',2i5)') kfx,km
      DO 115 K = 1,KFXX
C.....ABSORPTION INDEX
      ABIY(K) =  - 0.04 + EXP(-2.937 + 0.8445*FI(1,K) )
c****** Change made 6/2/99 by Greg Hand to fix absorption error
c       that was never fixed in VOACAP.
c******  per change to ICEPAC 3/18/93 by Frank Stewart to fix losses problem
c          this limits ABIY to a lower limit of .1 because that is the
c          lower limit of the curve.
c          IONCAP was wrong when it did not have the lower limit.
c          It was too optimistic with its prediction.
ccc      write(12,'(''abiy('',i1,'')='',f8.4)') k,abiy(k)
c          the "HAM" version does not do the absorption fix
c          to create the "HAM" version, edit the file:
c                  ..\database\version.w32
c             change the "W" to "H".
c          this was changed 2/8/2006 for Dean Straw
c          this was changed 4/18/2006 to "I" from "H". "I"=IONCAP
c          "I" or "a" is IONCAP. "a" is Alex's changes
      if(versn(8:8).eq.'I' .or. versn(8:8).eq.'a') then
      else
         ABIY(K)=MAX(0.1,ABIY(K))
      end if
C.....SYSTEM LOSS TABLE
      DUMMY1 = GLAT(LX(K))
      DUMMY2 = CLCK(LX(K))
      IDP=2
      IF(GCDKM.GT.2500.)IDP=5
      CALL SYSSY(DUMMY1,DUMMY2,IDP,ADJ,SU,SL,FMP,SUP,SLP)
ccc      write(luo,122) k,dummy1,dummy2,idp,adj,su,sl,fmp,sup,slp
ccc     +                ,abiy(k),fi(1,k)
ccc122   format('syssy=',i5,2e15.7,i5,8f10.4)
C.....MEDIAN
      XADJ = XADJ + ADJ
C.....UPPER DECILE
      XSU  = XSU  + SU
C.....LOWER DECILE
      XSL  = XSL  + SL
C.....PREDICTION ERROR, MEDIAN
      XFMP = XFMP + FMP
C.....UPPER DECILE
      XSUP = XSUP + SUP
C.....LOWER DECILE
      XSLP = XSLP + SLP
C.....E CRITICAL
      FEAV = FEAV + FI(1,K)
C.....ABSORBTION INDEX
      AC   = AC   + ABIY(K)
      ARTIC(K) = ADJ
C.....GEOMAGNETIC LATITUDE
      GLAV = GLAV + ABS(GLAT(LX(K)))
  115 CONTINUE
C.....BEGIN SECTION TO SAVE THE AVERAGE VALUES
      GLAV= GLAV/XKF
      AGLAT = GLAV
      ACAV= AC/XKF
      FEAV= FEAV/XKF
      ADJ = XADJ/XKF
ccc      write(luo,'(''adj='',4f11.4)') adj,xadj,xkf,acav
      SU  = XSU/XKF
      SL  = XSL/XKF
      ADS = XFMP/XKF
      SUS = XSUP/XKF
      SLS = XSLP/XKF
C.....END SECTION TO SAVE THE AVERAGE VALUES
      K = JMODE
C
C  PARAMETER FOR D-E ABSORTION LOSS.
C  ABIY( )  IS ABSORPTION INDEX. ACAV IS THE AVERAGE FOR THE PATH
C
C.....D - E REGION LOSS ADJUSTMENT FACTOR
      HTLOSS = 88.
      XVE    = XLIN( 90., HTRUE,FVERT,30, K )
      XVE    = XVE/FI(1,K)
      XNUZ = 63.07
      HNU  = 4.39
      IF(FEAV - 2.0) 125,125,120
C.....ADJUSTMENT TO CCIR 252 (HAYDON,LUCAS) LOSS EQUATION FOR E MODES
  120 AFE  = 1.359
      BFE  = 8.617
      GO TO 140
  125 IF(FEAV - 0. 5) 130,130,135
  130 AFE  = 0.0
      BFE  = 0.0
      GO TO 140
  135 AFE  = 1.359 *(FEAV - 0.5)/1.5
      BFE  = 8.617 *(FEAV - 0.5)/1.5
  140 CONTINUE
C
C  SET SIGNAL DISTRIBUTION TABLES. TABLE IS FOR FREQ = FTAB
C
C.....USE FOT, F2 LAYER
      GLAV = GLAV*R2D
      IF( GLAV - 40.)  145,145,150
  145 FTAB = YFOT(3)
      GO TO 165
  150 IF(GLAV  - 50.)  155,155,160
  155 FTAB = YFOT(3)
C.....INTERPOLATE IN BETWEEN
      FTAB = FTAB  - (GLAV - 40.)*(FTAB - 10.)/ 10.
      GO TO 165
C.....SET TO 10 MHZ (NEAR POLES)
  160 FTAB = 10.
  165 CONTINUE
      ESLSM = 0.0
      IF(YMUF(4).gt.0.) then
C   ES CONTRIBUTION TO TABLE( OBSCURTION LOSS).
C
C.....NOW REMOVE ES OBSCURATION AND F2 OVER-THE-MUF AT FTAB FROM TABLES
C.....AND WILL REPLACE FOR EACH MODE AND FREQUENCY AS NECESSARY
         DUMMY = YMUF(4)
         PES = PRBMUF(FTAB,DUMMY,DUMMY,4)
         PES =AMAX1(0.1,PES)
         PES =AMIN1(0.9,PES)
         ESLSM = -10.*ALOG10(1.-PES)
      end if
C  F2 OVER-THE-MUF CONTRIBUTION TO TABLE
C.....PROBABILITY CALCULATION
      DUMMY = YMUF(3)
      PF2 = PRBMUF(FTAB,DUMMY,DUMMY,3)
      PF2 = AMAX1(0.1, PF2)
ccc      PF2 = AMIN1(0.9, PF2)       !  this added 5/12/2005 GRH
      F2LSM = -10.*ALOG10(PF2)
      IF(IAND(INFO,2).GT.0)THEN
        WRITE(99,'(A,A,3F6.2)')' "SIGDIS" OTM(F2 NO SEC) LOSS E',
     1 'XCL FROM TAB 6 FTAB,YMUF,LOS=',FTAB,YMUF(3),F2LSM
      ENDIF
C RESIDUAL (AURORAL) LOSS ADJUSTMENT TO MEDIAN SIGNAL LEVEL.
      ASM = ADJ  - ESLSM - F2LSM
ccc      write(luo,123) asm,adj,eslsm,f2lsm,pf2
ccc123   format('ASM=',5f11.4)
      ASM = AMAX1(ASM, 0.0 )
C  UPPER DECILE
      PES = 0.0
      IF(YFOT(4).gt.0.) then
         DUMMY = YFOT(4)
         PES = PRBMUF(FTAB,DUMMY,DUMMY,4)
         PES = AMAX1( 0.1,PES)
         PES = AMIN1( 0.9,PES)
      end if
      DUMMY = YHPF(3)
      PF2 = PRBMUF(FTAB,DUMMY,DUMMY,3)
      PF2 = AMAX1(0.1, PF2)
ccc      PF2 = AMIN1(0.9, PF2)       !  this added 5/12/2005 GRH
C  UPPER DECILE SIGNAL LEVEL ADJUSTMENT TO MEDIAN
      DSU = 1.28* SL  -( +10.*ALOG10(1.- PES) + ESLSM)
     A                -( +10.*ALOG10(PF2) + F2LSM)
      DSU = AMAX1 ( DSU, 0.5)
C  LOWER DECILE.
      PES = 0.0
      IF(YHPF(4).gt.0.) then
         DUMMY = YHPF(4)
         PES = PRBMUF(FTAB,DUMMY,DUMMY,4)
         PES= AMAX1(0.1,PES)
         PES= AMIN1(0.9, PES)
      end if
      DUMMY = YFOT(3)
      PF2 = PRBMUF(FTAB,DUMMY,DUMMY,3)
      PF2 = AMAX1( 0.1,PF2)
ccc      PF2 = AMIN1(0.9, PF2)       !  this added 5/12/2005 GRH
C LOWER DECILE SIGNAL LEVEL ADJUSTMENT TO MEDIAN
      DSL = 1.28* SU - (-10.* ALOG10(1.-PES) - ESLSM)
     A               - (-10.* ALOG10(PF2) - F2LSM)
      DSL = AMAX1(1.,DSL)
      RETURN
      END
C--------------------------------

max_stars_repo_name: jawatson/voacapl

max_issues_repo_name: jawatson/voacapl

max_forks_repo_name: jawatson/voacapl


text:
LINK Tutoring Services offers K12 inhome tutoring in Sacramento and Yolo County Yolo counties.  They connect tutors and students in need of assistance. They claim to have the best prices and to be able to find tutors for virtually any subject at any level.

 Subjects
  Math
  Science
  English
  Reading
  Spanish
  French
  Test preparation

   Need help with classwork?  Davis Wiki has a list of tutoring tutoring resources!

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


           program demo_system_putenv
           use M_system, only : system_putenv
           use iso_c_binding
           implicit none
           integer :: ierr
              !
              write(*,'(a)')'no environment variables containing "GRU":'
              call execute_command_line('env|grep GRU')
              !
              call system_putenv('GRU=this is the value',ierr)
              write(*,'(a,i0)')'now "GRU" should be defined: ',ierr
              call execute_command_line('env|grep GRU')
              !
              call system_putenv('GRU2=this is the second value',ierr)
              write(*,'(a,i0)')'now "GRU" and "GRU2" should be defined: ',ierr
              call execute_command_line('env|grep GRU')
              !
              call system_putenv('GRU2',ierr)
              call system_putenv('GRU',ierr)
              write(*,'(a,i0)')'should be gone, varies with different putenv(3c): ',ierr
              call execute_command_line('env|grep GRU')
              write(*,'(a)')'system_unsetenv(3f) is a better way to remove variables'
              !
           end program demo_system_putenv


text:
      SUBROUTINE IFP4G (IBIT,FILE)
C
C     TURNS ON BIT -IBIT- IN TRAILER FOR DATA BLOCK -FILE-
C
      EXTERNAL    ORF
      INTEGER     ORF, TRAIL(7), FILE, TWO
      COMMON/TWO/ TWO(32)
C
      TRAIL(1) = FILE
      CALL RDTRL (TRAIL)
      I1 = (IBIT-1)/16 + 2
      I2 = IBIT - (I1-2)*16 + 16
      TRAIL(I1) = ORF(TRAIL(I1),TWO(I2))
      TRAIL(1) = FILE
      CALL WRTTRL (TRAIL)
      RETURN
      END


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine initialcfd(nef,ipkonf,kon,lakonf,co,coel,cofa,nface,
     &  ielfa,area,ipnei,neiel,xxn,xxi,xle,xlen,xlet,xrlfa,cosa,
     &  volume,neifa,xxj,cosb,vel,dmin)
!
!     calculating geometric variables of the cells and their faces
!
      implicit none
!
      character*8 lakonf(*)
!
      integer nef,ipkonf(*),kon(*),nface,ielfa(4,*),ipnei(*),neiel(*),
     &  ifaceq(8,6),i,j,k,indexe,kflag,index1,index2,j1,j2,nope,
     &  nodes(4),iel1,iel2,iel3,iface,indexf,neifa(*),nf(5),ifacet(7,4),
     &  ifacew(8,5),numfaces,ied4(2,6),ied6(2,9),ied8(2,12)
!
      real*8 co(3,*),coel(3,*),cofa(3,*),area(*),xxn(3,*),xxi(3,*),
     &  xle(*),xlen(*),xlet(*),xrlfa(3,*),cosa(*),xsj2(3),xi,et,
     &  shp2(7,4),xs2(3,7),xl2(3,8),xl13,volume(*),dxsj2,xl(3,8),
     &  xxj(3,*),cosb(*),dmin,vel(nef,0:5)
!
!     nodes belonging to the cell faces
!
      data ifaceq /4,3,2,1,11,10,9,12,
     &            5,6,7,8,13,14,15,16,
     &            1,2,6,5,9,18,13,17,
     &            2,3,7,6,10,19,14,18,
     &            3,4,8,7,11,20,15,19,
     &            4,1,5,8,12,17,16,20/
      data ifacet /1,3,2,7,6,5,11,
     &             1,2,4,5,9,8,12,
     &             2,3,4,6,10,9,13,
     &             1,4,3,8,10,7,14/
      data ifacew /1,3,2,9,8,7,0,0,
     &             4,5,6,10,11,12,0,0,
     &             1,2,5,4,7,14,10,13,
     &             2,3,6,5,8,15,11,14,
     &             4,6,3,1,12,15,9,13/
      data nf /3,3,4,4,4/
      data ied4 /1,2,2,3,3,1,1,4,2,4,3,4/
      data ied6 /1,2,2,3,3,1,4,5,5,6,6,4,1,4,2,5,3,6/
      data ied8 /1,2,2,3,3,4,4,1,5,6,6,7,7,8,8,1,1,5,2,6,3,7,4,8/
!
!     coordinates of the center of the cells
!
      do i=1,nef
         if(ipkonf(i).lt.0) cycle
         if(lakonf(i)(1:1).ne.'F') cycle
         indexe=ipkonf(i)
         if(lakonf(i)(4:4).eq.'8') then
            nope=8
         else if(lakonf(i)(4:4).eq.'6') then
            nope=6
         else
            nope=4
         endif
         do j=1,3
            do k=1,nope
               coel(j,i)=coel(j,i)+co(j,kon(indexe+k))
            enddo
            coel(j,i)=coel(j,i)/nope
         enddo
      enddo
!
      kflag=2
!
!     loop over all faces
!
      do i=1,nface
         iel1=ielfa(1,i)
         indexe=ipkonf(iel1)
         j1=ielfa(4,i)
         if(lakonf(iel1)(4:4).eq.'8') then
!
!           hexahedral element
!
!           coordinates of the face centers
!
            do j=1,4
               nodes(j)=kon(indexe+ifaceq(j,j1))
               do k=1,3
                  xl2(k,j)=co(k,nodes(j))
                  cofa(k,i)=cofa(k,i)+xl2(k,j)
               enddo
            enddo
            do k=1,3
               cofa(k,i)=cofa(k,i)/4.d0
            enddo
!
            xi=0.d0
            et=0.d0
            call shape4q(xi,et,xl2,xsj2,xs2,shp2,kflag)
!
!           area of the face
!
            dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
     &                    xsj2(3)*xsj2(3))
            area(i)=4.d0*dxsj2
!
!           normal and xi-vector on face viewed from cell 1
!
            index1=ipnei(iel1)+j1
            do k=1,3
               xxn(k,index1)=xsj2(k)/dxsj2
               xxi(k,index1)=cofa(k,i)-coel(k,iel1)
            enddo
!
!           distance from face center to the center of cell 1
!
            xle(index1)=dsqrt(xxi(1,index1)**2+xxi(2,index1)**2+
     &                         xxi(3,index1)**2)
            do k=1,3
               xxi(k,index1)=xxi(k,index1)/xle(index1)
            enddo
!
!           angle between the normal and the xi-vector
!
            cosa(index1)=xxn(1,index1)*xxi(1,index1)+
     &                   xxn(2,index1)*xxi(2,index1)+
     &                   xxn(3,index1)*xxi(3,index1)
!
            iel2=ielfa(2,i)
!
!           check whether there is an adjacent cell
!
            if(iel2.ne.0) then
               index2=ipnei(iel2)
               do j2=1,6
                  index2=index2+1
                  if(neiel(index2).eq.iel1) exit
               enddo
!
!              normal and xi-vector on face viewed from cell 2
!
               do k=1,3
                  xxi(k,index2)=cofa(k,i)-coel(k,iel2)
                  xxn(k,index2)=-xxn(k,index1)
               enddo
               xle(index2)=dsqrt(xxi(1,index2)**2+xxi(2,index2)**2+
     &              xxi(3,index2)**2)
               do k=1,3
                  xxi(k,index2)=xxi(k,index2)/xle(index2)
               enddo
!
!              angle between the normal and the xi-vector: xxn.xxi
!
               cosa(index2)=xxn(1,index2)*xxi(1,index2)+
     &                      xxn(2,index2)*xxi(2,index2)+
     &                      xxn(3,index2)*xxi(3,index2)
!
!              distance from the face center to the center of the
!              adjacent cell
!
               xlen(index1)=xle(index2)
               xlen(index2)=xle(index1)
!
               do k=1,3
                  xxj(k,index2)=coel(k,iel1)-coel(k,iel2)
               enddo
!
!              distance between the cell center and the center of the
!              adjacent cell
!
               xlet(index1)=dsqrt(xxj(1,index2)**2+xxj(2,index2)**2
     &                                            +xxj(3,index2)**2)
               xlet(index2)=xlet(index1)
!
!              xxj is the unit vector connecting neighboring cell centers
!
               do k=1,3
                  xxj(k,index2)=xxj(k,index2)/xlet(index2)
                  xxj(k,index1)=-xxj(k,index2)
               enddo
!
!              xxn.xxj
!
               cosb(index2)=xxn(1,index1)*xxj(1,index1)+
     &                      xxn(2,index1)*xxj(2,index1)+
     &                      xxn(3,index1)*xxj(3,index1)
               cosb(index1)=cosb(index2)
!
               xrlfa(1,i)=xle(index2)/(xle(index1)+xle(index2))
               xrlfa(2,i)=xle(index1)/(xle(index1)+xle(index2))
            else
!
!              xxi and xxj coincide
!
               do k=1,3
                  xxj(k,index1)=xxi(k,index1)
               enddo
               cosb(index1)=cosa(index1)
!
!              external face: determining the cell next to the
!              adjacent cell
!
               iel3=ielfa(3,i)
               if(iel3.eq.0) cycle
               xl13=dsqrt((coel(1,iel1)-coel(1,iel3))**2+
     &                    (coel(2,iel1)-coel(2,iel3))**2+
     &                    (coel(3,iel1)-coel(3,iel3))**2)
               xrlfa(1,i)=(xl13+xle(index1))/xl13
               xrlfa(3,i)=1.d0-xrlfa(1,i)
            endif
         else if(lakonf(iel1)(4:4).eq.'6') then
!
!           wedge element
!
!           coordinates of the face centers
!
            do j=1,nf(j1)
               nodes(j)=kon(indexe+ifacew(j,j1))
               do k=1,3
                  xl2(k,j)=co(k,nodes(j))
                  cofa(k,i)=cofa(k,i)+xl2(k,j)
               enddo
            enddo
            do k=1,3
               cofa(k,i)=cofa(k,i)/nf(j1)
            enddo
!
            xi=0.d0
            et=0.d0
            if(nf(j1).eq.3) then
               call shape3tri(xi,et,xl2,xsj2,xs2,shp2,kflag)
            else
               call shape4q(xi,et,xl2,xsj2,xs2,shp2,kflag)
            endif
!
!           area of the face
!
            dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
     &                    xsj2(3)*xsj2(3))
            if(nf(j1).eq.3) then
               area(i)=dxsj2/2.d0
            else
               area(i)=4.d0*dxsj2
            endif
!
!           normal and xi-vector on face viewed from cell 1
!
            index1=ipnei(iel1)+j1
            do k=1,3
               xxn(k,index1)=xsj2(k)/dxsj2
               xxi(k,index1)=cofa(k,i)-coel(k,iel1)
            enddo
!
!           distance from face center to the center of cell 1
!
            xle(index1)=dsqrt(xxi(1,index1)**2+xxi(2,index1)**2+
     &                         xxi(3,index1)**2)
            do k=1,3
               xxi(k,index1)=xxi(k,index1)/xle(index1)
            enddo
!
!           angle between the normal and the xi-vector
!
            cosa(index1)=xxn(1,index1)*xxi(1,index1)+
     &                   xxn(2,index1)*xxi(2,index1)+
     &                   xxn(3,index1)*xxi(3,index1)
!
            iel2=ielfa(2,i)
!
!           check whether there is an adjacent cell
!
            if(iel2.ne.0) then
               index2=ipnei(iel2)
               do j2=1,5
                  index2=index2+1
                  if(neiel(index2).eq.iel1) exit
               enddo
!
!              normal and xi-vector on face viewed from cell 2
!
               do k=1,3
                  xxi(k,index2)=cofa(k,i)-coel(k,iel2)
                  xxn(k,index2)=-xxn(k,index1)
               enddo
               xle(index2)=dsqrt(xxi(1,index2)**2+xxi(2,index2)**2+
     &              xxi(3,index2)**2)
               do k=1,3
                  xxi(k,index2)=xxi(k,index2)/xle(index2)
               enddo
!
!              angle between the normal and the xi-vector: xxn.xxi
!
               cosa(index2)=xxn(1,index2)*xxi(1,index2)+
     &                      xxn(2,index2)*xxi(2,index2)+
     &                      xxn(3,index2)*xxi(3,index2)
!
!              distance from the face center to the center of the
!              adjacent cell
!
               xlen(index1)=xle(index2)
               xlen(index2)=xle(index1)
!
               do k=1,3
                  xxj(k,index2)=coel(k,iel1)-coel(k,iel2)
               enddo
!
!              distance between the cell center and the center of the
!              adjacent cell
!
               xlet(index1)=dsqrt(xxj(1,index2)**2+xxj(2,index2)**2
     &                                            +xxj(3,index2)**2)
               xlet(index2)=xlet(index1)
!
!              xxj is the unit vector connecting neighboring cell centers
!
               do k=1,3
                  xxj(k,index2)=xxj(k,index2)/xlet(index2)
                  xxj(k,index1)=-xxj(k,index2)
               enddo
!
!              xxn.xxj
!
               cosb(index2)=xxn(1,index1)*xxj(1,index1)+
     &                      xxn(2,index1)*xxj(2,index1)+
     &                      xxn(3,index1)*xxj(3,index1)
               cosb(index1)=cosb(index2)
!
               xrlfa(1,i)=xle(index2)/(xle(index1)+xle(index2))
               xrlfa(2,i)=xle(index1)/(xle(index1)+xle(index2))
            else
!
!              xxi and xxj coincide
!
               do k=1,3
                  xxj(k,index1)=xxi(k,index1)
               enddo
               cosb(index1)=cosa(index1)
!
!              external face: determining the cell next to the
!              adjacent cell
!
               iel3=ielfa(3,i)
               if(iel3.eq.0) cycle
               xl13=dsqrt((coel(1,iel1)-coel(1,iel3))**2+
     &                    (coel(2,iel1)-coel(2,iel3))**2+
     &                    (coel(3,iel1)-coel(3,iel3))**2)
               xrlfa(1,i)=(xl13+xle(index1))/xl13
               xrlfa(3,i)=1.d0-xrlfa(1,i)
            endif
         else
!
!           tetrahedral element
!
!           coordinates of the face centers
!
            do j=1,3
               nodes(j)=kon(indexe+ifacet(j,j1))
               do k=1,3
                  xl2(k,j)=co(k,nodes(j))
                  cofa(k,i)=cofa(k,i)+xl2(k,j)
               enddo
            enddo
            do k=1,3
               cofa(k,i)=cofa(k,i)/3.d0
            enddo
!
            xi=0.d0
            et=0.d0
            call shape3tri(xi,et,xl2,xsj2,xs2,shp2,kflag)
!
!           area of the face
!
            dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
     &                    xsj2(3)*xsj2(3))
            area(i)=dxsj2/2.d0
!
!           normal and xi-vector on face viewed from cell 1
!
            index1=ipnei(iel1)+j1
            do k=1,3
               xxn(k,index1)=xsj2(k)/dxsj2
               xxi(k,index1)=cofa(k,i)-coel(k,iel1)
            enddo
!
!           distance from face center to the center of cell 1
!
            xle(index1)=dsqrt(xxi(1,index1)**2+xxi(2,index1)**2+
     &                         xxi(3,index1)**2)
            do k=1,3
               xxi(k,index1)=xxi(k,index1)/xle(index1)
            enddo
!
!           angle between the normal and the xi-vector
!
            cosa(index1)=xxn(1,index1)*xxi(1,index1)+
     &                   xxn(2,index1)*xxi(2,index1)+
     &                   xxn(3,index1)*xxi(3,index1)
!
            iel2=ielfa(2,i)
!
!           check whether there is an adjacent cell
!
            if(iel2.ne.0) then
               index2=ipnei(iel2)
               do j2=1,4
                  index2=index2+1
                  if(neiel(index2).eq.iel1) exit
               enddo
!
!              normal and xi-vector on face viewed from cell 2
!
               do k=1,3
                  xxi(k,index2)=cofa(k,i)-coel(k,iel2)
                  xxn(k,index2)=-xxn(k,index1)
               enddo
               xle(index2)=dsqrt(xxi(1,index2)**2+xxi(2,index2)**2+
     &              xxi(3,index2)**2)
               do k=1,3
                  xxi(k,index2)=xxi(k,index2)/xle(index2)
               enddo
!
!              angle between the normal and the xi-vector: xxn.xxi
!
               cosa(index2)=xxn(1,index2)*xxi(1,index2)+
     &                      xxn(2,index2)*xxi(2,index2)+
     &                      xxn(3,index2)*xxi(3,index2)
!
!              distance from the face center to the center of the
!              adjacent cell
!
               xlen(index1)=xle(index2)
               xlen(index2)=xle(index1)
!
               do k=1,3
                  xxj(k,index2)=coel(k,iel1)-coel(k,iel2)
               enddo
!
!              distance between the cell center and the center of the
!              adjacent cell
!
               xlet(index1)=dsqrt(xxj(1,index2)**2+xxj(2,index2)**2
     &                                            +xxj(3,index2)**2)
               xlet(index2)=xlet(index1)
!
!              xxj is the unit vector connecting neighboring cell centers
!
               do k=1,3
                  xxj(k,index2)=xxj(k,index2)/xlet(index2)
                  xxj(k,index1)=-xxj(k,index2)
               enddo
!
!              xxn.xxj
!
               cosb(index2)=xxn(1,index1)*xxj(1,index1)+
     &                      xxn(2,index1)*xxj(2,index1)+
     &                      xxn(3,index1)*xxj(3,index1)
               cosb(index1)=cosb(index2)
!
               xrlfa(1,i)=xle(index2)/(xle(index1)+xle(index2))
               xrlfa(2,i)=xle(index1)/(xle(index1)+xle(index2))
            else
!
!              xxi and xxj coincide
!
               do k=1,3
                  xxj(k,index1)=xxi(k,index1)
               enddo
               cosb(index1)=cosa(index1)
!
!              external face: determining the cell next to the
!              adjacent cell
!
               iel3=ielfa(3,i)
               if(iel3.eq.0) cycle
               xl13=dsqrt((coel(1,iel1)-coel(1,iel3))**2+
     &                    (coel(2,iel1)-coel(2,iel3))**2+
     &                    (coel(3,iel1)-coel(3,iel3))**2)
               xrlfa(1,i)=(xl13+xle(index1))/xl13
               xrlfa(3,i)=1.d0-xrlfa(1,i)
            endif
         endif
c         write(*,*)
      enddo
!
!     calculation of the volume of the elements
!
      do i=1,nef
         if(ipkonf(i).lt.0) cycle
         if(lakonf(i)(1:1).ne.'F') cycle
         indexf=ipnei(i)
         volume(i)=0.d0
         if(lakonf(i)(4:4).eq.'8') then
            numfaces=6
         elseif(lakonf(i)(4:4).eq.'6') then
            numfaces=5
         else
            numfaces=4
         endif
         do j=1,numfaces
            iface=neifa(indexf+j)
            volume(i)=volume(i)+
     &            area(iface)*cofa(1,iface)*xxn(1,indexf+j)
         enddo
c         write(*,*) 'initialcfd volume ',i,volume(i)
      enddo
!
!     calculation of the minimum length within the cells
!
      dmin=1.d30
      do i=1,nef
         indexe=ipkonf(i)
         read(lakonf(i)(4:4),'(i1)') nope
         do j=1,nope
            do k=1,3
               xl(k,j)=co(k,kon(indexe+j))
            enddo
         enddo
         if(nope.eq.4) then
            do j=1,6
               dmin=min(dmin,(xl(1,ied4(1,j))-xl(1,ied4(2,j)))**2+
     &                       (xl(2,ied4(1,j))-xl(2,ied4(2,j)))**2+
     &                       (xl(3,ied4(1,j))-xl(3,ied4(2,j)))**2)
            enddo
         elseif(nope.eq.6) then
            do j=1,9
               dmin=min(dmin,(xl(1,ied6(1,j))-xl(1,ied6(2,j)))**2+
     &                       (xl(2,ied6(1,j))-xl(2,ied6(2,j)))**2+
     &                       (xl(3,ied6(1,j))-xl(3,ied6(2,j)))**2)
            enddo
         else
            do j=1,12
               dmin=min(dmin,(xl(1,ied8(1,j))-xl(1,ied8(2,j)))**2+
     &                       (xl(2,ied8(1,j))-xl(2,ied8(2,j)))**2+
     &                       (xl(3,ied8(1,j))-xl(3,ied8(2,j)))**2)
            enddo
         endif
      enddo
      dmin=dsqrt(dmin)
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


! { dg-do compile }
!
! PR 40089: Public type with public component which has a private type
!
! Original test case by Juergen Reuter <reuter@physik.uni-freiburg.de>
! Adapted by Janus Weil <janus@gcc.gnu.org>

module m

  implicit none
  private

  public :: public_t

  type :: private_t
    integer :: i
  end type

  type :: public_t
     type(private_t), pointer :: public_comp_with_private_type
     procedure(ifc) , nopass, pointer :: ppc
  end type

  abstract interface
     integer function ifc ()
     end function
  end interface

end module m

program test
use m
implicit none
type(public_t) :: x
integer :: j
j = x%ppc()
end


! { dg-do compile }
module c_assoc_2
  use, intrinsic :: iso_c_binding, only: c_ptr, c_associated

contains
  subroutine sub0(my_c_ptr) bind(c)
    type(c_ptr), value :: my_c_ptr
    type(c_ptr), pointer :: my_c_ptr_2
    integer :: my_integer

    if(.not. c_associated(my_c_ptr)) then
       call abort()
    end if

    if(.not. c_associated(my_c_ptr, my_c_ptr)) then
       call abort()
    end if

    if(.not. c_associated(my_c_ptr, my_c_ptr, my_c_ptr)) then ! { dg-error "Too many arguments in call" }
       call abort()
    end if

    if(.not. c_associated()) then ! { dg-error "Missing actual argument" }
       call abort()
    end if

    if(.not. c_associated(my_c_ptr_2)) then
       call abort()
    end if

    if(.not. c_associated(my_integer)) then ! { dg-error "shall have the type TYPE.C_PTR. or TYPE.C_FUNPTR." }
       call abort()
    end if
  end subroutine sub0

end module c_assoc_2


text:
!<types>

!<typeblock>

  ! Structure which collects information for function evaluation routines.
  ! Remarks:
  ! a) If the structure describes a number of points where to evaluate
  !    without an underlying mesh, there is
  !        nelements = 1
  !        npointPerElement = #points
  !        p_Dpoints = coordinates
  !    and p_rtriangulation => null().
  ! b) If the structure describes only one point where to evaluate, there is
  !        nelements = 1
  !        npointPerElement = 1
  !        p_Dpoints = coordinates.
  type t_paramsSimEval

    ! <!-- ############################################### -->
    ! <!-- Mandatory information, must always be available -->
    ! <!-- ############################################### -->

    ! Number of elements.
    integer :: nelements = 0

    ! Number of points per element.
    integer :: npointsPerElement = 0

    ! Array of all points on all the elements where coefficients are needed.
    ! DIMENSION(#dimensions,npointsPerElement,nelements)
    real(DP), dimension(:,:,:), pointer :: p_Dpoints => null()

    ! <!-- ###################################### -->
    ! <!-- Optional information, mostly available -->
    ! <!-- ###################################### -->

    ! Pointer to the underlying boundary structure.
    ! If not available, this points to null().
    type(t_boundary), pointer :: p_rboundary => null()

    ! Pointer to the underlying triangulation
    ! If not available, this points to null().
    type(t_triangulation), pointer :: p_rtriangulation => null()

    ! List of elements.
    ! If not available, this points to null().
    integer, dimension(:), pointer :: p_Ielements => null()

    ! Coordinates of the evaluation points in the reference
    ! coordinate system.
    ! If this information is not available, this parameter points to null().
    ! DIMENSION(#dimensions,npointsPerElement,nelements)
    real(DP), dimension(:,:,:), pointer :: p_DpointsRef => null()

    ! The Jacobian matrix of the mapping between the reference and each
    ! real element, for all points on all elements in progress.
    ! If this information is not available, this parameter points to null().
    ! DIMENSION(#dimensions*#dimensions,npointsPerElement,nelements)
    real(DP), dimension(:,:,:), pointer :: p_Djac => null()

    ! The Jacobian determinant of the mapping of each point from the
    ! reference element to each real element in progress.
    ! If this information is not available, this parameter points to null().
    ! DIMENSION(npointsPerElement,nelements)
    real(DP), dimension(:,:), pointer :: p_Ddetj => null()

    ! <!-- ################################ -->
    ! <!-- Evaluation on the boundary in 2D -->
    ! <!-- ################################ -->

    ! Number of the current boundary component, or =0 if not available.
    ! All elements in p_Ielements refer to this boundary component.
    integer :: ibct = 0

    ! Parameter values of the points on the boundary
    ! If this information is not available, this parameter points to null().
    real(dp), dimension(:), pointer :: p_DpointsPar => null()

    ! Boundary region that is currently being processed.
    ! If this information is not available, this parameter points to null().
    type(t_boundaryRegion), pointer :: rboundaryRegion => null()

    ! Mesh region that is currently being processed.
    ! If this information is not available, this parameter points to null().
    type(t_meshRegion), pointer :: p_rmeshRegion => null()

    ! <!-- ########################## -->
    ! <!-- Evaluation of linear forms -->
    ! <!-- ########################## -->

    ! The bilinear form which is currently being evaluated.
    ! If this information is not available, this parameter points to null().
    type(t_linearForm), pointer :: p_rlinearForm => null()

    ! <!-- ############################ -->
    ! <!-- Evaluation of bilinear forms -->
    ! <!-- ############################ -->

    ! The bilinear form which is currently being evaluated.
    ! If this information is not available, this parameter points to null().
    type(t_bilinearForm), pointer :: p_rbilinearForm => null()

    ! <!-- ############################# -->
    ! <!-- Evaluation of trilinear forms -->
    ! <!-- ############################# -->

    ! The trilinear form which is currently being evaluated.
    ! If this information is not available, this parameter points to null().
    type(t_trilinearForm), pointer :: p_rtrilinearForm => null()
  end type

!</typeblock>

!</types>

!<constants>
!<constantblock description="Task identifiers for callback function evaluation">
  ! Do nothing - this is just a dummy.
  integer, parameter :: FUNC_TASK_NONE                  = 0

  ! Evaluate the function in a given set of points.
  ! Results must be saved to Dvalues(:,:,1).
  integer, parameter :: FUNC_TASK_EVAL_FUNC             = 1

  ! Evaluate one component of a derivative.
  ! The variable itag is a derivative quantifier (DER_DERIV_xxxx) that
  ! defines the component that should be evaluated.
  !   Dvalues(:,:,1)=derivatives identified by itag
  integer, parameter :: FUNC_TASK_EVAL_DERIV            = 2

  ! Evaluate the full 1st order derivative of the function.
  ! Results must be saved to:
  !   Dvalues(:,:,1)=x-derivatives,
  !   Dvalues(:,:,2)=y-derivatives (if exists),
  !   Dvalues(:,:,3)=z-derivatives (if exists).
  integer, parameter :: FUNC_TASK_EVAL_DERIV1FULL       = 3

  ! Evaluate the function in a given set of points on the boundary.
  ! Results must be saved to Dvalues(:,:,1).
  integer, parameter :: FUNC_TASK_EVAL_FUNC_BD          = 4

  ! Evaluate one component of a derivativeon the boundary.
  ! The variable itag is a derivative quantifier (DER_DERIV_xxxx) that
  ! defines the component that should be evaluated.
  !   Dvalues(:,:,1)=derivatives identified by itag
  integer, parameter :: FUNC_TASK_EVAL_DERIV_BD         = 5

  ! Evaluate the 1st derivative of the function on the boundary.
  ! Results must be saved to:
  !   Dvalues(:,:,1)=x-derivatives,
  !   Dvalues(:,:,2)=y-derivatives (if exists),
  !   Dvalues(:,:,3)=z-derivatives (if exists).
  integer, parameter :: FUNC_TASK_EVAL_DERIV1_BD        = 6

  ! Evaluate the coefficients in front of the terms of a linear form.
  ! Results must be saved to Dvalues(:,:,1..nterms) with nterms the number
  ! of terms in the linear form.
  integer, parameter :: FUNC_TASK_EVAL_LINFCOEFF        = 7

  ! Evaluate the coefficients in front of the terms of a bilinear form.
  ! Results must be saved to Dvalues(:,:,1..nterms) with nterms the number
  ! of terms in the bilinear form.
  integer, parameter :: FUNC_TASK_EVAL_BILFCOEFF        = 8

  ! Evaluate the coefficients in front of the terms of a trilinear form.
  ! Results must be saved to Dvalues(:,:,1..nterms) with nterms the number
  ! of terms in the trilinear form.
  integer, parameter :: FUNC_TASK_EVAL_TRILFCOEFF       = 9
!</constantblock>
!</constants>

  interface

  !<subroutine>

    subroutine ffunctionSimEval (ctask,itag,rparams,Dvalues,rcollection)

  !<description>
    ! This interface defines a general implementation of a callback routine
    ! used for the evaluation of functions in the kernel. It is used for matrix
    ! and vector assembly as well as for the calculation of itegrals.
    !
    ! The purpose of this routine is to calculate the values of a (usually
    ! analytically given) function in a number of points on a set of elements.
    ! The stricture rparams defines all parameters where to evaluate.
    ! The points where to evaluate can e.g. be accessed by:
    !    dx = rparams%p_Dpoints(1,ipoint,ielement)
    !    dy = rparams%p_Dpoints(2,ipoint,ielement)
    ! -> x- and y-coordinates of the point ipoint on element ielement.
    ! Depending on the purpose for which this routine is used,
    ! the parameters in rparams may provide additional information
    ! (like a triangulation, boundary, element list etc.).
    !
    ! The values that are calculated in this routine have to be written
    ! to the Dvalues(:,:,:) array. The last dimension of Dvalues defines
    ! here the number of a 'term'; for standard callback functions, 1 has
    ! to be used as an index here, e.g.
    !   Dvalues(:,:,1) = 16.0_DP * rparams%p_Dpoints(2,:,:) * (1.0_DP-rparams%p_Dpoints(2,:,:))
    ! returns a standard parabolic profile.
    !
    ! For more complicated function evaluations like the evaluation of a
    ! bilinear form, this routine may have to calculate multiple terms,
    ! like the coefficients of multiple additive integral terms; In this
    ! case, iterm in Dvalues(:,:,iterm) allows to access each term to
    ! calculate.
  !</description>

    use fsystem
    use collection

  !<input>
    ! A task identifier that specifies in which way the function should return
    ! its values. One of the FUNC_TASK_xxxx constants.
    ! The standard value is e.g. FUNC_TASK_EVAL_FUNC which tells the routine just to
    ! return function values.
    integer, intent(in) :: ctask
  !</input>

  !<inputoutput>
    ! A tag variable. The exact meaning depends on the assembly routine.
    ! If not used, 0 is passed here.
    integer, intent(inout) :: itag

    ! Parameter structure defining where to evaluate
    type(t_paramsSimEval), intent(inout) :: rparams

    ! Optional: A collection structure to provide additional
    ! information to the coefficient routine.
    type(t_collection), intent(inout), optional :: rcollection
  !</inputoutput>

  !<output>
    ! Result values. #terms varies, depending on ctask.
    ! DIMENSION(#points per element, #elements, #terms)
    real(DP), dimension(:,:,:), intent(out) :: Dvalues
  !</output>

  !</subroutine>

    end subroutine

  end interface

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


! Ensure that each parallel reduction variable as a copy or pcopy
! data clause.

! { dg-additional-options "-fdump-tree-gimple" }

program test
  implicit none
  integer :: v1, v2

  !$acc parallel reduction(+:v1,v2)
  !$acc end parallel

  !$acc parallel reduction(+:v1,v2) copy(v1,v2)
  !$acc end parallel

  !$acc parallel reduction(+:v1,v2) pcopy(v1,v2)
  !$acc end parallel

  !$acc parallel reduction(+:v1,v2) present(v1,v2)
  !$acc end parallel

  !$acc parallel reduction(+:v1,v2) copyin(v1,v2) ! { dg-warning "incompatible data clause" }
  !$acc end parallel

  !$acc parallel reduction(+:v1,v2) pcopyin(v1,v2) ! { dg-warning "incompatible data clause" }
  !$acc end parallel

  !$acc parallel reduction(+:v1,v2) copyout(v1,v2) ! { dg-warning "incompatible data clause" }
  !$acc end parallel

  !$acc parallel reduction(+:v1,v2) pcopyout(v1,v2) ! { dg-warning "incompatible data clause" }
  !$acc end parallel

  !$acc parallel reduction(+:v1,v2) create(v1,v2) ! { dg-warning "incompatible data clause" }
  !$acc end parallel

  !$acc parallel reduction(+:v1,v2) pcreate(v1,v2) ! { dg-warning "incompatible data clause" }
  !$acc end parallel
end program test

! { dg-final { scan-tree-dump-times "map.tofrom:v1" 8 "gimple" } }
! { dg-final { scan-tree-dump-times "map.tofrom:v2" 8 "gimple" } }
! { dg-final { scan-tree-dump-times "map.force_tofrom:v1" 1 "gimple" } }
! { dg-final { scan-tree-dump-times "map.force_tofrom:v2" 1 "gimple" } }
! { dg-final { scan-tree-dump-times "map.force_present:v1" 1 "gimple" } }
! { dg-final { scan-tree-dump-times "map.force_present:v2" 1 "gimple" } }


text:
      subroutine indices (low, high, datav, datar, numr)
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
!***********************************************************************
!D1
!D1 PURPOSE
!D1
!D1 Find indices of appropriate sigma and omega curves to be used in
!D1 interpolation (log value method).
!D1
!***********************************************************************
!D2
!D2 REVISION HISTORY
!D2
!D2 FEHM Version 2.10 [10086-STN-2.10-00]
!D2 Initial implementation: 05-FEB-99, Programmer: Z. Dash
!D2 Modification for data revision: 16-AUG-99
!D2
!D2 $Log:   /pvcs.config/fehm90/src/indices.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:43:16   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:08:00   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:09:56   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:30:58   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 12:03:10   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!***********************************************************************
!D3
!D3 REQUIREMENTS TRACEABILITY
!D3
!D3 2.3.5 Cell-based particle-tracking module
!D3 2.3.6 Streamline particle-tracking module
!D3
!***********************************************************************
!D4
!D4 SPECIAL COMMENTS AND REFERENCES
!D4
!D4 Requirements from SDN: 10086-RD-2.20-00
!D4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
!D4   FEHM Application Version 2.20
!D4
!***********************************************************************

      implicit none

      real*8 datav, datar(*)
      integer low, high, numr, midpoint

C Find the T vs C curve indices

      if (datav .le. datar(1)) then
         low = 1
         high = low
         return
      else if (datav .ge. datar(numr)) then
         low = numr
         high = low
         return
      else
         low = 1
         high = numr
 100     midpoint = (low + high) / 2
         if (datav .eq. datar(midpoint)) then
            low = midpoint
            high = low
            return
         else if (datav .lt. datar(midpoint)) then
            high = midpoint
         else
            low = midpoint
         endif
         if (high - low .eq. 1) return
         goto 100
      endif

      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
      FUNCTION VGG(B,OMC,IG,IP)
C
C     'VGG' IS A FUNCTION WHICH EVALUATES THE CENTRIFUGAL POTENTIAL
C     FOR THE GRAVITY GRADIENT BODY.
C
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 IG(3,3),IP
C
      DIMENSION B(3,3),OMC(3),OC(3),OCIG(3)
C
C
      CALL MULTM(B,OMC,OC,3,1,3)
      DO 10 I=1,3
      OCIG(I)=0.0D0
      DO 10 J=1,3
   10 OCIG(I)=OCIG(I) + OC(J)*IG(J,I)
C
      OCSQR=0.0D0
      SUM=0.0D0
C
      DO 20 I=1,3
      OCSQR=OCSQR + OC(I)*OC(I)
   20 SUM=SUM + OCIG(I)*OC(I)
      VGG=(SUM - OCSQR*IP)/2
C
      RETURN
      END

ext: for


text:
      SUBROUTINE XDXDCL(XDXD,ZLK,K,M)
C
C     'XDXDCL' CALCULATES THE SQUARES OF THE RELATIVE VELOCITIES OF
C     AN ELEMENT IN THE ELEMENT FRAME.
C
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER*4 P,Q,R
C
      COMMON/COMALP/ SZ02(10),SZ03(10),SZ04(10),SZ12(3,10),SZ13(3,10),
     .               SZ14(3,10),SZ15(3,10),SZ16(3,10),SZ21(9,10),
     .               SZ22(9,10),SZ23(9,10),SZ25(9,10),
     .               SZ26(9,10),SZ27(9,10),SZ28(9,10),SZ31(27,10),
     .               SZ32(27,10),SZ33(27,10),SZ34(27,10),SZ35(27,10),
     .               SZ41(81,10),SZ42(81,10),SZ43(81,10)
C
      COMMON/DEBUG3/ISWTCH
C
      COMMON/IPOOL1/ IGRAV,IDAMP,IK,K1,ITIM,IAB,IAPS,IBB,IBPS,NK(10),
     .               LK(10),LLK(10)
C
C
      DIMENSION XDXD(3,3),ZLK(10)
C
C
      IF(ISWTCH.EQ.0) WRITE(6,20002) K,M
      SUM7=0.0D0
      SUM8=0.0D0
      SUM11=0.0D0
      SUM12=0.0D0
C
      DO 40 N=1,M
      SUM5=0.0D0
      SUM6=0.0D0
      SUM9=0.0D0
      SUM10=0.0D0
C
      DO 30 P=1,M
      SUM3=0.0D0
      SUM4=0.0D0
C
      DO 20 Q=1,M
      SUM1=0.0D0
      SUM2=0.0D0
C
      DO 10 R=1,M
      A1=FUNA(K,K1,R)
      B1=FUNB(K,K1,R)
      IR=27*(N-1)+9*(P-1)+3*(Q-1)+R
      SUM1=SUM1+B1*SZ41(IR,K)
   10 SUM2=SUM2+A1*SZ41(IR,K)
C
      AD1=ADFUN(K,K1,Q)
      BD1=BDFUN(K,K1,Q)
      SUM3=SUM3 + BD1*SUM1
   20 SUM4=SUM4 + AD1*SUM2
C
      A1=FUNA(K,K1,P)
      B1=FUNB(K,K1,P)
      AD1=ADFUN(K,K1,P)
      BD1=BDFUN(K,K1,P)
      IP=3*(N-1)+P
C
      SUM5=SUM5 + A1*(SUM3 + SUM4)
      SUM6=SUM6 + B1*(SUM3 + SUM4)
      SUM9=SUM9+AD1*SZ21(IP,K)
   30 SUM10=SUM10+BD1*SZ21(IP,K)
C
      AD1=ADFUN(K,K1,N)
      BD1=BDFUN(K,K1,N)
C
      SUM7=SUM7  + AD1*SUM5
      SUM8=SUM8  + BD1*SUM6
      SUM11=SUM11 + AD1*SUM9
   40 SUM12=SUM12 + BD1*SUM10
C
      XDXD(1,1)=(SUM7 + SUM8)/(ZLK(K)*ZLK(K))
      XDXD(2,2)=SUM11
      XDXD(3,3)=SUM12
C
      IF(ISWTCH.EQ.0) WRITE(6,10000) ((XDXD(I,J),J=1,3),I=1,3)
C
      RETURN
C
10000 FORMAT('0XDX(K,I,J)',//3G15.5)
20002 FORMAT('0',5X,'XDXDCL ',2I4)
20010 FORMAT('0',2X,'J2 = SUM OF ZZT1')
C
      END


text:
C-----------------------------------------------------------------------
      SUBROUTINE SPFFTE(IMAX,INCW,INCG,KMAX,W,G,IDIR,AFFT)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:  SPFFTE     PERFORM MULTIPLE FAST FOURIER TRANSFORMS
C   PRGMMR: IREDELL       ORG: W/NMC23       DATE: 96-02-20
C
C ABSTRACT: THIS SUBPROGRAM PERFORMS MULTIPLE FAST FOURIER TRANSFORMS
C           BETWEEN COMPLEX AMPLITUDES IN FOURIER SPACE AND REAL VALUES
C           IN CYCLIC PHYSICAL SPACE.
C           SUBPROGRAM SPFFTE MUST BE INVOKED FIRST WITH IDIR=0
C           TO INITIALIZE TRIGONEMETRIC DATA.  USE SUBPROGRAM SPFFT1
C           TO PERFORM AN FFT WITHOUT PREVIOUS INITIALIZATION.
C           THIS VERSION INVOKES THE IBM ESSL FFT.
C
C PROGRAM HISTORY LOG:
C 1998-12-18  IREDELL
C 2012-11-12  MIRVIS -fixing hard-wired types problem on Intel/Linux

C USAGE:    CALL SPFFTE(IMAX,INCW,INCG,KMAX,W,G,IDIR,AFFT)
C
C   INPUT ARGUMENT LIST:
C     IMAX     - INTEGER NUMBER OF VALUES IN THE CYCLIC PHYSICAL SPACE
C                (SEE LIMITATIONS ON IMAX IN REMARKS BELOW.)
C     INCW     - INTEGER FIRST DIMENSION OF THE COMPLEX AMPLITUDE ARRAY
C                (INCW >= IMAX/2+1)
C     INCG     - INTEGER FIRST DIMENSION OF THE REAL VALUE ARRAY
C                (INCG >= IMAX)
C     KMAX     - INTEGER NUMBER OF TRANSFORMS TO PERFORM
C     W        - COMPLEX(INCW,KMAX) COMPLEX AMPLITUDES IF IDIR>0
C     G        - REAL(INCG,KMAX) REAL VALUES IF IDIR<0
C     IDIR     - INTEGER DIRECTION FLAG
C                IDIR=0 TO INITIALIZE TRIGONOMETRIC DATA
C                IDIR>0 TO TRANSFORM FROM FOURIER TO PHYSICAL SPACE
C                IDIR<0 TO TRANSFORM FROM PHYSICAL TO FOURIER SPACE
C     AFFT       REAL(8) (50000+4*IMAX) AUXILIARY ARRAY IF IDIR<>0
C
C   OUTPUT ARGUMENT LIST:
C     W        - COMPLEX(INCW,KMAX) COMPLEX AMPLITUDES IF IDIR<0
C     G        - REAL(INCG,KMAX) REAL VALUES IF IDIR>0
C     AFFT       REAL(8) (50000+4*IMAX) AUXILIARY ARRAY IF IDIR=0
C
C SUBPROGRAMS CALLED:
C   SCRFT        IBM ESSL COMPLEX TO REAL FOURIER TRANSFORM
C   DCRFT        IBM ESSL COMPLEX TO REAL FOURIER TRANSFORM
C   SRCFT        IBM ESSL REAL TO COMPLEX FOURIER TRANSFORM
C   DRCFT        IBM ESSL REAL TO COMPLEX FOURIER TRANSFORM
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C
C REMARKS:
C   THE RESTRICTIONS ON IMAX ARE THAT IT MUST BE A MULTIPLE
C   OF 1 TO 25 FACTORS OF TWO, UP TO 2 FACTORS OF THREE,
C   AND UP TO 1 FACTOR OF FIVE, SEVEN AND ELEVEN.
C
C   IF IDIR=0, THEN W AND G NEED NOT CONTAIN ANY VALID DATA.
C   THE OTHER PARAMETERS MUST BE SUPPLIED AND CANNOT CHANGE
C   IN SUCCEEDING CALLS UNTIL THE NEXT TIME IT IS CALLED WITH IDIR=0.
C
C   THIS SUBPROGRAM IS THREAD-SAFE.
C
C$$$
        IMPLICIT NONE
        INTEGER,INTENT(IN):: IMAX,INCW,INCG,KMAX,IDIR
        REAL,INTENT(INOUT):: W(2*INCW,KMAX)
        REAL,INTENT(INOUT):: G(INCG,KMAX)
        REAL(8),INTENT(INOUT):: AFFT(50000+4*IMAX)
        INTEGER:: INIT,INC2X,INC2Y,N,M,ISIGN,NAUX1,NAUX2,NAUX3
C ==EM==       ^(4)
        REAL:: SCALE
        REAL(8):: AUX2(20000+2*IMAX),AUX3
        INTEGER:: IACR,IARC
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        NAUX1=25000+2*IMAX
        NAUX2=20000+2*IMAX
        NAUX3=1
        IACR=1
        IARC=1+NAUX1
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  INITIALIZATION.
C  FILL AUXILIARY ARRAYS WITH TRIGONOMETRIC DATA
        SELECT CASE(IDIR)
          CASE(0)
            INIT=1
            INC2X=INCW
            INC2Y=INCG
            N=IMAX
            M=KMAX
            ISIGN=-1
            SCALE=1.
            IF(DIGITS(1.).LT.DIGITS(1._8)) THEN
              CALL SCRFT(INIT,W,INC2X,G,INC2Y,N,M,ISIGN,SCALE,
     &                   AFFT(IACR),NAUX1,AUX2,NAUX2,AUX3,NAUX3)
            ELSE
              CALL DCRFT(INIT,W,INC2X,G,INC2Y,N,M,ISIGN,SCALE,
     &                   AFFT(IACR),NAUX1,AUX2,NAUX2)
            ENDIF
            INIT=1
            INC2X=INCG
            INC2Y=INCW
            N=IMAX
            M=KMAX
            ISIGN=+1
            SCALE=1./IMAX
            IF(DIGITS(1.).LT.DIGITS(1._8)) THEN
              CALL SRCFT(INIT,G,INC2X,W,INC2Y,N,M,ISIGN,SCALE,
     &                   AFFT(IARC),NAUX1,AUX2,NAUX2,AUX3,NAUX3)
            ELSE
              CALL DRCFT(INIT,G,INC2X,W,INC2Y,N,M,ISIGN,SCALE,
     &                   AFFT(IARC),NAUX1,AUX2,NAUX2)
            ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  FOURIER TO PHYSICAL TRANSFORM.
          CASE(1:)
            INIT=0
            INC2X=INCW
            INC2Y=INCG
            N=IMAX
            M=KMAX
            ISIGN=-1
            SCALE=1.
            IF(DIGITS(1.).LT.DIGITS(1._8)) THEN
              CALL SCRFT(INIT,W,INC2X,G,INC2Y,N,M,ISIGN,SCALE,
     &                   AFFT(IACR),NAUX1,AUX2,NAUX2,AUX3,NAUX3)
            ELSE
              CALL DCRFT(INIT,W,INC2X,G,INC2Y,N,M,ISIGN,SCALE,
     &                   AFFT(IACR),NAUX1,AUX2,NAUX2)
            ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  PHYSICAL TO FOURIER TRANSFORM.
          CASE(:-1)
            INIT=0
            INC2X=INCG
            INC2Y=INCW
            N=IMAX
            M=KMAX
            ISIGN=+1
            SCALE=1./IMAX
            IF(DIGITS(1.).LT.DIGITS(1._8)) THEN
              CALL SRCFT(INIT,G,INC2X,W,INC2Y,N,M,ISIGN,SCALE,
     &                   AFFT(IARC),NAUX1,AUX2,NAUX2,AUX3,NAUX3)
            ELSE
              CALL DRCFT(INIT,G,INC2X,W,INC2Y,N,M,ISIGN,SCALE,
     &                   AFFT(IARC),NAUX1,AUX2,NAUX2)
            ENDIF
        END SELECT
      END SUBROUTINE

max_stars_repo_name: GMAO-SI-Team/jedi-stack

max_issues_repo_name: GMAO-SI-Team/jedi-stack

max_forks_repo_name: GMAO-SI-Team/jedi-stack


text:
 Take Back Our Democracy

Take Back Our Democracy is an ad hoc group of activists working locally to overturn the http://democracyisforpeople.org/supremecourtruling.php Citizens United v. FEC Supreme Court decision, the one opening the floodgates to unlimited corporate campaign spending.

Contact:  Mary Zhu, 1vox1vote@gmail.com

 Upcoming Events:

We are planning an event to coincide with nationwide actions protesting the third anniversary of the Citizens United decision.  The aim is to pressure politicians to pass a Constitutional amendment that will overrule the Supreme Court.  All concerned individuals are welcome to participate in the planning meeting:

When:  Saturday, December 15, 12:30 pm  2:30 pm
Where:  The Blanchard Community Room, Yolo County Library  Mary L. Stephens Branch, 315 East 14th Street

The meeting will include a national teleconference with http://www.citizen.org Public Citizen and Keith Ellison, cochair of the http://cpc.grijalva.house.gov/ Congressional Progressive Caucus.


text:
      subroutine linsol(b)
************************************************************************
*
*   LINSOL -- Linear equation solver.
*
************************************************************************
      include "parameters.h"

      integer nnod,ntri
      common /mfe010/ nnod,ntri

      double precision s
      integer bandw,xadj,adjncy
      common /mfe040/ s(smax),bandw,xadj(mnod+1),adjncy(2*medg)

      double precision b(*)

      if(.true.) then
*       ----------------------
*        Block banded solver.
*       ----------------------
        call bbslv(nnod,3,bandw,s,  b)

      else
C*       ------------------------
C*        Block envelope solver.
C*       ------------------------
C*       Input RHS vector to BENVPAK and solve.
C        call inrhs(b,  s)
C        call solve(s)
C
C*       Copy solution from S.
C        do 10 j=1,3*nnod
C   10   b(j) = s(j)
        endif

      return
      end

max_stars_repo_name: mfeproject/legacy-gwmfe

max_issues_repo_name: mfeproject/legacy-gwmfe

max_forks_repo_name: mfeproject/legacy-gwmfe


text:
!
! Copyright (C) 2002-2010 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
! ... This module contains functions to check if the code should
! ... be smoothly stopped.
! ... In particular the function check_stop_now returns .TRUE. if
! ... either the user has created a given file or if the
! ... elapsed time is larger than max_seconds
!
!------------------------------------------------------------------------------!
MODULE check_stop
!------------------------------------------------------------------------------!
  !
  USE kinds
  !
  IMPLICIT NONE
  !
  SAVE
  !
  REAL(DP) :: max_seconds = 1.E+7_DP
  !
  LOGICAL, PRIVATE :: tinit = .FALSE.
  !
  REAL(DP) :: init_second
  !
  CONTAINS
     !
     ! ... internal procedures
     !
     !-----------------------------------------------------------------------
     SUBROUTINE check_stop_init()
       !-----------------------------------------------------------------------
       !
       USE input_parameters, ONLY : max_seconds_ => max_seconds
       USE io_global,        ONLY : stdout
       USE io_files,         ONLY : prefix, exit_file
#if defined __TRAP_SIGUSR1
       USE set_signal,       ONLY : signal_trap_init
#endif

       !
       IMPLICIT NONE
       !
       REAL(DP), EXTERNAL :: cclock
       !
       IF ( tinit ) &
          WRITE( UNIT = stdout, &
                 FMT = '(/,5X,"WARNING: check_stop already initialized")' )
       !
       ! ... the exit_file name is set here
       !
       exit_file = TRIM( prefix ) // '.EXIT'
       !
       IF ( max_seconds_ > 0.0_DP ) max_seconds = max_seconds_
       !
       init_second = cclock()
       tinit   = .TRUE.
       !
#if defined __TRAP_SIGUSR1
       CALL signal_trap_init ( )
#endif
       !
       RETURN
       !
     END SUBROUTINE check_stop_init
     !
     !-----------------------------------------------------------------------
     FUNCTION check_stop_now( inunit )
       !-----------------------------------------------------------------------
       !
       USE mp,         ONLY : mp_bcast
       USE mp_global,  ONLY : intra_image_comm
       USE io_global,  ONLY : ionode, ionode_id, meta_ionode, stdout
       USE io_files,   ONLY : exit_file, iunexit
#if defined __TRAP_SIGUSR1
       USE set_signal, ONLY : signal_detected
#endif
       ! KNK_image
       ! USE mp_global, ONLY : mpime, root, world_comm
       !
       IMPLICIT NONE
       !
       INTEGER, OPTIONAL, INTENT(IN) :: inunit
       !
       INTEGER            :: unit
       LOGICAL            :: check_stop_now, tex
       LOGICAL            :: signaled
       REAL(DP)           :: seconds
       REAL(DP), EXTERNAL :: cclock
       !
       !
       ! ... cclock is a C function returning the elapsed solar
       ! ... time in seconds since the Epoch ( 00:00:00 1/1/1970 )
       !
       IF ( .NOT. tinit ) &
          CALL errore( 'check_stop_now', 'check_stop not initialized', 1 )
       !
       unit = stdout
       IF ( PRESENT( inunit ) ) unit = inunit
       !
       check_stop_now = .FALSE.
       !
       signaled = .FALSE.
       !
       ! KNK_image
       ! IF ( mpime == root ) THEN
       IF ( ionode ) THEN
          !
          INQUIRE( FILE = TRIM( exit_file ), EXIST = tex )
          !
          IF ( tex ) THEN
             !
             check_stop_now = .TRUE.
             !
             OPEN( UNIT = iunexit, FILE = TRIM( exit_file ) )
             CLOSE( UNIT = iunexit, STATUS = 'DELETE' )
             !
          ELSE
             !
             seconds = cclock() - init_second
             !
             check_stop_now = ( seconds  >  max_seconds )
             !
          END IF
          !
       END IF
       !
#if defined __TRAP_SIGUSR1
       signaled = signal_detected()
       check_stop_now = check_stop_now .OR. signaled
       tex = tex .OR. signaled
#endif
       !
       ! KNK_image
       ! CALL mp_bcast( check_stop_now, root, world_comm )
       CALL mp_bcast( check_stop_now, ionode_id, intra_image_comm )
       !
       IF ( check_stop_now .AND. meta_ionode ) THEN
          !
          IF ( tex ) THEN
             !
             WRITE( UNIT = unit, &
                    FMT = '(/,5X,"Program stopped by user request")' )
             !
          ELSE
             !
             WRITE( UNIT = unit, &
                    FMT = '(/,5X,"Maximum CPU time exceeded")' )
             WRITE( UNIT = unit, &
                    FMT = '(/,5X,"max_seconds     = ",F10.2)' ) max_seconds
             WRITE( UNIT = unit, &
                    FMT = '(5X,"elapsed seconds = ",F10.2)' ) seconds
             !
          END IF
          !
       END IF
       !
       RETURN
       !
     END FUNCTION check_stop_now
     !
END MODULE check_stop

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
!**********************************************************************
!   This is Steve White's rewrite of Mike Teter's integration routine.
!   Modified by J. Rehr for complex integration.
!   Modified by L. Campbell for Fortran 90
!   The following is a listing of the arguments in the initial function
!   statement:
!      fn    -- routine requires external function statement in MAIN
!      xmin  -- lower limit
!      xmax  -- upper limit
!      abr   -- absolute tolerable error
!      rlr   -- relative tolerable error
!      nsing -- number of singularities or regions requiring
!                   special attention
!      xsing -- array of locations of singularities or endpoints
!                   of special regions
!      error -- output for routine error messages
!      numcal-- the number of times fn was called
!      maxns -- the maximum number of regions being considered simultaneously

double precision function grater(fn,xmin,xmax,abr,rlr,nsing,xsing,error,numcal,maxns)
implicit none
integer :: nsing,numcal,maxns,mx,nstack,ii,jj,kk,icount
parameter (mx=1500)
double precision :: fn,value,valu,fval(3,mx)
double precision :: xmin,xmax,del,del1,dif,frac
double precision :: abr,rlr,error
double precision :: xleft(mx),dx(3),wt(3)
double precision :: wt9(9),xsing(30)
external fn
logical :: atsing
save dx,wt,wt9
data dx/0.1127016653792583,0.5,0.8872983346207417/
data wt/0.277777777777777778,0.4444444444444444444,0.2777777777777777778/
data wt9/0.0616938806304841571,0.108384229110206161,         &
&           0.0398463603260281088,0.175209035316976464,      &
&           0.229732989232610220,0.175209035316976464,       &
&           0.0398463603260281088,0.108384229110206161,      &
&           0.0616938806304841571  /

! nstack is the number of different intervals into which the
! integration region is currently divided. The number of regions can
! grow if more accuracy is needed by dividing the right-most region
! into three regions. The number of regions shrinks when the integral
! over the right-most region is accurate enough, in which case that
! integral is added to the total (stored in grater) and the region
! is removed from consideration (and a new region is the right-most).
  nstack=nsing+1
  maxns = nstack
  error=0.d0
  grater=0.d0
! The array xleft stores the boundary points of the regions.
! The singular points just govern the initial placement of the regions.
  xleft(1)=xmin
  xleft(nsing+2)=xmax
  if(nsing.gt.0) then
    do jj=1,nsing
      xleft(jj+1)=xsing(jj)
    enddo
  endif
! For each region, calculate the function and store at three selected points.
  do ii=1,nstack
    del=xleft(ii+1)-xleft(ii)
    do jj=1,3
      fval(jj,ii)=fn(xleft(ii)+del*dx(jj))
    enddo
  enddo
  numcal = nstack * 3
  do
    if(nstack+3.ge.mx) then
      write(*,*) 'TOO MANY REGIONS'
      stop 0006
    endif
! Divide the rightmost region into three subregions.
    del=xleft(nstack+1)-xleft(nstack)
    xleft(nstack+3)=xleft(nstack+1)
    xleft(nstack+1)=xleft(nstack)+del*dx(1)*2.
    xleft(nstack+2)=xleft(nstack+3)-del*dx(1)*2.
! The three data points already found for the region become the
! middle data points (number 2 in first index of fval) for each region.
    fval(2,nstack+2)=fval(3,nstack)
    fval(2,nstack+1)=fval(2,nstack)
    fval(2,nstack)=fval(1,nstack)
! Now do the integral over the right-most region in two different ways-
! a three point integral (valu) over each of the three subregions
! and a more accurate nine-point integral (value) over whole region.
! valu is used only for the error estimate.
    icount=0
    value=0.d0
    valu=0.d0
    do jj=nstack,nstack+2
      del1=xleft(jj+1)-xleft(jj)
!      print*, 'fn call 2'
      fval(1,jj)=fn(xleft(jj)+dx(1)*del1)
      fval(3,jj)=fn(xleft(jj)+dx(3)*del1)
!      print*, 'fn call 2'
      numcal = numcal + 2
      do kk=1,3
        icount=icount+1
        value=value+wt9(icount)*fval(kk,jj)*del
        valu=valu+fval(kk,jj)*wt(kk)*del1
      enddo
    enddo
    dif=abs(value-valu)
! If the following condition is true, add in this integral to the total,
! and reduce the number of regions under consideration.
    frac = del / (xmax - xmin)
    atsing = .false.
    if(frac .le. 1.0e-8) atsing = .true.
    if(dif .le. abr*frac .or. dif.le.rlr*abs(value) .or.           &
&       (atsing .and. (frac .le. 1.0e-15 .or. dif .le. abr*0.1  ))) then
      grater=grater+value
      error=error+abs(dif)
      nstack=nstack-1
! If no more regions, we are done.
      if(nstack.le.0) return
    else
! If the integration is insufficiently accurate, make each of the
! three subregions of the right-most region into regions.
! On next pass the right-most of these is the new current region.
      nstack=nstack+2
      maxns = max(maxns,nstack)
    endif
  enddo

end function grater

double complex function cgrater(fn,xmin,xmax,abr,rlr,nsing,xsing,error,numcal,maxns)
implicit none
integer :: nsing,numcal,maxns,mx,nstack,ii,jj,kk,icount
parameter (mx=1500)
double complex :: fn,value,valu,fval(3,mx)
double precision :: xmin,xmax,del,del1,dif,frac
double precision :: abr,rlr,error
double precision :: xleft(mx),dx(3),wt(3)
double precision :: wt9(9),xsing(20)
external fn
logical :: atsing
save dx,wt,wt9
data dx/0.1127016653792583,0.5,0.8872983346207417/
data wt/0.277777777777777778,0.4444444444444444444,0.2777777777777777778/
data wt9/0.0616938806304841571,0.108384229110206161,         &
&           0.0398463603260281088,0.175209035316976464,      &
&           0.229732989232610220,0.175209035316976464,       &
&           0.0398463603260281088,0.108384229110206161,      &
&           0.0616938806304841571  /

! nstack is the number of different intervals into which the
! integration region is currently divided. The number of regions can
! grow if more accuracy is needed by dividing the right-most region
! into three regions. The number of regions shrinks when the integral
! over the right-most region is accurate enough, in which case that
! integral is added to the total (stored in cgrater) and the region
! is removed from consideration (and a new region is the right-most).
  nstack=nsing+1
  maxns = nstack
  error=0.d0
  cgrater=(0.d0,0.d0)
! The array xleft stores the boundary points of the regions.
! The singular points just govern the initial placement of the regions.
  xleft(1)=xmin
  xleft(nsing+2)=xmax
  if(nsing.gt.0) then
    do jj=1,nsing
      xleft(jj+1)=xsing(jj)
    enddo
  endif
! For each region, calculate the function and store at three selected points.
  do ii=1,nstack
    del=xleft(ii+1)-xleft(ii)
    do jj=1,3
      fval(jj,ii)=fn(xleft(ii)+del*dx(jj))
    enddo
  enddo
  numcal = nstack * 3
  do
    if(nstack+3.ge.mx) then
      write(*,*) 'TOO MANY REGIONS'
      stop 0006
    endif
! Divide the rightmost region into three subregions.
    del=xleft(nstack+1)-xleft(nstack)
    xleft(nstack+3)=xleft(nstack+1)
    xleft(nstack+1)=xleft(nstack)+del*dx(1)*2.
    xleft(nstack+2)=xleft(nstack+3)-del*dx(1)*2.
! The three data points already found for the region become the
! middle data points (number 2 in first index of fval) for each region.
    fval(2,nstack+2)=fval(3,nstack)
    fval(2,nstack+1)=fval(2,nstack)
    fval(2,nstack)=fval(1,nstack)
! Now do the integral over the right-most region in two different ways-
! a three point integral (valu) over each of the three subregions
! and a more accurate nine-point integral (value) over whole region.
! valu is used only for the error estimate.
    icount=0
    value=(0.d0,0.d0)
    valu=(0.d0,0.d0)
    do jj=nstack,nstack+2
      del1=xleft(jj+1)-xleft(jj)
!      print*, 'fn call 2'
      fval(1,jj)=fn(xleft(jj)+dx(1)*del1)
      fval(3,jj)=fn(xleft(jj)+dx(3)*del1)
!      print*, 'fn call 2'
      numcal = numcal + 2
      do kk=1,3
        icount=icount+1
        value=value+wt9(icount)*fval(kk,jj)*del
        valu=valu+fval(kk,jj)*wt(kk)*del1
      enddo
    enddo
    dif=abs(value-valu)
! If the following condition is true, add in this integral to the total,
! and reduce the number of regions under consideration.
    frac = del / (xmax - xmin)
    atsing = .false.
    if(frac .le. 1.0e-8) atsing = .true.
    if(dif .le. abr*frac .or. dif.le.rlr*abs(value) .or.           &
&       (atsing .and. (frac .le. 1.0e-15 .or. dif .le. abr*0.1  ))) then
      cgrater=cgrater+value
      error=error+abs(dif)
      nstack=nstack-1
! If no more regions, we are done.
      if(nstack.le.0) return
    else
! If the integration is insufficiently accurate, make each of the
! three subregions of the right-most region into regions.
! On next pass the right-most of these is the new current region.
      nstack=nstack+2
      maxns = max(maxns,nstack)
    endif
  enddo

end function cgrater

! Copy of grater for two-dimensional integration
double precision function grater2(fn,xmin,xmax,abr,rlr,nsing,xsing,error,numcal,maxns)
implicit none
integer :: nsing,numcal,maxns,mx,nstack,ii,jj,kk,icount
parameter (mx=1500)
double precision :: fn,value,valu,fval(3,mx)
double precision :: xmin,xmax,del,del1,dif,frac
double precision :: abr,rlr,error
double precision :: xleft(mx),dx(3),wt(3)
double precision :: wt9(9),xsing(20)
external fn
logical :: atsing
save dx,wt,wt9
data dx/0.1127016653792583,0.5,0.8872983346207417/
data wt/0.277777777777777778,0.4444444444444444444,0.2777777777777777778/
data wt9/0.0616938806304841571,0.108384229110206161,         &
&           0.0398463603260281088,0.175209035316976464,      &
&           0.229732989232610220,0.175209035316976464,       &
&           0.0398463603260281088,0.108384229110206161,      &
&           0.0616938806304841571  /

! nstack is the number of different intervals into which the
! integration region is currently divided. The number of regions can
! grow if more accuracy is needed by dividing the right-most region
! into three regions. The number of regions shrinks when the integral
! over the right-most region is accurate enough, in which case that
! integral is added to the total (stored in grater2) and the region
! is removed from consideration (and a new region is the right-most).
  nstack=nsing+1
  maxns = nstack
  error=0.d0
  grater2=0.d0
! The array xleft stores the boundary points of the regions.
! The singular points just govern the initial placement of the regions.
  xleft(1)=xmin
  xleft(nsing+2)=xmax
  if(nsing.gt.0) then
    do jj=1,nsing
      xleft(jj+1)=xsing(jj)
    enddo
  endif
! For each region, calculate the function and store at three selected points.
  do ii=1,nstack
    del=xleft(ii+1)-xleft(ii)
    do jj=1,3
      fval(jj,ii)=fn(xleft(ii)+del*dx(jj))
    enddo
  enddo
  numcal = nstack * 3
  do
    if(nstack+3.ge.mx) then
      write(*,*) 'TOO MANY REGIONS'
      stop 0006
    endif
! Divide the rightmost region into three subregions.
    del=xleft(nstack+1)-xleft(nstack)
    xleft(nstack+3)=xleft(nstack+1)
    xleft(nstack+1)=xleft(nstack)+del*dx(1)*2.
    xleft(nstack+2)=xleft(nstack+3)-del*dx(1)*2.
! The three data points already found for the region become the
! middle data points (number 2 in first index of fval) for each region.
    fval(2,nstack+2)=fval(3,nstack)
    fval(2,nstack+1)=fval(2,nstack)
    fval(2,nstack)=fval(1,nstack)
! Now do the integral over the right-most region in two different ways-
! a three point integral (valu) over each of the three subregions
! and a more accurate nine-point integral (value) over whole region.
! valu is used only for the error estimate.
    icount=0
    value=0.d0
    valu=0.d0
    do jj=nstack,nstack+2
      del1=xleft(jj+1)-xleft(jj)
!      print*, 'fn call 2'
      fval(1,jj)=fn(xleft(jj)+dx(1)*del1)
      fval(3,jj)=fn(xleft(jj)+dx(3)*del1)
!      print*, 'fn call 2'
      numcal = numcal + 2
      do kk=1,3
        icount=icount+1
        value=value+wt9(icount)*fval(kk,jj)*del
        valu=valu+fval(kk,jj)*wt(kk)*del1
      enddo
    enddo
    dif=abs(value-valu)
! If the following condition is true, add in this integral to the total,
! and reduce the number of regions under consideration.
    frac = del / (xmax - xmin)
    atsing = .false.
    if(frac .le. 1.0e-8) atsing = .true.
    if(dif .le. abr*frac .or. dif.le.rlr*abs(value) .or.           &
&       (atsing .and. (frac .le. 1.0e-15 .or. dif .le. abr*0.1  ))) then
      grater2=grater2+value
      error=error+abs(dif)
      nstack=nstack-1
! If no more regions, we are done.
      if(nstack.le.0) return
    else
! If the integration is insufficiently accurate, make each of the
! three subregions of the right-most region into regions.
! On next pass the right-most of these is the new current region.
      nstack=nstack+2
      maxns = max(maxns,nstack)
    endif
  enddo

end function grater2

! Copy of grater for three-dimensional integration
double precision function grater3(fn,xmin,xmax,abr,rlr,nsing,xsing,error,numcal,maxns)
implicit none
integer :: nsing,numcal,maxns,mx,nstack,ii,jj,kk,icount
parameter (mx=1500)
double precision :: fn,value,valu,fval(3,mx)
double precision :: xmin,xmax,del,del1,dif,frac
double precision :: abr,rlr,error
double precision :: xleft(mx),dx(3),wt(3)
double precision :: wt9(9),xsing(20)
external fn
logical :: atsing
save dx,wt,wt9
data dx/0.1127016653792583,0.5,0.8872983346207417/
data wt/0.277777777777777778,0.4444444444444444444,0.2777777777777777778/
data wt9/0.0616938806304841571,0.108384229110206161,         &
&           0.0398463603260281088,0.175209035316976464,      &
&           0.229732989232610220,0.175209035316976464,       &
&           0.0398463603260281088,0.108384229110206161,      &
&           0.0616938806304841571  /

! nstack is the number of different intervals into which the
! integration region is currently divided. The number of regions can
! grow if more accuracy is needed by dividing the right-most region
! into three regions. The number of regions shrinks when the integral
! over the right-most region is accurate enough, in which case that
! integral is added to the total (stored in grater3) and the region
! is removed from consideration (and a new region is the right-most).
  nstack=nsing+1
  maxns = nstack
  error=0.d0
  grater3=0.d0
! The array xleft stores the boundary points of the regions.
! The singular points just govern the initial placement of the regions.
  xleft(1)=xmin
  xleft(nsing+2)=xmax
  if(nsing.gt.0) then
    do jj=1,nsing
      xleft(jj+1)=xsing(jj)
    enddo
  endif
! For each region, calculate the function and store at three selected points.
  do ii=1,nstack
    del=xleft(ii+1)-xleft(ii)
    do jj=1,3
      fval(jj,ii)=fn(xleft(ii)+del*dx(jj))
    enddo
  enddo
  numcal = nstack * 3
  do
    if(nstack+3.ge.mx) then
      write(*,*) 'TOO MANY REGIONS'
      stop 0006
    endif
! Divide the rightmost region into three subregions.
    del=xleft(nstack+1)-xleft(nstack)
    xleft(nstack+3)=xleft(nstack+1)
    xleft(nstack+1)=xleft(nstack)+del*dx(1)*2.
    xleft(nstack+2)=xleft(nstack+3)-del*dx(1)*2.
! The three data points already found for the region become the
! middle data points (number 2 in first index of fval) for each region.
    fval(2,nstack+2)=fval(3,nstack)
    fval(2,nstack+1)=fval(2,nstack)
    fval(2,nstack)=fval(1,nstack)
! Now do the integral over the right-most region in two different ways-
! a three point integral (valu) over each of the three subregions
! and a more accurate nine-point integral (value) over whole region.
! valu is used only for the error estimate.
    icount=0
    value=0.d0
    valu=0.d0
    do jj=nstack,nstack+2
      del1=xleft(jj+1)-xleft(jj)
!      print*, 'fn call 2'
      fval(1,jj)=fn(xleft(jj)+dx(1)*del1)
      fval(3,jj)=fn(xleft(jj)+dx(3)*del1)
!      print*, 'fn call 2'
      numcal = numcal + 2
      do kk=1,3
        icount=icount+1
        value=value+wt9(icount)*fval(kk,jj)*del
        valu=valu+fval(kk,jj)*wt(kk)*del1
      enddo
    enddo
    dif=abs(value-valu)
! If the following condition is true, add in this integral to the total,
! and reduce the number of regions under consideration.
    frac = del / (xmax - xmin)
    atsing = .false.
    if(frac .le. 1.0e-8) atsing = .true.
    if(dif .le. abr*frac .or. dif.le.rlr*abs(value) .or.           &
&       (atsing .and. (frac .le. 1.0e-15 .or. dif .le. abr*0.1  ))) then
      grater3=grater3+value
      error=error+abs(dif)
      nstack=nstack-1
! If no more regions, we are done.
      if(nstack.le.0) return
    else
! If the integration is insufficiently accurate, make each of the
! three subregions of the right-most region into regions.
! On next pass the right-most of these is the new current region.
      nstack=nstack+2
      maxns = max(maxns,nstack)
    endif
  enddo

end function grater3

max_stars_repo_name: dinisAbranches/nwchem

max_issues_repo_name: dinisAbranches/nwchem

max_forks_repo_name: dinisAbranches/nwchem


text:
##||<bgcolor#E0E0FF>Menus||
##||https://themelt.com/menu Menu||
## remove the ## characters if adding a menu

The Melt is a grilled cheese and soup Restaurants restaurant (as http://www.davisenterprise.com/business/newgrilledcheese/ noted by Wendy Wetzel). It fills the spot previously known as Ben & Jerrys.

The Melt has a simple menu, reminiscent of InNOut Burger. The food is prepared behind a protective glass screen in front of customers. When ordering at the futuremodern looking white cash registers, they ask for your initials which are displayed alongside your order on a large screen in the upper right hand corner of the kitchen so you can see where you are in the order queue. If you order online, you pick up your order in store by scanning your digital receipt at the built in scanner to the right of the cash registers. Ordering online allows you to skip any lines and your order goes straight to the top of the order queue.

Students will always get a free soda by showing a student ID and purchasing a combo (grilled cheese + side of soup or salad). The soda is made with cane sugar rather than http://en.wikipedia.org/wiki/High_fructose_corn_syrup HFCS and comes in only three flavors: cola, rootbeer, & lemonlime.

They have a fenced off exterior patio, opposite the path of Mikuni Mikunis patio, which more than doubles their available seating.

I got an informal confirmation today that they are continuing their training/softopening tomorrow (Nov. 8) for the dinner hour (5pm8pm), which probably means free sandwiches again. Users/jefftolentino
They gave away free grilled cheese on November 12th, a day before their grand opening.

20131101 12:47:27 nbsp Well, they have open wifi without being open for business.  I approve.  (So far, anyways)Users/DavidBarnum

20131106 13:17:25 nbsp Happened to catch them on thier staff training day, when they were giving away everything for free.  Got an Italian Job with onions and bacon, with a green salad and a soda pop.  The Italian was made with garlic bread, and was highly satisfing.  The green salad was equally good, and included apples and goat cheese.  The fact that it was free was a great incentive to try this place, but I will probably go back after they officially open and pay actually money for more.  The menu shows sandwiches for about $5, with a combo option for an additional $2.   Users/jefftolentino

20131107 14:05:42 nbsp Came back today and got another free lunch.  Got the Mac Daddy, which was basically a grilled cheese where the cheese was actually macaroni and cheese.  Ill say that again.  I got a grilled cheese were the cheese was actually macaroni and cheese.  Yeah, basically the two most kickass comfort foods in one and it was fricking good.  I got the smoked bacon added in to make it extra special.  It was good and Id get it again.  On the down side they overdid it a little on the seasoning that they sprinkle on the outside of the bread.  It was a lot, like more than reasonable.  Plus it was a processed seasoning mix, and not super impressive. But it was free, so I wasnt too worried about it.  I got a cesaer salad on the side that was typical fair.  They had croutons that look like they were made from grilled cheese cubes which were visually interesting, but not much different from normal croutons tastewise.  Personlly I preferred the green salad with goat cheese and apples that I got yesterday.  Moving on, I also ordered chips since everything was free, but they forgot them.  It was also crowded, and again free, so Im not going to complain too much.  If it were something Id paid for, I probably would have just asked for my chips and been done with it.  Given all that, Im rating these guys an 8/10.  Its not 100% prime ingredients, but they do have some really good grilled cheese sandwiches and the Mac Daddy was pretty awesome. Users/jefftolentino

20131209 14:13:07 nbsp Had lunch here today, and it was delicious! The Mission grilled cheese with tomato, spicy black bean soup, and a pickle for $6.95a great price for a filling meal (+free beverage for UCD undergrads, grads, and staff). Definitely going back. Users/zombiek

20140327 21:53:16 nbsp I went a few weeks ago when I was just getting over a cold, so Id be open to try them again sometime. I had their jalapeno melt. It was pretty good. The italian soda was nice as well. I hope they stick around. They fit a unique niche in the dinning options in Davis. Users/CarlosOverstreet