chem/KPP/kpp/kpp-2.1/int/oldies/qssa1.f

   1 C -- QSSA WITH STEADY STATE APPROXIMATION --
   2 C    For plain QSSA (to remove the steady state assumption)
   3 C    modify slow -> 0, fast -> 1e20
   4 C
   5
   6       SUBROUTINE INTEGRATE( TIN, TOUT  )
   7
   8       INCLUDE 'KPP_ROOT_params.h'
   9       INCLUDE 'KPP_ROOT_global.h'
  10
  11 C TIN - Start Time
  12       KPP_REAL TIN
  13 C TOUT - End Time
  14       KPP_REAL TOUT
  15
  16 C Local variables
  17       KPP_REAL P_VAR(NVAR), D_VAR(NVAR), V1(NVAR), V2(NVAR)
  18       LOGICAL IsReject
  19       KPP_REAL T, Tnext, STEP, STEPold, Told, SUP
  20       KPP_REAL ERR, ERRold, ratio, factor, facmax, tmp
  21       INTEGER I
  22       KPP_REAL slow, fast
  23
  24       T = TIN
  25       Tnext = TOUT
  26       STEP = DMAX1(STEPMIN,1.d-10)
  27       Told = T
  28       SUP  = 1e-14
  29       IsReject = .false.
  30       ERR = 1.d0
  31       ERRold = 1.d0
  32       slow = 0.01
  33       fast = 10.
  34
  35  10    continue
  36        Tplus = T + STEP
  37        if ( Tplus .gt. Tnext ) then
  38           STEP = Tnext - T
  39           Tplus = Tnext
  40        end if
  41
  42
  43         TITI = TIME
  44         TIME = T
  45         CALL Update_SUN()
  46         CALL Update_RCONST()
  47         TIME = TITI
  48         CALL FSPLIT_VAR ( VAR,  P_VAR, D_VAR )
  49
  50         do i=1,NVAR
  51           IF ( D_VAR(i)*step .lt. slow) THEN       ! SLOW SPECIES
  52              XXX = STEP * (P_VAR(I) - D_VAR(I)*VAR(I))
  53              V1(I) = VAR(I) + XXX
  54              V2(I) = VAR(I) + 0.5*XXX
  55           ELSE IF ( D_VAR(i)*step .GT. fast) THEN  ! FAST SPECIES
  56              V1(I) = P_VAR(I)/D_VAR(I)
  57              V2(I) = V1(I)
  58           ELSE                                    ! MEDIUM LIVED
  59             if ( abs(D_VAR(i)).gt.SUP ) then
  60               ratio = P_VAR(i)/D_VAR(i)
  61               tmp = exp(-D_VAR(i)*STEP*0.5)
  62               V1(i) = tmp * tmp * (VAR(i) - ratio) + ratio
  63               V2(i) = tmp * (VAR(i) - ratio) + ratio
  64             else
  65               tmp = D_VAR(i) * STEP * 0.5
  66               V1(i) = VAR(i) + P_VAR(i) * STEP * ( 1 - tmp *
  67      *             ( 1 - 2.0 / 3.0 * tmp ) )
  68               V2(i) = VAR(i) + P_VAR(i) * 0.5 * STEP*( 1-0.5*tmp*
  69      *             ( 1 - 1.0 / 3.0 * tmp ) )
  70               end if
  71           END IF
  72         end do
  73
  74         TITI = TIME
  75         TIME = T + 0.5*STEP
  76         CALL Update_SUN()
  77         CALL Update_RCONST()
  78         TIME = TITI
  79         CALL FSPLIT_VAR ( V2,  P_VAR, D_VAR )
  80
  81         do i=1,NVAR
  82           IF ( D_VAR(i)*step .lt. slow) THEN       ! SLOW SPECIES
  83              XXX = STEP * (P_VAR(I) - D_VAR(I)*VAR(I))
  84              V2(I) = V2(I) + 0.5*XXX
  85           ELSE IF ( D_VAR(i)*step .GT. fast) THEN  ! FAST SPECIES
  86              V2(I) = P_VAR(I)/D_VAR(I)
  87           ELSE                                    ! MEDIUM LIVED
  88             if ( abs(D_VAR(i)).gt.SUP ) then
  89               ratio = P_VAR(i)/D_VAR(i)
  90               tmp = exp(-D_VAR(i)*STEP*0.5)
  91               V2(i) = tmp * (V2(i) - ratio) + ratio
  92             else
  93               tmp = D_VAR(i) * STEP * 0.5
  94               V2(i) = V2(i) + P_VAR(i) * 0.5 * STEP * ( 1 - 0.5 * tmp *
  95      *             ( 1 - 1.0 / 3.0 * tmp ) )
  96             end if
  97           END IF
  98         end do
  99
 100 C ===== Extrapolation and error estimation ========
 101
 102         ERRold=ERR
 103         ERR=0.0D0
 104           do i=1,NVAR
 105              ERR = ERR + ((V2(i)-V1(i))/(ATOL(i) + RTOL(i)*V2(i)))**2
 106           end do
 107         ERR = DSQRT( ERR/NVAR )
 108         STEPold=STEP
 109
 110 C ===== choosing the stepsize =====================
 111
 112         factor = 0.9*ERR**(-0.35)*ERRold**0.2
 113         if (IsReject) then
 114             facmax=1.
 115         else
 116             facmax=8.
 117         end if
 118         factor = DMAX1( 1.25D-1, DMIN1(factor,facmax) )
 119         STEP = DMIN1( STEPMAX, DMAX1(STEPMIN,factor*STEP) )
 120
 121 C===================================================
 122
 123         if ( (ERR.gt.1).and.(STEPold.gt.STEPMIN) ) then
 124           IsReject = .true.
 125         else
 126           IsReject = .false.
 127           do 140 i=1,NVAR
 128              VAR(i)  = DMAX1(V2(i), 0.d0)
 129  140      continue
 130           T = Tplus
 131         end if
 132       if ( T .lt. Tnext ) go to 10
 133
 134       TIME = Tnext
 135       RETURN
 136       END