text:
C Subroutine EXSECT Calculates electron impact cross sections
C
C This software is part of the GLOW model.  Use is governed by the Open Source
C Academic Research License Agreement contained in the file glowlicense.txt.
C For more information see the file glow.txt.
C
C Adapted from Banks & Nagy 2-stream code by Stan Solomon, 1988
C Added high-energy relativistic cross section correction, SCS, 1999
C Updated comments, SCS, 2002
C Included in GLOW v. 0.97, SCS, 2005
C Replaced common blocks with use-associated variables defined in module cglow, Ben Foster, 2015
C Reduced O(3p5P) (7774) cross section A0 from 0.817 to 0.327, SCS, 1/2017
C
C Definitions:
C SIGS   elastic cross sections for each species, energy; cm2
C PE     elastic backscatter probabilities for each species, energy
C PIN    inelastic  "
C SIGA   energy loss cross section for each species, loss, energy; cm2
C SEC    secondary production xsect for species, Esec, Epri; cm2
C SIGEX  excitation xsect for each state, species, energy; cm2
C     O states:   1D,   1S, 3s5S, 3s3S, 3p5P, 3p3P, 3d3D, 3s'3D
C     O2 states:   a,    b, AA'c,    B,  9.9, Ryds,  vib
C     N2 states: ABW,   B',    C, aa'w,  1Pu,   b', Ryds,  vib
C SIGIX  ionization xsect for each state, species, energy; cm2
C     O states:   4S,  2Do,  2Po
C     O2 states:   X,    a,    A,    b,    B,   c,  37eV
C     N2 states:   X,    A,    B,    D,    C, 40eV
C IIMAX  number of bins for secondary production for each primary energy
C WW     energy threshold for each excited state, species; eV
C WW, AO, OMEG, ANU, BB: revised excitation cross section parameters,
C        from Green & Stolarski (1972) formula (W, A, omega, nu, gamma)
C AUTO   autoionization coefs (= 0 as autoion. included in ion xsects)
C THI    energy threshold for each ionized state, species; eV
C AK, AJ, TS, TA, TB, GAMS, GAMB:  Jackman et al (1977) ioniz. params
C ENER   energy grid; eV
C DEL    energy grid spacing; eV
C NNN    number of excited states for each species
C NINN   number of ionized states for each species
C NUM    number of points on elastic data trid for each species
C EC     data energy grid of elastic xsects and backscatter ratios
C        for each species; eV
C CC     elastic xsects on data grid for each species, cm2
C CE     elastic backscat. probs on data grid for each species; cm2
C CI     inelastic "
C
C Array dimensions:
C NBINS  number of energy levels
C NMAJ   number of major species
C NEI    number of slots for excited and ionized states
C
C
      SUBROUTINE EXSECT (ENER, DEL)
C
      use cglow,only: NMAJ,NEI,NBINS
      use cglow,only: WW,AO,OMEG,ANU,BB,AUTO,THI,AK,AJ,TS,TA,TB,GAMS, ! /CXPARS/
     |                GAMB
      use cglow,only: SIGS,PE,PIN,SIGEX,SIGIX,SIGA,SEC,SIGA,IIMAXX    ! /CXSECT/
C
      implicit none
C
      real,intent(in) :: ENER(NBINS), DEL(NBINS)
C
      real ::   SIGI(NBINS), T12(NBINS), RATIO(NBINS),
     |          EC(31,NMAJ), CC(31,NMAJ), CE(31,NMAJ), CI(31,NMAJ)
      integer :: NNNI(NMAJ),NINN(NMAJ),NUM(NMAJ),NNN(NMAJ)
      integer :: ij,iv,ii,i,k,j,i1,i2,i3,jy,kuk,kuk1,iee,kk,ie,ibz,ml,
     |  itmax
      real :: EX,FAC,WE,AE,GAMMA,T0,ETJ,DETJ,SIGG,ETA,FF,WAG,TMT,
     |  E1,E2,TMAX,WTH1
      integer,external :: INV
      real,external :: SIGION
C
      real,parameter :: QQN = 6.51E-14
C
      NNN  = (/8,7,8/)
      NINN = (/3,7,6/)
      NUM  = (/31,28,28/)
