integ.f

Go to the documentation of this file.

!                   ****************
                    SUBROUTINE integ
!                   ****************
!
     &( a , b , iein , npoin)
!
!***********************************************************************
! SISYPHE   V6P1                                   21/07/2011
!***********************************************************************
!
!brief    COMPUTES THE EINSTEIN INTEGRAL
!+                FOR SUSPENDED TRANSPORT.
!code
!+                              A-1     /1       A
!+                             B       |  / 1-Y  \     33Y
!+           I = 0.216*1.83 * ------ * | | ----- | LN --- DY
!+                                 A   |  \  Y   /      B
!+                            (1-B)    /B
!+
!+                          WS                    KR
!+      WITH    A = -------------------  AND  B = --
!+                  KAPPA * BETA * U*CW           H
!
!note     VALUES - CALCULATED BY SIMPSON - ARE TABULATED
!+         AND LINEARLY INTERPOLATED.
!
!history
!+        2000
!+        V5P5
!+   MODIFICATIONS EBS/TB
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A              |-->| ROUSE NUMBER
!| B              |-->| RATIO BEDLOAD LAYER/WATER DEPTH
!| IEIN           |<--| INTEGRAL VALUE I
!| NPOIN          |-->| NUMBER OF POINTS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NPOIN
      DOUBLE PRECISION, INTENT(IN)    :: A(npoin)
      DOUBLE PRECISION, INTENT(INOUT) :: IEIN(npoin), B(npoin)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, PARAMETER :: LNA = 15, lnb = 42
      SAVE einj3,aein,bein
!
      DOUBLE PRECISION EINJ3(lnb,lna)
      DOUBLE PRECISION AEIN(lna),BEIN(lnb)
!
      DOUBLE PRECISION A1, A2, B1, B2, INT1, INT2
!
      INTEGER I, IT, JT, INF, SUP, MIL
!
! **********************************************************************
!        TABLES
! **********************************************************************
!
! RATIO QSS/QSC
! COMPUTED USING EINSTEIN INTEGRALS
! AND GENERATED BY MAPLE V
!
! DATA FOR A=  .0
!
      einj3(:,1) = (/
     & .0000d+00, .4806d+00, .9944d+00, .1535d+01, .2386d+01, .3889d+01,
     & .7128d+01, .9708d+01, .1422d+02, .1798d+02, .2055d+02, .2383d+02,
     & .2813d+02, .3399d+02, .4243d+02, .5549d+02, .7811d+02, .1257d+03,
     & .1628d+03, .2269d+03, .2797d+03, .3155d+03, .3609d+03, .4202d+03,
     & .5005d+03, .6152d+03, .7913d+03, .1093d+04, .1721d+04, .2206d+04,
     & .3037d+04, .3716d+04, .4176d+04, .4756d+04, .5511d+04, .6531d+04,
     & .7982d+04, .1020d+05, .1398d+05, .4627d+05, .9803d+05, .5538d+06/
     & )
!
! DATA FOR A=  .1
!
      einj3(:,2) = (/
     & .0000d+00, .4298d+00, .8772d+00, .1338d+01, .2048d+01, .3271d+01,
     & .5813d+01, .7778d+01, .1113d+02, .1386d+02, .1569d+02, .1801d+02,
     & .2101d+02, .2505d+02, .3077d+02, .3944d+02, .5407d+02, .8379d+02,
     & .1063d+03, .1442d+03, .1747d+03, .1951d+03, .2207d+03, .2536d+03,
     & .2977d+03, .3595d+03, .4526d+03, .6082d+03, .9202d+03, .1154d+04,
     & .1545d+04, .1857d+04, .2066d+04, .2326d+04, .2660d+04, .3105d+04,
     & .3727d+04, .4659d+04, .6207d+04, .1844d+05, .3649d+05, .1758d+06/
     & )
