Original WRF subgrid support version from John Michalakes without fire
[wrffire.git] / wrfv2_fire / chem / KPP / kpp / kpp-2.1 / int / oldies / exqssa.f
blob925d4a7896d3933e2939321c04c0b4ba4f362d1a
1 C -- EXTRAPOLATED QSSA WITH STEADY STATE APPROXIMATION --
2 C For extrapolated plain QSSA (to remove the steady state assumption)
3 C modify slow -> 0, fast -> 1e20
6 SUBROUTINE INTEGRATE( TIN, TOUT )
8 INCLUDE 'KPP_ROOT_params.h'
9 INCLUDE 'KPP_ROOT_global.h'
11 C TIN - Start Time
12 KPP_REAL TIN
13 C TOUT - End Time
14 KPP_REAL TOUT
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
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.
35 10 continue
36 Tplus = T + STEP
37 if ( Tplus .gt. Tnext ) then
38 STEP = Tnext - T
39 Tplus = Tnext
40 end if
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 )
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
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 )
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
100 C ===== Extrapolation and error estimation ========
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
110 C ===== choosing the stepsize =====================
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) )
121 C===================================================
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(2*V2(i)-V1(i), 0.d0)
129 140 continue
130 T = Tplus
131 end if
132 if ( T .lt. Tnext ) go to 10
134 TIME = Tnext
135 RETURN
136 END