C
C Set arrays in module cglow (formerly COMMON/CXSECT/)
C (all are dimensioned (nei,nmaj) (10,3))
C
      WW(1:nei,1)=
     |  (/1.96, 4.17, 9.29, 9.53,10.76,10.97,12.07,12.54, 0., 0./)
      WW(1:nei,2)=
     |  (/0.98, 1.64, 4.50, 8.44, 9.90,13.50, 0.25, 0.00, 0., 0./)
      WW(1:nei,3)=
     |  (/6.17, 8.16,11.03, 8.40,12.85,14.00,13.75, 1.85, 0., 0./)

      AO(1:nei,1)=
     |  (/.0100,.0042,.1793,.3565,.0327,.0245,.0293,.1221, 0.,0./)
      AO(1:nei,2)=
     |  (/.0797,.0211,.0215,.3400,.0657,1.110,3.480, 0.00, 0.,0./)
      AO(1:nei,3)=
     |  (/2.770,.1140,.1790,.0999,.8760,.6010,1.890,1.350, 0.,0./)

      OMEG(1:nei,1)=
     |  (/1.00, 1.00, 3.00, 0.75, 3.00, 0.85, 0.75, 0.75, 0.,0./)
      OMEG(1:nei,2)=
     |  (/2.00, 2.00, 1.15, 0.75, 0.75, 0.75, 7.00, 0.00, 0.,0./)
      OMEG(1:nei,3)=
     |  (/3.00, 3.00, 3.00, 1.00, 0.75, 0.75, 0.75, 8.00, 0.,0./)

      ANU(1:nei,1)=
     |  (/2.00, 1.04, 2.53, 0.54, 2.43, 2.87, 0.93, 0.72, 0.,0./)
      ANU(1:nei,2)=
     |  (/6.18, 4.14, 1.00, 1.05, 1.60, 3.00,10.87, 0.00, 0.,0./)
      ANU(1:nei,3)=
     |  (/4.53, 4.78, 4.32, 4.05, 1.47, 1.27, 3.00, 1.58, 0.,0./)

      BB(1:nei,1)=
     |  (/1.00, 0.50, 1.02, 0.01, 4.19, 4.88, 0.66, 0.17, 0.,0./)
      BB(1:nei,2)=
     |  (/0.53, 0.51, 0.98, 0.99, 1.86, 1.00, 1.00, 0.00, 0.,0./)
      BB(1:nei,3)=
     |  (/1.42, 3.54,12.70, 5.20, 0.86, 0.45, 1.00, 1.00, 0.,0./)

      AUTO(1:nei,1) = (/0.,0.,0.,0.,0.,0.,0.,0.,0.,0./)
      AUTO(1:nei,2) = (/0.,0.,0.,0.,0.,0.,0.,0.,0.,0./)
      AUTO(1:nei,3) = (/0.,0.,0.,0.,0.,0.,0.,0.,0.,0./)

      THI(1:nei,1)=
     |  (/13.60,16.90,18.50, 0.00, 0.00, 0.00, 0.00, 0.,0.,0./)
      THI(1:nei,2)=
     |  (/12.10,16.10,16.90,18.20,20.00,23.00,37.00, 0.,0.,0./)
      THI(1:nei,3)=
     |  (/15.58,16.73,18.75,22.00,23.60,40.00, 0.00, 0.,0.,0./)

      AK(1:nei,1)=
     |  (/ 1.13, 1.25, 0.67, 0.00, 0.00, 0.00, 0.00, 0.,0.,0./)
      AK(1:nei,2)=
     |  (/0.47, 1.13, 1.13, 1.01, 0.65, 0.95, 0.59, 0.,0.,0./)
      AK(1:nei,3)=
     |  (/2.42, 1.06, 0.55, 0.37, 0.37, 0.53, 0.00, 0.,0.,0./)

      AJ(1:nei,1)=
     |  (/1.81, 1.79, 1.78, 0.00, 0.00, 0.00, 0.00, 0.,0.,0./)
      AJ(1:nei,2)=
     |  (/3.76, 3.76, 3.76, 3.76, 3.76, 3.76, 3.76, 0.,0.,0./)
      AJ(1:nei,3)=
     |  (/1.74, 1.74, 1.74, 1.74, 1.74, 1.74, 0.00, 0.,0.,0./)

      TS(1:nei,1)=
     |  (/6.41, 6.41, 6.41, 0.00, 0.00, 0.00, 0.00, 0.,0.,0./)
      TS(1:nei,2)=
     |  (/1.86, 1.86, 1.86, 1.86, 1.86, 1.86, 1.86, 0.,0.,0./)
      TS(1:nei,3)=
     |  (/4.71, 4.71, 4.71, 4.71, 4.71, 4.71, 0.00, 0.,0.,0./)

      TA(1:nei,1)=
     |  (/3450.,3450.,3450.,   0.,   0.,   0.,   0., 0.,0.,0./)
      TA(1:nei,2)=
     |  (/1000.,1000.,1000.,1000.,1000.,1000.,1000., 0.,0.,0./)
      TA(1:nei,3)=
     |  (/1000.,1000.,1000.,1000.,1000.,1000.,   0., 0.,0.,0./)

      TB(1:nei,1)=
     |  (/162.00,162.0,162.0, 0.00, 0.00, 0.00, 0.00, 0.,0.,0./)
      TB(1:nei,2)=
     |  (/24.20,32.20,33.80,36.40,40.60,46.00,74.00, 0.,0.,0./)
      TB(1:nei,3)=
     |  (/31.16,33.46,37.50,44.00,47.20,80.00, 0.00, 0.,0.,0./)

      GAMS(1:nei,1)=
     |  (/13.00,13.0,13.00, 0.00, 0.00, 0.00, 0.00, 0.,0.,0./)
      GAMS(1:nei,2)=
     |  (/18.50,18.5,18.50,18.50,18.50,18.50,18.50, 0.,0.,0./)
      GAMS(1:nei,3)=
     |  (/13.80,13.8,13.80,13.80,13.80,13.80, 0.00, 0.,0.,0./)

      GAMB(1:nei,1)=
     |  (/-.815,-.815,-.815,0.00, 0.00, 0.00, 0.00, 0.,0.,0./)
      GAMB(1:nei,2)=
     |  (/12.10,16.10,16.90,18.2,20.30,23.00,37.00, 0.,0.,0./)
      GAMB(1:nei,3)=
     |  (/15.58,16.73,18.75,22.0,23.60,40.00, 0.00, 0.,0.,0./)

      EC(1:31,1) = (/
     |               1.00,     2.00,     4.00,     6.00,     8.00,
     |              10.00,    12.00,    14.00,    16.00,    18.00,
     |              20.00,    30.00,    40.00,    50.00,    60.00,
     |              70.00,    80.00,    90.00,   100.00,   150.00,
     |             200.00,   300.00,   500.00,  1000.00,  2000.00,
     |            3000.00,  5000.00, 10000.00, 20000.00, 40000.00,
     |           50000.00/)

      EC(1:31,2) = (/
     |               1.00,     2.00,     3.00,     5.00,     7.00,
     |              10.00,    15.00,    20.00,    30.00,    40.00,
     |              50.00,    70.00,   100.00,   150.00,   200.00,
     |             300.00,   400.00,   500.00,   600.00,   700.00,
     |            1000.00,  2000.00,  3000.00,  5000.00, 10000.00,
     |           20000.00, 40000.00, 50000.00,     0.00,     0.00,
     |               0.00/)

      EC(1:31,3) = (/
     |               1.00,     2.00,     2.50,     3.00,     4.00,
     |               5.00,     6.00,     8.00,    10.00,    15.00,
     |              20.00,    30.00,    40.00,    50.00,    70.00,
     |             100.00,   200.00,   300.00,   500.00,   700.00,
     |            1000.00,  2000.00,  3000.00,  5000.00, 10000.00,
     |           20000.00, 40000.00, 50000.00,     0.00,     0.00, 0.0/)

      CC(1:31,1) = (/
     |           5.00E-16, 6.00E-16, 7.50E-16, 7.60E-16, 7.70E-16,
     |           7.80E-16, 7.50E-16, 7.20E-16, 6.90E-16, 6.70E-16,
     |           6.50E-16, 5.60E-16, 4.60E-16, 4.00E-16, 3.50E-16,
     |           3.20E-16, 2.90E-16, 2.70E-16, 2.50E-16, 1.90E-16,
     |           1.50E-16, 1.20E-16, 8.00E-17, 5.00E-17, 3.02E-17,
     |           1.99E-17, 1.20E-17, 6.08E-18, 3.06E-18, 1.55E-18,
     |           1.24E-18/)

      CC(1:31,2) = (/
     |           5.50E-16, 6.90E-16, 7.50E-16, 8.50E-16, 9.60E-16,
     |           1.00E-15, 1.00E-15, 9.00E-16, 8.30E-16, 7.70E-16,
     |           6.90E-16, 5.70E-16, 4.40E-16, 3.30E-16, 2.70E-16,
     |           2.10E-16, 1.80E-16, 1.60E-16, 1.40E-16, 1.30E-16,
     |           1.10E-16, 7.00E-17, 5.00E-17, 3.00E-17, 1.53E-17,
     |           7.72E-18, 3.90E-18, 3.13E-18, 0.00E+00, 0.00E+00,
     |           0.00E+00/)

      CC(1:31,3) = (/
     |           9.00E-16, 2.27E-15, 2.52E-15, 1.93E-15, 1.32E-15,
     |           1.15E-15, 1.16E-15, 1.17E-15, 1.18E-15, 1.14E-15,
     |           1.13E-15, 9.50E-16, 8.60E-16, 7.30E-16, 5.90E-16,
     |           4.70E-16, 3.30E-16, 2.50E-16, 1.60E-16, 1.30E-16,
     |           1.10E-16, 6.35E-17, 4.18E-17, 2.54E-17, 1.28E-17,
     |           6.44E-18, 3.27E-18, 2.62E-18, 0.00E+00, 0.00E+00, 0.0/)

      CE(1:31,1) = (/
     |            0.50000,  0.49500,  0.46800,  0.43600,  0.42000,
     |            0.40500,  0.37000,  0.36000,  0.34000,  0.33000,
     |            0.32000,  0.27000,  0.24000,  0.22000,  0.20000,
     |            0.18000,  0.17000,  0.16000,  0.15000,  0.13000,
     |            0.11500,  0.09000,  0.06800,  0.04600,  0.02400,
     |            0.01660,  0.01000,  0.00510,  0.00255,  0.00125,
     |            0.00100/)

      CE(1:31,2) = (/
     |            0.50000,  0.50000,  0.49000,  0.44500,  0.42700,
     |            0.40500,  0.36800,  0.34300,  0.31600,  0.28900,
     |            0.25800,  0.22000,  0.18400,  0.16400,  0.13300,
     |            0.11000,  0.10000,  0.09200,  0.08500,  0.08000,
     |            0.06800,  0.03700,  0.02600,  0.01600,  0.00800,
     |            0.00400,  0.00200,  0.00160,  0.00000,  0.00000,
     |            0.00000/)

      CE(1:31,3) = (/
     |            0.50000,  0.50000,  0.50000,  0.49000,  0.46800,
     |            0.44500,  0.43600,  0.42000,  0.40500,  0.36800,
     |            0.34300,  0.31600,  0.28900,  0.25800,  0.22000,
     |            0.18400,  0.14000,  0.11000,  0.08400,  0.07400,
     |            0.06300,  0.03400,  0.02400,  0.01500,  0.00740,
     |            0.00370,  0.00180,  0.00140,  0.00000,  0.00000, 0.0/)

      CI(1:31,1) = (/
     |            0.60000,  0.60000,  0.60000,  0.60000,  0.60000,
     |            0.60000,  0.55000,  0.46000,  0.40000,  0.36000,
     |            0.32000,  0.22000,  0.15000,  0.10000,  0.08200,
     |            0.07000,  0.06100,  0.05400,  0.05000,  0.04400,
     |            0.03800,  0.02800,  0.02000,  0.01050,  0.00600,
     |            0.00400,  0.00250,  0.00130,  0.00060,  0.00030,
     |            0.00025/)

      CI(1:31,2) = (/
     |            0.50000,  0.50000,  0.50000,  0.50000,  0.48000,
     |            0.44000,  0.36000,  0.28000,  0.20000,  0.14000,
     |            0.10000,  0.07000,  0.05000,  0.04600,  0.04300,
     |            0.03700,  0.03200,  0.02800,  0.02400,  0.02100,
     |            0.01600,  0.00900,  0.00620,  0.00400,  0.00200,
     |            0.00100,  0.00050,  0.00040,  0.00000,  0.00000,
     |            0.00000/)

      CI(1:31,3) = (/
     |            0.50000,  0.50000,  0.50000,  0.50000,  0.50000,
     |            0.50000,  0.50000,  0.50000,  0.50000,  0.50000,
     |            0.44000,  0.30000,  0.20000,  0.13000,  0.09000,
     |            0.06000,  0.05000,  0.04200,  0.03200,  0.02500,
     |            0.02000,  0.01100,  0.00800,  0.00500,  0.00250,
     |            0.00120,  0.00060,  0.00050,  0.00000,  0.00000, 0.0/)