!
! DATA FOR A=  .2
!
      einj3(:,3) = (/
     & .0000d+00, .3881d+00, .7823d+00, .1180d+01, .1782d+01, .2792d+01,
     & .4821d+01, .6344d+01, .8874d+01, .1089d+02, .1223d+02, .1390d+02,
     & .1604d+02, .1888d+02, .2282d+02, .2869d+02, .3835d+02, .5729d+02,
     & .7124d+02, .9408d+02, .1121d+03, .1239d+03, .1386d+03, .1573d+03,
     & .1820d+03, .2160d+03, .2663d+03, .3481d+03, .5065d+03, .6219d+03,
     & .8095d+03, .9561d+03, .1052d+04, .1171d+04, .1322d+04, .1521d+04,
     & .1793d+04, .2193d+04, .2841d+04, .7579d+04, .1400d+05, .5756d+05/
     & )
!
! DATA FOR A=  .4
!
      einj3(:,4) = (/
     & .0000d+00, .3241d+00, .6390d+00, .9459d+00, .1393d+01, .2110d+01,
     & .3460d+01, .4418d+01, .5934d+01, .7092d+01, .7840d+01, .8753d+01,
     & .9895d+01, .1137d+02, .1335d+02, .1619d+02, .2063d+02, .2877d+02,
     & .3441d+02, .4320d+02, .4981d+02, .5405d+02, .5919d+02, .6557d+02,
     & .7375d+02, .8468d+02, .1001d+03, .1241d+03, .1674d+03, .1970d+03,
     & .2429d+03, .2771d+03, .2989d+03, .3253d+03, .3580d+03, .3997d+03,
     & .4551d+03, .5332d+03, .6534d+03, .1408d+04, .2273d+04, .6813d+04/
     & )
!
! DATA FOR A=  .6
!
      einj3(:,5) = (/
     & .0000d+00, .2773d+00, .5367d+00, .7818d+00, .1128d+01, .1659d+01,
     & .2603d+01, .3238d+01, .4199d+01, .4903d+01, .5346d+01, .5877d+01,
     & .6525d+01, .7340d+01, .8403d+01, .9866d+01, .1205d+02, .1579d+02,
     & .1823d+02, .2186d+02, .2446d+02, .2609d+02, .2802d+02, .3036d+02,
     & .3327d+02, .3703d+02, .4215d+02, .4969d+02, .6239d+02, .7059d+02,
     & .8262d+02, .9120d+02, .9653d+02, .1028d+03, .1104d+03, .1199d+03,
     & .1320d+03, .1484d+03, .1724d+03, .3023d+03, .4277d+03, .9411d+03/
     & )
!
! DATA FOR A=  .8
!
      einj3(:,6) = (/
     & .0000d+00, .2419d+00, .4607d+00, .6618d+00, .9376d+00, .1346d+01,
     & .2033d+01, .2474d+01, .3114d+01, .3564d+01, .3841d+01, .4166d+01,
     & .4555d+01, .5031d+01, .5635d+01, .6435d+01, .7576d+01, .9407d+01,
     & .1054d+02, .1213d+02, .1323d+02, .1390d+02, .1467d+02, .1559d+02,
     & .1670d+02, .1808d+02, .1990d+02, .2244d+02, .2644d+02, .2887d+02,
     & .3226d+02, .3457d+02, .3597d+02, .3758d+02, .3948d+02, .4177d+02,
     & .4462d+02, .4832d+02, .5346d+02, .7757d+02, .9715d+02, .1601d+03/
     & )
!
! DATA FOR A= 1.0
!
      einj3(:,7) = (/
     & .0000d+00, .2142d+00, .4023d+00, .5711d+00, .7968d+00, .1120d+01,
     & .1638d+01, .1956d+01, .2400d+01, .2702d+01, .2884d+01, .3093d+01,
     & .3338d+01, .3631d+01, .3992d+01, .4455d+01, .5086d+01, .6039d+01,
     & .6595d+01, .7344d+01, .7837d+01, .8129d+01, .8461d+01, .8844d+01,
     & .9296d+01, .9844d+01, .1053d+02, .1145d+02, .1281d+02, .1358d+02,
     & .1461d+02, .1528d+02, .1567d+02, .1611d+02, .1662d+02, .1722d+02,
     & .1794d+02, .1884d+02, .2003d+02, .2487d+02, .2817d+02, .3657d+02/
     & )
