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

   1       SUBROUTINE INTEGRATE( TIN, TOUT )
   2
   3       INCLUDE 'KPP_ROOT_params.h'
   4       INCLUDE 'KPP_ROOT_global.h'
   5
   6 C TIN - Start Time
   7       KPP_REAL TIN
   8 C TOUT - End Time
   9       KPP_REAL TOUT
  10
  11       INTEGER    INFO(5)
  12
  13       EXTERNAL   FUNC_CHEM, JAC_CHEM
  14
  15       INFO(1) = Autonomous
  16
  17         CALL ROS3(NVAR,TIN,TOUT,STEPMIN,STEPMAX,
  18      +                   STEPMIN,VAR,ATOL,RTOL,
  19      +                   Info,FUNC_CHEM,JAC_CHEM)
  20
  21       RETURN
  22       END
  23
  24
  25       SUBROUTINE ROS3(N,T,Tnext,Hmin,Hmax,Hstart,
  26      +                   y,AbsTol,RelTol,
  27      +                   Info,FUNC_CHEM,JAC_CHEM)
  28       IMPLICIT NONE
  29       INCLUDE 'KPP_ROOT_params.h'
  30       INCLUDE 'KPP_ROOT_sparse.h'
  31
  32 C       L-stable Rosenbrock 3(2), with
  33 C     strongly A-stable embedded formula for error control.
  34 C
  35 C     All the arguments aggree with the KPP syntax.
  36 C
  37 C  INPUT ARGUMENTS:
  38 C     y = Vector of (NVAR) concentrations, contains the
  39 C         initial values on input
  40 C     [T, Tnext] = the integration interval
  41 C     Hmin, Hmax = lower and upper bounds for the selected step-size.
  42 C          Note that for Step = Hmin the current computed
  43 C          solution is unconditionally accepted by the error
  44 C          control mechanism.
  45 C     AbsTol, RelTol = (NVAR) dimensional vectors of
  46 C          componentwise absolute and relative tolerances.
  47 C     FUNC_CHEM = name of routine of derivatives. KPP syntax.
  48 C          See the header below.
  49 C     JAC_CHEM = name of routine that computes the Jacobian, in
  50 C          sparse format. KPP syntax. See the header below.
  51 C     Info(1) = 1  for  autonomous   system
  52 C             = 0  for nonautonomous system
  53 C
  54 C  OUTPUT ARGUMENTS:
  55 C     y = the values of concentrations at Tend.
  56 C     T = equals Tend on output.
  57 C     Info(2) = # of FUNC_CHEM calls.
  58 C     Info(3) = # of JAC_CHEM calls.
  59 C     Info(4) = # of accepted steps.
  60 C     Info(5) = # of rejected steps.
  61 C
  62 C     Adrian Sandu, April 1996
  63 C     The Center for Global and Regional Environmental Research
  64
  65       KPP_REAL K1(NVAR), K2(NVAR), K3(NVAR)
  66       KPP_REAL F1(NVAR), JAC(LU_NONZERO)
  67       KPP_REAL Hmin,Hmax,Hnew,Hstart,ghinv,uround
  68       KPP_REAL y(NVAR), ynew(NVAR)
  69       KPP_REAL AbsTol(NVAR), RelTol(NVAR)
  70       KPP_REAL T, Tnext, Tplus, H, elo
  71       KPP_REAL ERR, factor, facmax
  72       KPP_REAL gam, c21, c31, c32, b1, b2, b3
  73       KPP_REAL d1, d2, d3, a21, a31, a32
  74       KPP_REAL alpha2, alpha3, g1, g2, g3
  75       KPP_REAL tau, x1, x2, x3, dround, ytol
  76       INTEGER    n,nfcn,njac,Naccept,Nreject,i,j,ier
  77       INTEGER    Info(5)
  78       LOGICAL    IsReject,Autonomous
  79       EXTERNAL   FUNC_CHEM, JAC_CHEM
  80
  81       gam=   .43586652150845899941601945119356d+00
  82       c21=  -.10156171083877702091975600115545d+01
  83       c31=   .40759956452537699824805835358067d+01
  84       c32=   .92076794298330791242156818474003d+01
  85        b1=   .10000000000000000000000000000000d+01
  86        b2=   .61697947043828245592553615689730d+01
  87        b3=  -.42772256543218573326238373806514d+00
  88        d1=   .50000000000000000000000000000000d+00
  89        d2=  -.29079558716805469821718236208017d+01
  90        d3=   .22354069897811569627360909276199d+00
  91        a21 = 1.d0
  92        a31 = 1.d0
  93        a32 = 0.d0
  94        alpha2 = gam
  95        alpha3 = gam
  96        g1=   .43586652150845899941601945119356d+00
  97        g2=   .24291996454816804366592249683314d+00
  98        g3=   .21851380027664058511513169485832d+01
  99
 100
 101 c     Initialization of counters, etc.
 102       Autonomous = Info(1) .EQ. 1
 103       uround = 1.d-15
 104       dround = DSQRT(uround)
 105       IF (Hmax.le.0.D0) THEN
 106           Hmax = DABS(Tnext-T)
 107       END IF
 108       H = DMAX1(1.d-8, Hstart)
 109       Tplus = T
 110       IsReject = .false.
 111       Naccept  = 0
 112       Nreject  = 0
 113       Nfcn     = 0
 114       Njac     = 0
 115
 116 C === Starting the time loop ===
 117  10    continue
 118
 119        Tplus = T + H
 120        if ( Tplus .gt. Tnext ) then
 121           H = Tnext - T
 122           Tplus = Tnext
 123        end if
 124
 125        CALL JAC_CHEM(NVAR, T, y, JAC)
 126        Njac = Njac+1
 127        gHinv = -1.0d0/(gam*H)
 128        do 15 j=1,LU_NONZERO
 129          JAC(j) = -JAC(j)
 130  15    continue
 131        do 20 j=1,NVAR
 132          JAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) - gHinv
 133  20    continue
 134        CALL KppDecomp (JAC, ier)
 135
 136        if (ier.ne.0) then
 137          if ( H.gt.Hmin) then
 138             H = 5.0d-1*H
 139             go to 10
 140          else
 141             print *,'IER <> 0, H=',H
 142             stop
 143          end if
 144        end if
 145
 146        CALL FUNC_CHEM(NVAR, T, y, F1)
 147
 148 C ====== NONAUTONOMOUS CASE ===============
 149        IF (.not. Autonomous) THEN
 150          tau = DSIGN(dround*DMAX1( 1.0d-6, DABS(T) ), T)
 151          CALL FUNC_CHEM(NVAR, T+tau, y, K2)
 152          nfcn=nfcn+1
 153          do 30 j = 1,NVAR
 154            K3(j) = ( K2(j)-F1(j) )/tau
 155  30      continue
 156
 157 C ----- STAGE 1 (NONAUTONOMOUS) -----
 158          x1 = g1*H
 159          do 35 j = 1,NVAR
 160            K1(j) =  F1(j) + x1*K3(j)
 161  35      continue
 162          CALL KppSolve (JAC, K1)
 163
 164 C ----- STAGE 2 (NONAUTONOMOUS) -----
 165        do 40 j = 1,NVAR
 166          ynew(j) = y(j) + K1(j)
 167  40    continue
 168        CALL FUNC_CHEM(NVAR, T+gam*H, ynew, F1)
 169        nfcn=nfcn+1
 170        x1 = c21/H
 171        x2 = g2*H
 172        do 45 j = 1,NVAR
 173          K2(j) = F1(j) + x1*K1(j) + x2*K3(j)
 174  45    continue
 175        CALL KppSolve (JAC, K2)
 176
 177 C ----- STAGE 3  (NONAUTONOMOUS) -----
 178        x1 = c31/H
 179        x2 = c32/H
 180        x3 = g3*H
 181        do 50 j = 1,NVAR
 182          K3(j) = F1(j) + x1*K1(j) + x2*K2(j) + x3*K3(j)
 183  50    continue
 184        CALL KppSolve (JAC, K3)
 185
 186
 187 C ====== AUTONOMOUS CASE ===============
 188        ELSE
 189
 190 C ----- STAGE 1 (AUTONOMOUS) -----
 191          do 60 j = 1,NVAR
 192            K1(j) =  F1(j)
 193  60      continue
 194          CALL KppSolve (JAC, K1)
 195
 196 C ----- STAGE 2 (AUTONOMOUS) -----
 197        do 65 j = 1,NVAR
 198          ynew(j) = y(j) + K1(j)
 199  65    continue
 200        CALL FUNC_CHEM(NVAR, T + gam*H, ynew, F1)
 201        nfcn=nfcn+1
 202          x1 = c21/H
 203          do 70 j = 1,NVAR
 204            K2(j) = F1(j) + x1*K1(j)
 205  70      continue
 206          CALL KppSolve (JAC, K2)
 207
 208 C ----- STAGE 3  (AUTONOMOUS) -----
 209        x1 = c31/H
 210        x2 = c32/H
 211        do 90 j = 1,NVAR
 212          K3(j) = F1(j) + x1*K1(j) + x2*K2(j)
 213  90    continue
 214        CALL KppSolve (JAC, K3)
 215
 216        END  IF ! Autonomousous
 217
 218 C ---- The Solution ---
 219
 220        do 120 j = 1,NVAR
 221          ynew(j) = y(j) + b1*K1(j) + b2*K2(j) + b3*K3(j)
 222  120   continue
 223
 224
 225 C ====== Error estimation ========
 226
 227         ERR=0.d0
 228         do 130 i=1,NVAR
 229            ytol = AbsTol(i) + RelTol(i)*DMAX1(DABS(y(i)),DABS(ynew(i)))
 230            ERR=ERR+((d1*K1(i)+d2*K2(i)+d3*K3(i))/ytol)**2
 231  130    continue
 232         ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )
 233
 234 C ======= Choose the stepsize ===============================
 235
 236         elo    = 3.0D0 ! estimator local order
 237         factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
 238         Hnew   = DMIN1(Hmax,DMAX1(Hmin, H/factor))
 239
 240 C ======= Rejected/Accepted Step ============================
 241
 242         IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
 243           IsReject = .true.
 244           H = DMIN1(H/10,Hnew)
 245           Nreject  = Nreject+1
 246         ELSE
 247           DO 140 i=1,NVAR
 248              y(i)  = ynew(i)
 249  140      CONTINUE
 250           T = Tplus
 251           IF (.NOT.IsReject) THEN
 252               H = Hnew   ! Do not increase stepsize if previos step was rejected
 253           END IF
 254           IsReject = .false.
 255           Naccept = Naccept+1
 256         END IF
 257
 258
 259 C ======= End of the time loop ===============================
 260       if ( T .lt. Tnext ) go to 10
 261
 262
 263 C ======= Output Information =================================
 264       Info(2) = Nfcn
 265       Info(3) = Njac
 266       Info(4) = Naccept
 267       Info(5) = Nreject
 268       Hstart  = H
 269
 270       RETURN
 271       END
 272
 273
 274
 275       SUBROUTINE FUNC_CHEM(N, T, Y, P)
 276       INCLUDE 'KPP_ROOT_params.h'
 277       INCLUDE 'KPP_ROOT_global.h'
 278       INTEGER N
 279       KPP_REAL   T, Told
 280       KPP_REAL   Y(NVAR), P(NVAR)
 281       Told = TIME
 282       TIME = T
 283       CALL Update_SUN()
 284       CALL Update_RCONST()
 285       CALL Fun( Y,  FIX, RCONST, P )
 286       TIME = Told
 287       RETURN
 288       END
 289
 290       SUBROUTINE JAC_CHEM(N, T, Y, J)
 291       INCLUDE 'KPP_ROOT_params.h'
 292       INCLUDE 'KPP_ROOT_global.h'
 293       INTEGER N
 294       KPP_REAL   Told, T
 295       KPP_REAL   Y(NVAR), J(LU_NONZERO)
 296       Told = TIME
 297       TIME = T
 298       CALL Update_SUN()
 299       CALL Update_RCONST()
 300       CALL Jac_SP( Y,  FIX, RCONST, J )
 301       TIME = Told
 302       RETURN
 303       END
 304