C
C
C Interpolate elastic cross sections and backscatter ratios:
C
      DO 90 IJ=1,NMAJ
        DO  80 IV=1,NBINS
          EX=ENER(IV)
          DO 50 II=1,NUM(IJ)

            IF (EC(II,IJ) .GT. EX) GOTO 60
   50     CONTINUE
          SIGS(IJ,IV)=CC(NUM(IJ),IJ)*(EC(NUM(IJ),IJ)/EX)**0.8
          IF(IJ.EQ.1) SIGS(IJ,IV)=CC(NUM(IJ),IJ)*(EC(NUM(IJ),IJ)/EX)**2

          PE(IJ,IV) = CE(NUM(IJ),IJ)* (EC(NUM(IJ),IJ)/EX)
          PIN(IJ,IV) = CI(NUM(IJ),IJ)* (EC(NUM(IJ),IJ)/EX)
          GOTO 80
   60     I=II-1
          IF (I .LE. 0) THEN
            SIGS(IJ,IV)=CC(II,IJ)
            PE(IJ,IV)=CE(II,IJ)
            PIN(IJ,IV)=CI(II,IJ)
          ELSE
            FAC = ALOG (EX/EC(I,IJ)) / ALOG (EC(II,IJ)/EC(I,IJ))
            SIGS(IJ,IV) = EXP (ALOG (CC(I,IJ))
     >                         + ALOG (CC(II,IJ)/CC(I,IJ)) * FAC)
            PE(IJ,IV) = EXP (ALOG (CE(I,IJ))
     >                       + ALOG (CE(II,IJ)/CE(I,IJ)) * FAC)
            PIN(IJ,IV) = EXP (ALOG (CI(I,IJ))
     >                       + ALOG (CI(II,IJ)/CI(I,IJ)) * FAC)
          ENDIF
   80   CONTINUE
   90 CONTINUE
C
C
C Calculate electron impact excitation and ionization cross sections:
C
      DO 140 I=1,NMAJ
        DO 140 K=1,NEI
          DO 140 J=1,NBINS
          IF (ENER(J).GT.WW(K,I) .AND. WW(K,I).GT.0.001) THEN
            WE = WW(K,I) / ENER(J)
            SIGEX(K,I,J) = QQN * AO(K,I)
     >                       * (WE**OMEG(K,I) / WW(K,I)**2)
     >                       * (1.0 - WE**BB(K,I)) ** ANU(K,I)
            IF (SIGEX(K,I,J) .LT. 1.E-30) SIGEX(K,I,J) = 0.0
          ELSE
            SIGEX(K,I,J) = 0.0
          ENDIF
          IF (ENER(J).GT.THI(K,I) .AND. THI(K,I).GT.0.001) THEN
            AE = AK(K,I)/ENER(J) * ALOG(ENER(J)/AJ(K,I))
            GAMMA = GAMS(K,I) * ENER(J) / (ENER(J)+GAMB(K,I))
            T0 = TS(K,I) - (TA(K,I)/(ENER(J)+TB(K,I)))
            SIGIX(K,I,J) = 1.E-16 * AE * GAMMA
     >                    * ( ATAN(((ENER(J)-THI(K,I))/2.-T0)/GAMMA)
     >                       +ATAN(T0/GAMMA) )
            IF (SIGIX(K,I,J) .LT. 1.E-30) SIGIX(K,I,J) = 0.0
          ELSE
            SIGIX(K,I,J) = 0.0
          ENDIF
  140 CONTINUE