!
! DATA FOR A= 1.2
!
      einj3(:,8) = (/
     & .0000d+00, .1919d+00, .3563d+00, .5006d+00, .6892d+00, .9512d+00,
     & .1354d+01, .1592d+01, .1912d+01, .2122d+01, .2246d+01, .2387d+01,
     & .2549d+01, .2738d+01, .2965d+01, .3247d+01, .3615d+01, .4141d+01,
     & .4432d+01, .4807d+01, .5044d+01, .5180d+01, .5332d+01, .5504d+01,
     & .5701d+01, .5933d+01, .6215d+01, .6572d+01, .7066d+01, .7332d+01,
     & .7668d+01, .7877d+01, .7997d+01, .8129d+01, .8277d+01, .8447d+01,
     & .8644d+01, .8880d+01, .9178d+01, .1024d+02, .1084d+02, .1207d+02/
     & )
!
! DATA FOR A= 1.5
!
      einj3(:,9) = (/
     & .0000d+00, .1658d+00, .3032d+00, .4205d+00, .5694d+00, .7685d+00,
     & .1058d+01, .1220d+01, .1428d+01, .1560d+01, .1635d+01, .1718d+01,
     & .1812d+01, .1918d+01, .2041d+01, .2188d+01, .2369d+01, .2609d+01,
     & .2732d+01, .2881d+01, .2970d+01, .3020d+01, .3074d+01, .3133d+01,
     & .3198d+01, .3271d+01, .3356d+01, .3456d+01, .3583d+01, .3645d+01,
     & .3719d+01, .3763d+01, .3786d+01, .3812d+01, .3839d+01, .3869d+01,
     & .3903d+01, .3941d+01, .3986d+01, .4117d+01, .4174d+01, .4259d+01/
     & )
!
! DATA FOR A= 2.0
!
      einj3(:,10) = (/
     & .0000d+00, .1349d+00, .2418d+00, .3298d+00, .4372d+00, .5737d+00,
     & .7582d+00, .8542d+00, .9708d+00, .1040d+01, .1078d+01, .1120d+01,
     & .1164d+01, .1213d+01, .1268d+01, .1329d+01, .1399d+01, .1483d+01,
     & .1522d+01, .1566d+01, .1590d+01, .1603d+01, .1616d+01, .1630d+01,
     & .1645d+01, .1661d+01, .1678d+01, .1697d+01, .1718d+01, .1727d+01,
     & .1737d+01, .1742d+01, .1745d+01, .1748d+01, .1751d+01, .1754d+01,
     & .1757d+01, .1760d+01, .1764d+01, .1772d+01, .1774d+01, .1777d+01/
     & )
!
! DATA FOR A= 2.5
!
      einj3(:,11) = (/
     & .0000d+00, .1135d+00, .2004d+00, .2699d+00, .3524d+00, .4531d+00,
     & .5821d+00, .6459d+00, .7199d+00, .7621d+00, .7849d+00, .8089d+00,
     & .8344d+00, .8615d+00, .8907d+00, .9222d+00, .9568d+00, .9955d+00,
     & .1012d+01, .1031d+01, .1040d+01, .1045d+01, .1050d+01, .1055d+01,
     & .1060d+01, .1066d+01, .1072d+01, .1077d+01, .1084d+01, .1086d+01,
     & .1089d+01, .1090d+01, .1091d+01, .1091d+01, .1092d+01, .1093d+01,
     & .1093d+01, .1094d+01, .1095d+01, .1096d+01, .1097d+01, .1097d+01/
     & )
!
! DATA FOR A= 3.0
!
      einj3(:,12) = (/
     & .0000d+00, .9783d-01, .1708d+00, .2278d+00, .2939d+00, .3724d+00,
     & .4689d+00, .5148d+00, .5664d+00, .5950d+00, .6101d+00, .6259d+00,
     & .6424d+00, .6596d+00, .6777d+00, .6969d+00, .7173d+00, .7392d+00,
     & .7484d+00, .7580d+00, .7630d+00, .7655d+00, .7681d+00, .7707d+00,
     & .7733d+00, .7759d+00, .7786d+00, .7813d+00, .7841d+00, .7853d+00,
     & .7864d+00, .7870d+00, .7872d+00, .7875d+00, .7878d+00, .7881d+00,
     & .7884d+00, .7887d+00, .7890d+00, .7896d+00, .7897d+00, .7898d+00/
     & )
