LCOV - code coverage report
Current view: top level - math - outint.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 88 88 100.0 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             :       MODULE m_outint
       2             :       CONTAINS
       3  1103541822 :       SUBROUTINE outint(
       4      225766 :      >                  msh,e,fkap,cs,cis,s,vr,z,rn,rnot,h,d,
       5      225766 :      <                  a0,b0,a,b,
       6             :      <                  ki,nodes)
       7             : c
       8             : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
       9             : c-----x  perform the outward integration.   this routine is a     x----
      10             : c-----x  derivative of start1/diff1.   much unnecessary indexing  x----
      11             : c-----x  has been removed among other things.    dale koelling    x----
      12             : 
      13             : c     with adams' predictor and corrector. more stable than milne's
      14             : c       for GGA.
      15             : cc      Handbook of math.func. p.896. T.Asada Feb.20,1998.  
      16             : c-----x----x----x----x----x----x----x----x----x----x----x----x----x----
      17             : 
      18             :       IMPLICIT NONE
      19             : 
      20             : c     .. Arguments ..
      21             :       INTEGER,INTENT(IN)  :: msh
      22             :       INTEGER,INTENT(OUT) :: ki,nodes
      23             :       REAL   ,INTENT(IN)  :: e,fkap,cs,cis,s,z,rn,rnot,h,d
      24             :       REAL   ,INTENT(IN)  :: vr(msh)
      25             :       REAL   ,INTENT(OUT) :: a0(5),b0(5),a(msh),b(msh)
      26             : c     ..
      27             : c     .. local scalars ..
      28             :       REAL astr,atk,aw,bstr,btk,bw,da1,da2,da3,da4,db1,db2,db3,db4,g,gm,
      29             :      +     gp,h3,q11,q22,r,rp12,rp21,test,vt
      30             :       INTEGER i,k,kim,n
      31             :       REAL :: da(5),db(5)
      32             : c     ..
      33             : c     .. intrinsic functions ..
      34             :       INTRINSIC abs,sign,sqrt
      35             : c     ..
      36             : c.....------------------------------------------------------------------
      37             : c     .. equivalences ..
      38             :       EQUIVALENCE (da1,da(1)),(da2,da(2)),(da3,da(3)),(da4,da(4)),
      39             :      &  (da5,da(5))
      40             :       EQUIVALENCE (db1,db(1)),(db2,db(2)),(db3,db(3)),(db4,db(4)),
      41             :      &  (db5,db(5))
      42             : 
      43             :       REAL t18,t58,t98,t198,t378,t558,t598, da5,db5
      44             : c.....------------------------------------------------------------------
      45             : 
      46      225766 :       h3 = h/3.0e0
      47      225766 :       n = msh
      48             : cta+
      49      225766 :       t558=55.e0/8.e0
      50      225766 :       t598=59.e0/8.e0
      51      225766 :       t378=37.e0/8.e0
      52      225766 :       t98=9.e0/8.e0
      53             : 
      54      225766 :       t198=19.e0/8.e0
      55      225766 :       t58=5.e0/8.e0
      56      225766 :       t18=1.e0/8.e0
      57             : cta-
      58      225766 :       g = sqrt(fkap**2- (cis*z)**2)
      59      225766 :       gp = g + s*fkap
      60      225766 :       gm = g - s*fkap
      61      225766 :       q11 = -gm
      62      225766 :       q22 = -gp
      63             : c     25 april,1978
      64             : c     simple self-consistent starting procedure:
      65      225766 :       astr = sqrt(abs(gp))
      66      225766 :       bstr = sqrt(abs(gm))
      67      225766 :       r = rnot
      68     1354596 :       do 10 i = 1,5
      69     1128830 :          a(i) = astr
      70     1128830 :          b(i) = bstr
      71     1128830 :          rp21 = cis* (e*r-vr(i))
      72     1128830 :          rp12 = -rp21 - 2*cs*r
      73     1128830 :          a0(i) = h* (q11*astr+rp12*bstr)
      74     1128830 :          b0(i) = h* (q22*bstr+rp21*astr)
      75     1128830 :          r = r*d
      76      225766 :    10 continue
      77     2660716 :       test = 1.e-6
      78     2886482 :    20 vt = 0
      79     2886482 :       r = rnot
      80     2886482 :       r = rnot*d
      81    14432410 :       do 30 i = 2,5
      82    11545928 :          aw = a(i-1) + .5e0* (a0(i-1)+a0(i))
      83    11545928 :          bw = b(i-1) + .5e0* (b0(i-1)+b0(i))
      84    11545928 :          vt = vt + abs(aw-a(i))/astr + abs(bw-b(i))/bstr
      85    11545928 :          a(i) = 0.5e0* (a(i)+aw)
      86    11545928 :          b(i) = 0.5e0* (b(i)+bw)
      87    11545928 :          rp21 = cis* (e*r-vr(i))
      88    11545928 :          rp12 = -rp21 - 2*cs*r
      89    11545928 :          a0(i) = h* (q11*a(i)+rp12*b(i))
      90    11545928 :          b0(i) = h* (q22*b(i)+rp21*a(i))
      91    11545928 :          r = r*d
      92     2886482 :    30 continue
      93     2886482 :       if(vt.gt.test) go to 20
      94             :       r = rnot
      95     2483426 :       do 40 i = 1,5
      96     1128830 :          rp21 = cis* (e*r-vr(i))
      97     1128830 :          rp12 = -rp21 - 2*cs*r
      98     1128830 :          da(i) = h3* (q11*a(i)+rp12*b(i))
      99     1128830 :          db(i) = h3* (q22*b(i)+rp21*a(i))
     100     1128830 :          r = r*d
     101      225766 :    40 continue
     102             : c     end of starting procedure
     103             : 
     104      225766 :       nodes = 0
     105      225766 :       r = rn
     106             : 
     107             : c..... find the classical turning point........
     108             : 
     109    37491011 :       ki = n
     110    37491011 :    50 ki = ki - 1
     111    37491011 :       if(ki.le.10) go to 60
     112    37491011 :       r = r/d
     113    37491011 :       if(e*r.lt.vr(ki)) go to 50
     114      225766 :    60 ki = ki + 1
     115      225766 :       if(ki.ge.n) ki = ki - 1
     116      225766 :       kim = ki + 1
     117      225766 :       if(kim.gt.n) kim = n
     118      225766 :       r = rnot* (d**3)
     119      225766 :       k = 4
     120   137914507 :    70 k = k + 1
     121   137914507 :       r = r*d
     122   137914507 :       rp21 = cis* (e*r-vr(k))
     123   137914507 :       rp12 = -rp21 - 2.0e0*cs*r
     124             : 
     125             : c     adams' extrapolation formula for predictor.
     126   689572535 :       atk=a(k-1)+ t558*da4-t598*da3+t378*da2-t98*da1
     127   689572535 :       btk=b(k-1)+ t558*db4-t598*db3+t378*db2-t98*db1
     128             : 
     129   137914507 :       da1 = da2
     130   137914507 :       da2 = da3
     131   137914507 :       da3 = da4
     132             : 
     133   137914507 :       db1 = db2
     134   137914507 :       db2 = db3
     135   137914507 :       db3 = db4
     136             : 
     137   137914507 :       da4 = h3* (q11*atk+rp12*btk)
     138   137914507 :       db4 = h3* (q22*btk+rp21*atk)
     139             : 
     140             : c     adams interpolation formula for corrector.
     141   137914507 :       a(k)=a(k-1)+ t98*da4+t198*da3-t58*da2+t18*da1
     142   137914507 :       b(k)=b(k-1)+ t98*db4+t198*db3-t58*db2+t18*db1
     143             : 
     144   137914507 :       da4=h3*(q11*a(k)+rp12*b(k))
     145   137914507 :       db4=h3*(q22*b(k)+rp21*a(k))
     146             : 
     147   137914507 :       if(k.ge.kim) go to 80
     148             : 
     149   137688741 :       nodes = nodes + 0.501e0*abs(sign(1.0e0,a(k))-sign(1.0e0,a(k-1)))
     150   137688741 :       go to 70
     151             : 
     152             :    80 continue
     153             : 
     154      225766 :       RETURN
     155             :       END SUBROUTINE
     156             :       END

Generated by: LCOV version 1.13