C
C
C Obtain high-energy correction factors:
C
      CALL HEXC(ENER,SIGIX,RATIO)
      DO 142 J=1,NBINS
        DO 142 I=1,NMAJ
          DO 142 K=1,NEI
            SIGIX(K,I,J)=SIGIX(K,I,J)/RATIO(J)
  142 CONTINUE
C
C
C Zero energy loss xsect and secondary production xsect arrays:
C
      DO 145 I1=1,NMAJ
      DO 145 I2=1,NBINS
      DO 145 I3=1,NBINS
        SIGA(I1,I2,I3)=0.0
        SEC(I1,I2,I3)=0.0
  145 CONTINUE
C
C
C Loop over energy:
C
      DO 500 JY=1,NBINS
C
      KUK=0
      KUK1=0
      ETJ=ENER(JY)
      DETJ=DEL(JY)
C
C
C Loop over species:
C
      DO 400 I=1,NMAJ
C
C
C Loop over excited states:
C
      DO 200 J=1,NNN(I)
C
C
C Calculate energy loss from JY to J-K for each species.
C The cross secton is divided proportionally between bin INV and bin
C INV-1, the two bins closest to J-K:
C
      SIGG = SIGEX(J,I,JY)
      ETA = ETJ - WW(J,I)
      IF (ETA .GT. 0.) THEN
        IE = INV (ETA,JY,ENER)
        IEE = IE - 1
        IF (IEE .LT. 1) IEE=IE
        K = JY - IE
        KK = JY - IEE
        IF (KK .GE. KUK) KUK = KK
        IF (IE .EQ. JY) THEN
          IF (JY .EQ. 1) THEN
            SIGA(I,1,JY) = SIGA(I,1,JY) + SIGG
          ELSE
            SIGA(I,1,JY) = SIGA(I,1,JY) + SIGG * (DETJ/DEL(JY-1))
     >                    * WW(J,I) / (ENER(JY)-ENER(JY-1))
          ENDIF
        ELSE
          IF (IE .EQ. 1) THEN
C           SIGA(I,K,JY) = SIGA(I,K,JY) + SIGG * 2.*ETA*DETJ/DEL(1)**2
            SIGA(I,K,JY) = SIGA(I,K,JY) + SIGG * DETJ/DEL(1)
          ELSE
            FF = (ENER(IE)-ETA) / (ENER(IE)-ENER(IEE))
            FF = 1.0 - ABS(FF)
            SIGA(I,K,JY) = SIGA(I,K,JY) + SIGG * FF * DETJ/DEL(IE)
            SIGA(I,KK,JY)= SIGA(I,KK,JY) + SIGG * (1.0-FF)*DETJ/DEL(IEE)
          ENDIF
          WAG=WW(J,I)-THI(1,I)
          IF (WAG .GT. 0. .AND. AUTO(J,I) .GT. 0.) THEN
            IBZ = INV (WAG,JY,ENER)
            SEC(I,IBZ,JY) = SEC(I,IBZ,JY)+SIGG*(DETJ/DEL(IBZ))*AUTO(J,I)
            IF(IBZ.GE.KUK1)KUK1=IBZ
          ENDIF
        ENDIF
      ENDIF
  200 CONTINUE
C
C
C Loop over ion states:
C
      DO 300 ML=1,NINN(I)
C
      DO 210 II=1,NBINS
      SIGI(II) = 0.0
      T12(II) = 0.0
  210 CONTINUE
C
C
C Calculate cross-section for production of secondaries into each
C bin from 1 to ITMAX and store in SIGI(II).  Apply relativistic correction.
C Also store the average energy of the secondaries in T12(II):
C
      WAG = THI(ML,I)
      TMAX = (ETJ-WAG) / 2.
      if (tmax .gt. 1.e6) tmax=1.e6
      IF (TMAX .LE. 0.) GOTO 300
      ITMAX = INV (TMAX,JY,ENER)
      IF (ITMAX .GE. KUK1) KUK1 = ITMAX + 1
      TMT = ENER(1) + DEL(1) / 2.0
      IF (TMAX .LT. TMT) TMT = TMAX
      SIGI(1) = SIGION(I,ML,ETJ,0.0,TMT,T12(1)) / RATIO(JY)
      TMT = ENER(1) + DEL(1) / 2.
      IF (TMAX .GT. TMT) THEN
        IF (TMAX .LE. ENER(2)) ITMAX = 2
        DO 220 II=2,ITMAX
        E1 = ENER(II) - DEL(II) / 2.
        E2 = E1 + DEL(II)
        IF (E2 .GT. TMAX) E2 = TMAX
        IF(E1 .LE. E2)SIGI(II)=SIGION(I,ML,ETJ,E1,E2,T12(II))/RATIO(JY)
  220   CONTINUE
      ENDIF
C
C
C Add the secondary production cross-section to SEC; calculate
C the ionization energy loss cross-section and add to SIGA:
C
      DO 250 II=1,ITMAX
      SEC(I,II,JY) = SEC(I,II,JY) + SIGI(II) * DETJ / DEL(II)
      WTH1 = T12(II) + WAG
      ETA = ETJ - WTH1
      IF (ETA .GT. 0.) THEN
        IE = INV (ETA,JY,ENER)
        IEE = IE - 1
        K = JY - IE
        KK = JY - IEE
        IF (IEE .LT. 1) IEE = IE
        IF (IE .EQ. JY) THEN
          SIGA(I,1,JY) = SIGA(I,1,JY) + SIGI(II) * (DETJ/DEL(JY-1))
     >                  * WTH1 / (ENER(JY)-ENER(JY-1))
        ELSE
          IF (KK .GE. KUK) KUK = KK
          IF (IE .EQ. 1) THEN
C           SIGA(I,K,JY) = SIGA(I,K,JY)+2.*ETA*SIGI(II)*DETJ/DEL(1)**2
            SIGA(I,K,JY) = SIGA(I,K,JY)+SIGI(II)*DETJ/DEL(1)
          ELSE
            FF = (ENER(IE)-ETA) / (ENER(IE)-ENER(IEE))
            FF = 1.0 - ABS(FF)
            SIGA(I,K,JY) = SIGA(I,K,JY) + FF*SIGI(II)*DETJ/DEL(IE)
            SIGA(I,KK,JY)=SIGA(I,KK,JY)+(1.0-FF)*SIGI(II)*DETJ/DEL(IEE)
          ENDIF
        ENDIF
      ENDIF
C
  250 CONTINUE