!
! DATA FOR A= 4.0
!
      einj3(:,13) = (/
     & .0000d+00, .7657d-01, .1314d+00, .1730d+00, .2196d+00, .2727d+00,
     & .3345d+00, .3624d+00, .3927d+00, .4088d+00, .4172d+00, .4258d+00,
     & .4346d+00, .4437d+00, .4530d+00, .4626d+00, .4725d+00, .4828d+00,
     & .4870d+00, .4913d+00, .4935d+00, .4946d+00, .4957d+00, .4968d+00,
     & .4979d+00, .4990d+00, .5001d+00, .5012d+00, .5023d+00, .5028d+00,
     & .5033d+00, .5035d+00, .5036d+00, .5037d+00, .5038d+00, .5039d+00,
     & .5041d+00, .5042d+00, .5043d+00, .5045d+00, .5046d+00, .5046d+00/
     & )
!
! DATA FOR A= 7.0
!
      einj3(:,14) = (/
     & .0000d+00, .4618d-01, .7721d-01, .9958d-01, .1235d+00, .1493d+00,
     & .1771d+00, .1890d+00, .2013d+00, .2076d+00, .2108d+00, .2141d+00,
     & .2174d+00, .2207d+00, .2240d+00, .2274d+00, .2308d+00, .2343d+00,
     & .2357d+00, .2371d+00, .2378d+00, .2381d+00, .2385d+00, .2388d+00,
     & .2392d+00, .2396d+00, .2399d+00, .2403d+00, .2406d+00, .2408d+00,
     & .2409d+00, .2410d+00, .2410d+00, .2410d+00, .2411d+00, .2411d+00,
     & .2412d+00, .2412d+00, .2412d+00, .2413d+00, .2413d+00, .2413d+00/
     & )
!
! DATA FOR A=15.0
!
      einj3(:,15) = (/
     & .0000d+00, .2236d-01, .3656d-01, .4639d-01, .5654d-01, .6702d-01,
     & .7786d-01, .8230d-01, .8681d-01, .8909d-01, .9023d-01, .9138d-01,
     & .9254d-01, .9369d-01, .9486d-01, .9602d-01, .9720d-01, .9837d-01,
     & .9884d-01, .9932d-01, .9955d-01, .9967d-01, .9979d-01, .9991d-01,
     & .1000d+00, .1001d+00, .1003d+00, .1004d+00, .1005d+00, .1005d+00,
     & .1006d+00, .1006d+00, .1006d+00, .1006d+00, .1007d+00, .1007d+00,
     & .1007d+00, .1007d+00, .1007d+00, .1007d+00, .1007d+00, .1007d+00/
     & )
!
! VALUES OF TABLE BEIN
!
      bein = (/
     &      1.d0  ,.75d0  ,.6d0   ,.5d0   ,.4d0    ,.3d0    ,
     &     .2d0   ,.16d0  ,.12d0  ,.1d0   ,.09d0   ,.08d0   ,
     &     .07d0  ,.06d0  ,.05d0  ,.04d0  ,.03d0   ,.02d0   ,
     &     .016d0 ,.012d0 ,.01d0  ,.009d0 ,.008d0  ,.007d0  ,
     &     .006d0 ,.005d0 ,.004d0 ,.003d0 ,.002d0  ,.0016d0 ,
     &     .0012d0,.001d0 ,.0009d0,.0008d0,.0007d0 ,.0006d0 ,
     &     .0005d0,.0004d0,.0003d0,.0001d0,.00005d0,.00001d0 /
     & )
!
! VALUES OF TABLE AEIN
!
      aein =   (/ 0.d0,.1d0,.2d0,.4d0,.6d0,.8d0,1.d0,1.2d0,1.5d0,
     &            2.d0 ,2.5d0,3.d0,4.d0,7.d0,15.d0 /)
!
      DO i=1,npoin
