source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/oldies/qssa1.f @ 4653

Last change on this file since 4653 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

File size: 3.9 KB
Line 
1C -- QSSA WITH STEADY STATE APPROXIMATION --
2C    For plain QSSA (to remove the steady state assumption)
3C    modify slow -> 0, fast -> 1e20
4C
5
6      SUBROUTINE INTEGRATE( TIN, TOUT  ) 
7
8      INCLUDE 'KPP_ROOT_params.h'
9      INCLUDE 'KPP_ROOT_global.h'
10
11C TIN - Start Time
12      KPP_REAL TIN
13C TOUT - End Time
14      KPP_REAL TOUT
15
16C 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
100C ===== 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
110C ===== 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
121C===================================================
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       
Note: See TracBrowser for help on using the repository browser.