C
  300 CONTINUE
C
  400 CONTINUE
C
      IIMAXX(JY) = KUK1
C
  500 CONTINUE
C
      RETURN

      END SUBROUTINE EXSECT
C
!-----------------------------------------------------------------------
C
C Function SIGION calculates ionization cross section for species I,
C state ML, primary energy E, secondary energy from E1 to E2
C
      real FUNCTION SIGION(I,ML,E,E1,E2,T12)
      use cglow,only: THI,AK,AJ,TS,TA,TB,GAMS,GAMB ! /CXPARS/
      implicit none
C
      integer,intent(in) :: i,ml
      real,intent(in) :: E,E1
      real,intent(out) :: T12
      real,intent(inout) :: E2
C
      DOUBLE PRECISION ABB, ABC, ABD
      real :: QQ
      DATA QQ/1.E-16/
      real :: AK1,AJ1,TS1,TA1,TB1,GAMS1,GAMB1,S,A,TZ,GG,TTL,AL2,AL1,
     |  TTL1
C
C
      IF (E .LE. THI(ML,I)) GOTO 30
C
      AK1=AK(ML,I)
      AJ1=AJ(ML,I)
      TS1=TS(ML,I)
      TA1=TA(ML,I)
      TB1=TB(ML,I)
      GAMS1=GAMS(ML,I)
      GAMB1=GAMB(ML,I)
      S=QQ*AK1*ALOG(E/AJ1)
      A=S/E
      TZ=TS1-TA1/(E+TB1)
      GG=(GAMS1*E)/(E+GAMB1)
      TTL=(E-THI(ML,I))/2.0
      TTL1=TTL-0.01
      IF(E1.GE.TTL1)GO TO 30
      IF(E2.GE.TTL)E2=TTL
      ABB=(E2-TZ)/GG
      ABC=(E1-TZ)/GG
      AL2=GG*GG*(ABB*ABB+1.0)
      AL1=GG*GG*(ABC*ABC+1.0)
      ABD=DATAN(ABB)-DATAN(ABC)
      T12=TZ+0.5*GG*(ALOG(AL2)-ALOG(AL1))/ABD
      SIGION=A*GG*ABD
      RETURN
C
   30 SIGION=0.0
      RETURN
C
      END FUNCTION SIGION
C
!-----------------------------------------------------------------------
C
C Function INV finds the bin number closest to energy ETA on grid ENER.
C Bin INV or INV-1 will contain ETA.
C
      INTEGER FUNCTION INV (ETA, JY, ENER)
      use cglow,only: NBINS
      implicit none
C
      real,intent(in) :: ETA,ENER(NBINS)
      integer,intent(in) :: JY
C
      integer :: iv
C
      IF (ETA .LT. 0.) THEN
        INV = -1
      ELSE
        DO 30 IV=1,JY
          IF (ETA .LE. ENER(IV)) GOTO 40
   30   CONTINUE
        IV = JY
   40   INV = IV
      ENDIF
C
      RETURN

      END FUNCTION INV
C
!-----------------------------------------------------------------------
C
C Subroutine HEXC
C
C High Energy Cross Section Correction
C Calculates ratio of low energy (non-relativistic) to high energy
C (relativistic) ionization cross sections, based on N2.
C Extends to 1 GeV.
C
C Originally coded by Ann Windnagel, 11/98
C Re-written by Stan Solomon, 2/99
C Re-designed with table lookup, SCS, 4/99
C Updated comments, SCS, 4/02
C References:
C   Porter et al., J. Chem. Phys., 65, 154, 1976.
C   Rieke and Prepejchal, Phys. Rev. A, 6, 1507, 1990.
C   Saksena et al., Int. Jour. of Mass Spec. & Ion Proc., 171, L1, 1997.

      SUBROUTINE HEXC(ENER,SIGIX,RATIO)
      use cglow,only: NEI,NMAJ,NBINS
      implicit none
C
      real,intent(in) :: ENER(NBINS),SIGIX(NEI,NMAJ,NBINS)
      real,intent(out) :: RATIO(NBINS)
C
      real ::  TOTX(NBINS), TOTNEW(NBINS), EGR(13), SGR(13)
      DATA EGR/1.E4,      2.E4,      5.E4,      1.E5,      2.E5,
     >         3.E5,      5.E5,      1.E6,      2.E6,      5.E6,
     >         1.E7,      1.E8,      1.E9/
      DATA SGR/1.20E-17,  7.03E-18,  3.37E-18,  1.96E-18,  1.26E-18,
     >         1.05E-18,  9.50E-19,  9.00E-19,  9.00E-19,  9.40E-19,
     >         1.00E-18,  1.26E-18,  1.59E-18/
      integer :: k,i,kg
      real,external :: TERPOO

C Calculate total low-energy cross section for N2:

      DO 20 K = 1,NBINS
        TOTX(K) = 0.
        DO 20 I = 1,NEI
          TOTX(K) = TOTX(K) + SIGIX(I,3,K)
   20 CONTINUE

C Calculate high-energy cross section for N2, using tabulated values:

      DO 70 K=1,NBINS
        IF (ENER(K) .GE. EGR(1)) THEN
        DO 60 KG=1,12
          IF (ENER(K) .GE. EGR(KG) .AND. ENER(K) .LT. EGR(KG+1))
     >    TOTNEW(K)=TERPOO(ENER(K),EGR(KG),EGR(KG+1),SGR(KG),SGR(KG+1))
   60   CONTINUE
        ELSE
          TOTNEW(K)=TOTX(K)
        ENDIF
   70 CONTINUE

C Calculate ratio (=1 < 10 keV):

      DO 90 K = 1,NBINS
        IF (ENER(K) .GE. EGR(1)) THEN
          RATIO(K) = TOTX(K)/TOTNEW(K)
C         IF (RATIO(K) .GT. 1.) RATIO(K) = 1.
        ELSE
          RATIO(K) = 1.
        ENDIF
   90 CONTINUE

      RETURN

      END SUBROUTINE HEXC

!-----------------------------------------------------------------------

      REAL FUNCTION TERPOO(X,X1,X2,Y1,Y2)

      implicit none

      real,intent(in) :: x,x1,x2,y1,y2

      TERPOO = EXP ( ALOG(Y1) + ALOG(X/X1)*ALOG(Y2/Y1)/ALOG(X2/X1) )

      RETURN

      END FUNCTION TERPOO

max_stars_repo_name: scivision/NCAR-GLOW

max_issues_repo_name: space-physics/NCAR-GLOW

max_forks_repo_name: space-physics/NCAR-GLOW


