Numérotation des lignes :

C FUNSHO    SOURCE    BECC      09/12/07    21:15:35     6579      subroutine funsho(nordpo, Tmaxcv, acv_r, Rgas_r,     &     P_r, T_r, u_r, logr,     &     T, P, rho, u, d)CC     INPUTCC     nordpo           = order of polynomial for cp and cv (see alsoC                        Tmaxcv)CC     Tmaxcv           = maximum temperature for cv polynomial expansionC                        cv(T) = cv(Tmaxcv) if T > TmaxcvCC     acv_r            = to compute cv and ether for r (and rs);C                        vector such thatC                        cv = \sum_{i=1,nordpo+1} acv(i) T^{i-1}C                        ether = \sum_{i=1,nordpo+1} acv(i) T^{i} / (i)CC     Rgas_r           = gas constant for r and (rs)CC     p_r, T_r, u_r    = pressure, temperature and velocities for rCC     logr             =  .true. right travelling shockC                      =  .false. left travelling shockCC     T                =  temperature for rsCC     OUTPUTCC     p, uC     rho, d           = pressure, velociy, density andC                        shock speed for rsC                        NB. for a non-er shock, d = u + cC                            with c = speed sound in rsCc      implicit none      integer nordpo      real*8  rindr      real*8 Tmaxcv, acv_r(1:(nordpo+1)), Rgas_r     &     , P_r, T_r, u_r     &     , T, P, rho, u, d      real*8 un_r, rho_r     &     , et, cv, et_r, cv_r, ht, ht_r, b, c, csi, b2      logical logrC      if (logr) then         rindr = 1.0D0      else         rindr = -1.0D0      endif      un_r = rindr * u_rC      rho_r = P_r / (Rgas_r * T_r)      call prith1(nordpo, acv_r, Tmaxcv, T_r, et_r, cv_r)      ht_r = et_r + (Rgas_r * T_r)      call prith1(nordpo, acv_r, Tmaxcv, T, et, cv)      ht = et + (Rgas_r * T)C      b = ((ht - ht_r) / Rgas_r) + (0.5d0 * (T_r - T))      c = T_r      b2 = b * b      if (((c * T) .lt. (1d-4 * b2)) .and. (b .lt. 0.0D0)) then         csi = (0.5d0 * c / b2) - (0.125d0 * (T * ((c / b2)**2))) +     &        (0.0625d0 *  (T*T * ((c / b2)**3)))         csi = csi * abs(b)      else         csi = sqrt(b2 + (c * T)) + b         csi = csi / T      endifC      rho = rho_R * csi      P = Rgas_r * rho * T      u = abs (2.0D0 * (ht - ht_r) *     &     (csi - 1.0D0) / (csi + 1.0D0))      u = un_r + (sign(1.0D0,(csi - 1.0D0)) * sqrt(u))C      if (T .lt. T_r) then         d = (cv + Rgas_r) / cv         d = d * Rgas_r * T         d = u + sqrt(d)      else         if (abs (T - T_r) .lt. (1d-5 * T_r)) then            d = 2.0D0 * (cv_r + Rgas_r) / cv_r            d = d * Rgas_r * T_r * csi * csi / (csi + 1.0D0)         else            d = 2.0D0 * (ht - ht_r) * csi * csi / ((csi * csi) - 1.0D0)         endif         d = abs(d)         d = un_r + sqrt(d)      endifC      u = rindr * u      d = rindr * dC      return      end

© Cast3M 2003 - Tous droits réservés.
Mentions légales