!
!
        IF (b(i).LT.bein(lnb)) b(i)=bein(lnb)
! TREATS B AND A
! =======================
!
! CASE WHERE SUSPENSION IS NEGLIGIBLE
!
!
        IF (a(i).GT.aein(lna).OR.(b(i).GE.1.d0)) THEN
          iein(i) = 0.d0
!
! VALUES OUTSIDE OF THE TABLE
!
        ELSEIF ((b(i).GT.1.d0).OR.(b(i).LT.bein(lnb)).OR.
     &          (a(i).LT.aein(1))) THEN
!
 100      FORMAT(/,1x,'TBIJFR : B EST SUPERIEUR A 1',/,
     &             3x,'B     = ',g10.4,/,
     &             3x,'NOEUD = ',i5,/,
     &             3x,'VERIFIER LA PERTINENCE DE LA FORMULE DE BIJKER'
     &               ,' DANS VOTRE CAS',/)
 101      FORMAT(/,1x,'TBIJFR : B EST INFERIEUR A ',g10.4,
     &                ' LIMITE DE LA ZONE TABULEE',/,
     &             3x,'B     = ',g10.4,/,
     &             3x,'NOEUD = ',i5,/,
     &             3x,'VERIFIER LA PERTINENCE DE LA FORMULE DE BIJKER'
     &               ,' DANS VOTRE CAS',/)
 102      FORMAT(/,1x,'TBIJFR : A EST INFERIEUR A ',g10.4,
     &                ' LIMITE DE LA ZONE TABULEE',/,
     &             3x,'A     = ',g10.4,/,
     &             3x,'NOEUD = ',i5,/,
     &             3x,'VERIFIER LA PERTINENCE DE LA FORMULE DE BIJKER'
     &               ,' DANS VOTRE CAS',/)
!
          IF (b(i).GT.1.d0) THEN
            WRITE(lu,100) b(i),i
          ELSEIF (b(i).LT.bein(lnb)) THEN
            WRITE(lu,101) bein(lnb),b(i),i
          ELSEIF (a(i).LT.aein(1)) THEN
            WRITE(lu,102) aein(1),a(i),i
          ENDIF
          CALL plante(0)
          stop
!
! VALUES IN THE TABLE
!
        ELSE
!
! LOOKS FOR THE BOUNDS FOR B
! ===============================
!
          inf = 1
          sup = lnb
!----START OF THE DO WHILE LOOP----
 51       mil = (inf + sup ) / 2
            IF (bein(mil).GT.b(i)) THEN
              inf = mil
            ELSE
              sup = mil
            ENDIF
          IF ((sup - inf).NE.1)  GOTO 51
!----END OF THE DO WHILE LOOP----
          jt=sup
          b1 = bein(inf)
          b2 = bein(sup)
!
! LOOKS FOR THE BOUNDS FOR A
! ===============================
!
          inf = 1
          sup = lna
!----START OF THE DO WHILE LOOP----
          DO
          mil = (inf + sup ) / 2
            IF (aein(mil).LT.a(i)) THEN
              inf = mil
            ELSE
              sup = mil
            ENDIF
          IF ((sup - inf).EQ.1) EXIT
          ENDDO
!----END OF THE DO WHILE LOOP----
          it=sup
          a1 = aein(inf)
          a2 = aein(sup)
!
! COMPUTES THE IEIN INTEGRAL
! ==========================
!
! INTERPOLATION WITH CONSTANT A1
          int1 = (einj3(jt-1,it-1)-einj3(jt,it-1))/(b1-b2)*
     &           (b(i)-b1) + einj3(jt-1,it-1)
! INTERPOLATION WITH CONSTANT A2
          int2 = (einj3(jt-1,it)  -einj3(jt,it)  )/(b1-b2)*
     &           (b(i)-b1) + einj3(jt-1,it)
! RE-INTERPOLATION BETWEEN THE TWO PRECEDING VALUES
          iein(i) = (int1-int2)/(a1-a2)*(a(i)-a1) + int1
!
        ENDIF
!
      ENDDO !I
!
      RETURN
      END SUBROUTINE integ