program exercise_2_7
   implicit none

   character(*), parameter    :: input_file = "../data/input.txt", output_file = "output.txt"
   character(:), allocatable  :: fmt  !формат для ввода-вывода
   integer                    :: In = 0, Out = 0
   real(4)                      :: X , Y , Z

   open (file=input_file, newunit=In)
      read (In, *) X, Y
   close (In)

   open (file=output_file, newunit=Out)
      fmt = "(a, '=', F0.2, '  ', a, '=', F0.2 )"
      write (Out, fmt) "x", X , "y", Y
  ! close (Out)  !На самом деле можно оставить файл открытым и дальнейшие действия выполнять при открытом файле
   !Много раз открывать-закрывать файл ввода-вывода тоже не хорошо

   open (file=output_file, newunit=Out, position='append')
   if(Y>X) then
      z = (ABS(x)+ ABS(y))/2  !|11|+|-33| /2 =22
      write (Out,'(a)') "x<y"
   else if( X>Y.or.Y==X) then
      z = 1+ 2*ABS(x)    !1+2*|-11|=23
      write (Out,'(a)') "x>=y"
   end if
      write (Out,* ) "z=", Z

   close (Out)

end program exercise_2_7


text:
!>
!! @file
!! This file is part of SeisSol.
!!
!! @author Michael Dumbser (michael.dumbser AT unitn.it, https://www5.unitn.it/People/it/Web/Persona/PER0029602#INFO)
!! @author Sebastian Rettenberger (sebastian.rettenberger AT tum.de, http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger)
!!
!! @section LICENSE
!! Copyright (c) 2007-2016, SeisSol Group
!! All rights reserved.
!!
!! Redistribution and use in source and binary forms, with or without
!! modification, are permitted provided that the following conditions are met:
!!
!! 1. Redistributions of source code must retain the above copyright notice,
!!    this list of conditions and the following disclaimer.
!!
!! 2. Redistributions in binary form must reproduce the above copyright notice,
!!    this list of conditions and the following disclaimer in the documentation
!!    and/or other materials provided with the distribution.
!!
!! 3. Neither the name of the copyright holder nor the names of its
!!    contributors may be used to endorse or promote products derived from this
!!    software without specific prior written permission.
!!
!! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
!! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
!! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
!! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
!! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!! POSSIBILITY OF SUCH DAMAGE.

#include <Initializer/preProcessorMacros.fpp>

MODULE ini_SeisSol_mod
  !--------------------------------------------------------------------------
  USE TypesDef
  !--------------------------------------------------------------------------
  IMPLICIT NONE
  PRIVATE
  !----------------------------------------------------------------------------
  INTERFACE ini_SeisSol
     MODULE PROCEDURE ini_SeisSol
  END INTERFACE

  !----------------------------------------------------------------------------
  PUBLIC  :: ini_SeisSol
  !----------------------------------------------------------------------------

CONTAINS

  SUBROUTINE ini_SeisSol(time,timestep,pvar,cvar,EQN,IC,MESH,MPI,      &
       SOURCE,DISC,BND,OptionalFields,IO, &
       programTitle) !
    !--------------------------------------------------------------------------
    USE COMMON_readpar_mod
    USE ini_OptionalFields_mod
    USE dg_setup_mod
    USE ini_MODEL_mod
 !   USE DGSponge_mod
    USE calc_deltaT_mod

    use MeshReaderCBinding
#ifdef HDF
    USE HDF5
#endif
    use jacobiNormal_mod
    use iso_c_binding, only: c_loc, c_null_char
    use f_ftoc_bind_interoperability
    !--------------------------------------------------------------------------
    IMPLICIT NONE                                                              !
    !--------------------------------------------------------------------------
#ifdef PARALLEL
    INCLUDE 'mpif.h'
#endif
    REAL                           :: time,x,y,Variable(8),k1,k2               !
    INTEGER                        :: timestep                                 !
    REAL,POINTER                   :: pvar(:,:)                                !
    REAL,POINTER                   :: cvar(:,:)                                !
    TYPE (tEquations)              :: EQN                                      !
    TYPE (tInitialCondition)       :: IC                                       !
    TYPE (tUnstructMesh)           :: MESH                                     !
    TYPE (tMPI), OPTIONAL          :: MPI                                      !
    TYPE (tSource)                 :: SOURCE                                   !
    TYPE (tDiscretization)         :: DISC                                     !
    TYPE (tUnstructOptionalFields) :: OptionalFields                           !
    TYPE (tInputOutput)            :: IO                                       !
    TYPE (tBoundary)               :: BND                                      !
    CHARACTER(LEN=100)             :: programTitle                             !
    ! local variable declaration                                               !
    INTEGER                        :: iElem, iSide, iVtx, i                    !
    INTEGER                        :: counter
    INTEGER                        :: allocStat,iErr                           !
    INTEGER                        :: dummy(4)                                 !
    INTEGER                        :: nGraphVertex                             !
    INTEGER                        :: LocPoly                                  !
	INTEGER						   :: IntegrationMask(1:9)					   !
    INTEGER, POINTER               :: MetisWeight(:)                           !
    REAL                           :: AnelasticFactor(7,10)                    ! Factor of additional cost due to anelasticity
    REAL                           :: OrderFactor(7)                           ! Factor of additional cost due to order
    CHARACTER(LEN=256)             :: outfile
    CHARACTER(LEN=5)               :: cmyrank
    CHARACTER(LEN=600)             :: Filename, InFileName, OutFileName
    LOGICAL                        :: exists
    REAL, ALLOCATABLE              :: MaterialVals(:)
    REAL, ALLOCATABLE              :: WaveSpeeds(:)
#ifdef HDF
! HDF5 variables
     CHARACTER(LEN=100) :: file_name = "" ! File name
     CHARACTER(LEN=8) :: dsetname   ! Dataset name

     INTEGER(HID_T) :: file_id       ! File identifier
     INTEGER(HID_T) :: dset_id       ! Dataset identifier
     INTEGER(HID_T) :: filespace     ! Dataspace identifier in file
     INTEGER(HID_T) :: memspace      ! Dataspace identifier in memory
     INTEGER(HID_T) :: plist_id      ! Property list identifier

     INTEGER(HSIZE_T), DIMENSION(2) :: dimsf ! Dataset dimensions.
     INTEGER(HSIZE_T), DIMENSION(2) :: dimsfi
     INTEGER        :: info   ! Number of variables to output
     INTEGER(HSIZE_T), DIMENSION(2) :: count
     INTEGER(HSSIZE_T), DIMENSION(2) :: offset
     REAL (4), ALLOCATABLE :: data (:,:)  ! Data to write
     INTEGER :: rank = 2 ! Dataset rank
#endif
     INTEGER :: error, error_n  ! Error flags
     INTEGER  :: t1, t2, clock_max      ! needed for dwalltime
     REAL     :: clock_rate = 0.0d0         ! needed for dwalltime
! HDF5 variables end
     integer :: l_ruptureFace, l_side, l_elementId, l_localFaceId ! loop counter of the dynamic rupture faces
     real :: l_jInv
     real, dimension(9, 9) :: l_A, l_B
     real, dimension(3, 3) :: l_attenuation

    !--------------------------------------------------------------------------
    INTENT(IN)                     :: programTitle                             !
    INTENT(INOUT)                  :: EQN,DISC,IO                              ! Some values are set in the TypesDef
    INTENT(OUT)                    :: IC,MESH,SOURCE,BND, OptionalFields       !
    INTENT(OUT)                    :: time,timestep             !

    ! register ini_seissol function
    EPIK_FUNC_REG("ini_SeisSol")
    SCOREP_USER_FUNC_DEFINE()
    ! register read_mesh/compute mesh region
    EPIK_USER_REG(r_read_compute_mesh, "read_compute_mesh")
    SCOREP_USER_REGION_DEFINE( r_read_compute_mesh )
    !--------------------------------------------------------------------------

    ! start epik/scorep function ini_SeisSol
    EPIK_FUNC_START()
    SCOREP_USER_FUNC_BEGIN( "ini_SeisSol" )

    time                     = 0.0                                             ! Initialize
    timestep                 = 0                                               ! Global iteration number, incl. restarts
    DISC%iterationstep       = 0                                               ! local iteration number (0 after restart)
    DISC%LoopCPUTime         = 0.0                                             ! Initialize CPU-time-measurement
    EQN%pi                   = ACOS(-1.)                                       ! Initialize
    MESH%MaxElementsOnVertex = 20                                              ! Initialize
    MESH%nZones              = 1                                               ! Default values
    !                                                                          !
    IO%ContourFile = 'contour'                                                 ! File for plot body contour
    !
    !
    CALL readpar(                                         &                    ! read the parameter file
         EQN          = EQN                             , &                    !
         IC           = IC                              , &                    !
         usMESH       = MESH                            , &                    !
         DISC         = DISC                            , &                    !
         SOURCE       = SOURCE                          , &                    !
         BND          = BND                             , &                    !
         IO           = IO                              , &                    !
         programTitle = programTitle                    , &                    !
         MPI          = MPI                               )

    !
    IF(MPI%nCPU.GT.1) THEN
     DISC%IterationCriterion = 3 ! Associated to local time step
    ENDIF
    !

    ! Set ouput and checkpoint parameters
    ! This has to be done before the LTS setup!
    if( io%format .eq. 6 ) then
      call c_interoperability_enableWaveFieldOutput( i_waveFieldInterval = io%outInterval%timeInterval, &
                                                     i_waveFieldFilename = trim(io%OutputFile) // c_null_char )
    endif

    if( io%checkpoint%interval .gt. 0 ) then
        call c_interoperability_enableCheckPointing( i_checkPointInterval = io%checkpoint%interval, &
                                                     i_checkPointFilename = trim(io%checkpoint%filename) // c_null_char, &
                                                     i_checkPointBackend = trim(io%checkpoint%backend) // c_null_char )
    endif

#ifdef INTEGRATE_QUANTITIES
	do i = 1,9
		if ( io%IntegrationMask(i) ) then
			IntegrationMask(i) = 1
		else
			IntegrationMask(i) = 0
		end if
	end do

	call c_interoperability_getIntegrationMask( i_integrationMask = IntegrationMask(1:9) )
#endif

    ! Start mesh reading/computing section
    EPIK_USER_START(r_read_compute_mesh)
    SCOREP_USER_REGION_BEGIN( r_read_compute_mesh, "read_compute_mesh", SCOREP_USER_REGION_TYPE_COMMON )
    call read_mesh_fast(IO,EQN,DISC,MESH,BND,MPI)
    !                                                                          !
    ! output neighbour list
    !OPEN(UNIT=999,FILE='Neighbourhood.dat')
    !   DO iElem = 1,MESH%nElem
    !         WRITE(999,*) MESH%ELEM%SideNeighbor(:,iElem)
    !   ENDDO
    !CLOSE(999)

    logInfo(*) '<--------------------------------------------------------->' !
    logInfo(*) '<           Calling DG Initialization level 1             >' !
    logInfo(*) '<--------------------------------------------------------->' !

    CALL iniGalerkin3D_us_level1_new(                 &              !
        EQN    = EQN                                , &              !
        DISC   = DISC                               , &              !
        MESH   = MESH                               , &              !
        BND    = BND                                , &              !
        IC     = IC                                 , &              !
        SOURCE = SOURCE                             , &              !
        OptionalFields = OptionalFields             , &              !
        MPI    = MPI                                , &              !
        IO     = IO                                   )              !

    ! End mesh reading/computing section
    EPIK_USER_END(r_read_compute_mesh)
    SCOREP_USER_REGION_END( r_read_compute_mesh )

    logInfo(*) 'Generated Kernels: Checking boundary conditions'

    do iElem = 1, mesh%nElem
      do iSide = 1,4
        ! write an error message in the case of not supported boundary conditions
        if(       ( mesh%elem%reference( iSide, iElem ) .ne. 0 )& ! no boundary conditions: insides the computational domain
            .and. ( mesh%elem%reference( iSide, iElem ) .ne. 6 )& ! periodic boundary conditions
            .and. ( mesh%elem%reference( iSide, iElem ) .ne. 5 )& ! absorbing boundary conditions
            .and. ( mesh%elem%reference( iSide, iElem ) .ne. 1 )& ! free surface boundary conditions
            .and. ( mesh%elem%reference( iSide, iElem ) .ne. 3 )& ! dynamic rupture boundary conditions
          ) then
          logError(*) 'boundary conditions not supported'
          stop
        endif
      enddo
    enddo

    !                                                                               !
    DISC%NGP = 1                                                               !
    !                                                                          !
    !                                                                          !
    ALLOCATE(                                          &                       !
         pvar( MESH%nElem , EQN%nVar),                 &                       ! primitive variables
         cvar( MESH%nElem , EQN%nVar),                 &                       ! conservative var.
         STAT=allocStat                                )                       !
    !                                                                          !
    IF (allocStat .NE. 0) THEN                                                 ! Error Handling
       logError(*) 'could not allocate all variables!'      !
       STOP                                                                    ! STOP
    END IF                                                                     !
    !                                                                          !
    pvar(:,:)  = 0.                                                            ! Initialize
    cvar(:,:)  = 0.                                                            ! Initialize
    !                                                                          !
    !                                                                          ! Compute the number of background variables
    CALL ini_OptionalFields(                              &                    ! Allociert: OptionalFields%dt
         OptionalFields = OptionalFields                , &                    !            OptionalFields%RK_k
         SOURCE         = SOURCE                        , &                    !            OptionalFields%Residual
         EQN            = EQN                           , &                    !            OptionalFields%Backgroundvalue
         MESH           = MESH                          , &                    !
         DISC           = DISC                          , &                    !
         IO             = IO                              )                    !
    !                                                                          !
    !
    !T. Ulrich 08.2015 Read 2D basis Function for Vr output
    CALL Read2dGF(DISC,IO)
    !
    IF (EQN%linearized) THEN                                                   !
       CALL ini_MODEL(                                      &                  ! Initialize Local Linearized calculation
            MaterialVal    = OptionalFields%BackgroundValue,&                  ! Initialize Local Linearized calculation
            EQN            = EQN                          , &                  ! Initialize Local Linearized calculation
            MESH           = MESH                         , &                  ! Initialize Local Linearized calculation
            IO             = IO                           , &                  ! Initialize Local Linearized calculation
            DISC           = DISC                         , &                  ! Initialize Local Linearized calculation
            BND            = BND                            )                  ! Initialize Local Linearized calculation
    END IF                                                                     !
    logInfo(*) '<--------------------------------------------------------->' !
    logInfo(*) '<           Calling DG Initialization level 2             >' !
    logInfo(*) '<--------------------------------------------------------->' !
    !                                                                       !
    CALL iniGalerkin3D_us_level2_new(                    &              !
         EQN    = EQN                                , &              !
         DISC   = DISC                               , &              !
         MESH   = MESH                               , &              !
         BND    = BND                                , &              !
         IC     = IC                                 , &              !
         SOURCE = SOURCE                             , &              !
         OptionalFields = OptionalFields             , &              !
         MPI    = MPI                                , &              !
         IO     = IO                                   )              !
    !                                                                       !
    logInfo(*) 'Galerkin module initialized correctly.'    !
    !
    !aheineck: Metisweighs are not used, @TODO we should delete this
#if 0
    IF(IO%METISWeights) THEN ! not yet done for hybrids
        ! Compute Metis weights
        logInfo(*) 'Computing Metis weights ...'

        IF(BND%periodic.NE.0)THEN
           logError(*) 'MetisWeights can only be computed for non-periodic boundary conditions!'
           STOP
        ENDIF

        ALLOCATE(MetisWeight(MESH%nElem))

        !OrderFactor(:)       = (/4.34, 5.36, 8.39, 16.69, 37.08, 83.69, 180.70/)
        OrderFactor(:)       = (/1.00, 1.07, 1.29, 1.60, 2.62, 5.07, 10.29/)

        AnelasticFactor(1,:) = (/1.38, 1.48, 1.62, 1.68, 1.82, 1.88, 1.95, 2.08, 2.18, 2.27/)
        AnelasticFactor(2,:) = (/1.38, 1.48, 1.62, 1.68, 1.82, 1.88, 1.95, 2.08, 2.18, 2.27/)
        AnelasticFactor(3,:) = (/1.49, 1.70, 1.91, 2.07, 2.28, 2.49, 2.66, 2.82, 2.99, 3.20/)
        AnelasticFactor(4,:) = (/1.77, 2.07, 2.43, 2.70, 2.99, 3.34, 3.69, 4.03, 4.47, 4.67/)
        AnelasticFactor(5,:) = (/1.92, 2.30, 2.83, 3.27, 3.84, 4.38, 4.83, 5.40, 5.77, 6.32/)
        AnelasticFactor(6,:) = (/1.93, 2.44, 3.00, 3.55, 4.13, 4.65, 5.06, 5.82, 6.44, 6.87/)
        AnelasticFactor(7,:) = (/1.93, 2.44, 3.00, 3.55, 4.13, 4.65, 5.06, 5.82, 6.44, 6.87/)

        DO iElem = 1,MESH%nElem
           LocPoly = INT(DISC%Galerkin%LocPoly(iElem))
!           IF(EQN%LocAnelastic(iElem).EQ.1)THEN
!                  MetisWeight(iElem) = CEILING( (OrderFactor(LocPoly+1)                            &
!                                              * (AnelasticFactor(LocPoly+1,EQN%nMechanisms))) *100 )
!           ELSE
                  MetisWeight(iElem) = CEILING( OrderFactor(LocPoly+1) *100 )
!           ENDIF

        ENDDO

        !MetisWeight(:) = CEILING( REAL(MetisWeight(:) / minval(MetisWeight(:))) )
        logInfo(*) 'Minimum Metis weight:  ', minval(MetisWeight(:))
        logInfo(*) 'Maximum Metis weight:  ', maxval(MetisWeight(:))

        ! Read Metis graph file and output weighted Metis graph file
        InFileName  = TRIM(IO%MetisFile) // '.dgraph'
        OutFileName = TRIM(IO%MetisFile) // '.dgraph.weighted'

        OPEN(UNIT=998,FILE=TRIM(InFileName))
        OPEN(UNIT=999,FILE=TRIM(OutFileName))
        READ(998,*)  dummy(1), dummy(2)
        WRITE(999,*) dummy(1), dummy(2), 10, 1
        nGraphVertex = dummy(1)
        IF(nGraphVertex.EQ.MESH%nElem)THEN
           DO iElem = 1,MESH%nElem
              dummy(:) = 0
              counter  = 0
              DO iSide = 1,MESH%GlobalElemType
                 IF(MESH%ELEM%SideNeighbor(iSide,iElem).NE.MESH%nElem+1)THEN
                     counter = counter+1
                 ENDIF
              ENDDO
              READ(998,*)  dummy(1:counter)
              WRITE(999,*) MetisWeight(iElem), dummy(1:counter)
           ENDDO
        ELSE
           PRINT *, 'Number of Vertices in Graph File not equal to number of element of the mesh!'
        ENDIF

        CLOSE(998)
        CLOSE(999)

        logInfo(*) 'Metis weights computed successfully!'
        logInfo(*) 'weighted graph file written to ',TRIM(IO%MetisFile) // '.dgraph.weighted'
        logInfo(*) 'Terminating Programm !'
        STOP
    ENDIF
#endif

    ! TODO do we need to call this function when we read a checkpoint?? (SR)
    CALL icGalerkin3D_us_new(EQN, DISC, MESH, IC, SOURCE, IO)
    ! Pre-compute basis functions at Gauss points of non-conforming boundaries
 !   CALL NonConformingGPEvaluation3D(DISC, MESH, IO, BND)
 !   CALL DGSponge(0., EQN, DISC, MESH, BND, IC, SOURCE, IO) ! not yet done for hybrids/unified version
    !

    logInfo(*) 'Galerkin3D module initialized correctly.'     !
    logInfo(*) ' '                                             !
    !                                                                          !
    logInfo(*) '<--------------------------------------------------------->'  !
    logInfo(*) '<     ini_SeisSol finished successfully                   >'  !
    logInfo(*) '<--------------------------------------------------------->'  !
    !                                                                          !

    ! end epik function ini_SeisSol
    EPIK_FUNC_END()
    SCOREP_USER_FUNC_END()
  END SUBROUTINE ini_SeisSol                                                    !

END MODULE ini_SeisSol_mod

max_stars_repo_name: ivotron/SeisSol

max_issues_repo_name: ivotron/SeisSol

max_forks_repo_name: ivotron/SeisSol