chem/KPP/kpp/kpp-2.1/int/oldies/rodas3.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         CALL RODAS3(NVAR,TIN,TOUT,STEPMIN,STEPMAX,STEPMIN,
  17      +                   VAR,ATOL,RTOL,
  18      +                   Info,FUNC_CHEM,JAC_CHEM)
  19
  20       RETURN
  21       END
  22
  23
  24
  25       SUBROUTINE RODAS3(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     Stiffly accurate Rosenbrock 3(2), with
  33 C     stiffly accurate 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, March 1996
  63 C     The Center for Global and Regional Environmental Research
  64
  65       KPP_REAL   K1(NVAR), K2(NVAR), K3(NVAR), K4(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, H, Hold, Tplus
  71       KPP_REAL   ERR, factor, facmax
  72       KPP_REAL   c43, tau, x1, x2, ytol, elo
  73
  74       INTEGER    n,nfcn,njac,Naccept,Nreject,i,j,ier
  75       INTEGER    Info(5)
  76       LOGICAL    IsReject,Autonomous
  77       EXTERNAL   FUNC_CHEM, JAC_CHEM
  78
  79 c     Initialization of counters, etc.
  80       Autonomous = Info(1) .EQ. 1
  81       uround = 1.d-15
  82       c43 = - 8.d0/3.d0
  83       H = DMAX1(1.d-8, Hstart)
  84       Hmin = DMAX1(Hmin,uround*(Tnext-T))
  85       Hmax = DMIN1(Hmax,Tnext-T)
  86       Tplus = T
  87       IsReject = .false.
  88       Naccept  = 0
  89       Nreject  = 0
  90       Nfcn     = 0
  91       Njac     = 0
  92
  93
  94 C === Starting the time loop ===
  95  10    continue
  96        Tplus = T + H
  97        if ( Tplus .gt. Tnext ) then
  98           H = Tnext - T
  99           Tplus = Tnext
 100        end if
 101
 102        CALL JAC_CHEM(NVAR, T, y, JAC)
 103        Njac = Njac+1
 104        gHinv = -2.0d0/H
 105        do 20 j=1,NVAR
 106          JAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) + gHinv
 107  20    continue
 108        CALL KppDecomp (JAC, ier)
 109
 110        if (ier.ne.0) then
 111          if ( H.gt.Hmin) then
 112             H = 5.0d-1*H
 113             go to 10
 114          else
 115             print *,'IER <> 0, H=',H
 116             stop
 117          end if
 118        end if
 119
 120        CALL FUNC_CHEM(NVAR, T, y, F1)
 121
 122 C ====== NONAUTONOMOUS CASE ===============
 123        IF (.not. Autonomous) THEN
 124          tau = DSQRT( uround*DMAX1( 1.0d-5, DABS(T) ) )
 125          CALL FUNC_CHEM(NVAR, T+tau, y, K2)
 126          nfcn=nfcn+1
 127          do 30 j = 1,NVAR
 128            K3(j) = ( K2(j)-F1(j) )/tau
 129  30      continue
 130
 131 C ----- STAGE 1 (NONAUTONOMOUS) -----
 132          x1 = 0.5*H
 133          do 40 j = 1,NVAR
 134            K1(j) =  F1(j) + x1*K3(j)
 135  40      continue
 136          CALL KppSolve (JAC, K1)
 137
 138 C ----- STAGE 2 (NONAUTONOMOUS) -----
 139          x1 = 4.d0/H
 140          x2 = 1.5d0*H
 141          do 50 j = 1,NVAR
 142            K2(j) = F1(j) - x1*K1(j) + x2*K3(j)
 143  50      continue
 144          CALL KppSolve (JAC, K2)
 145
 146 C ====== AUTONOMOUS CASE ===============
 147        ELSE
 148 C ----- STAGE 1 (AUTONOMOUS) -----
 149          do 60 j = 1,NVAR
 150            K1(j) =  F1(j)
 151  60      continue
 152          CALL KppSolve (JAC, K1)
 153
 154 C ----- STAGE 2 (AUTONOMOUS) -----
 155          x1 = 4.d0/H
 156          do 70 j = 1,NVAR
 157            K2(j) = F1(j) - x1*K1(j)
 158  70      continue
 159          CALL KppSolve (JAC, K2)
 160        END IF
 161
 162 C ----- STAGE 3 -----
 163        do 80 j = 1,NVAR
 164          ynew(j) = y(j) - 2.0d0*K1(j)
 165  80    continue
 166        CALL FUNC_CHEM(NVAR, T+H, ynew, F1)
 167        nfcn=nfcn+1
 168        do 90 j = 1,NVAR
 169          K3(j) = F1(j) + ( -K1(j) + K2(j) )/H
 170  90    continue
 171        CALL KppSolve (JAC, K3)
 172
 173 C ----- STAGE 4 -----
 174        do 100 j = 1,NVAR
 175          ynew(j) = y(j) - 2.0d0*K1(j) - K3(j)
 176  100   continue
 177        CALL FUNC_CHEM(NVAR, T+H, ynew, F1)
 178        nfcn=nfcn+1
 179        do 110 j = 1,NVAR
 180          K4(j) = F1(j) + ( -K1(j) + K2(j) - C43*K3(j)  )/H
 181  110   continue
 182        CALL KppSolve (JAC, K4)
 183
 184 C ---- The Solution ---
 185
 186        do 120 j = 1,NVAR
 187          ynew(j) = y(j) - 2.0d0*K1(j) - K3(j) - K4(j)
 188  120   continue
 189
 190
 191 C ====== Error estimation ========
 192
 193         ERR=0.d0
 194         do 130 i=1,NVAR
 195            ytol = AbsTol(i) + RelTol(i)*DABS(ynew(i))
 196            ERR = ERR + ( K4(i)/ytol )**2
 197  130    continue
 198         ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )
 199
 200 C ======= Choose the stepsize ===============================
 201         elo    = 3.0D0 ! estimator local order
 202         factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
 203         Hnew   = DMIN1(Hmax,DMAX1(Hmin, H/factor))
 204
 205 C ======= Rejected/Accepted Step ============================
 206
 207         IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
 208           IsReject = .true.
 209           H = DMIN1(H/10,Hnew)
 210           Nreject  = Nreject+1
 211         ELSE
 212           DO 140 i=1,NVAR
 213              y(i)  = ynew(i)
 214  140      CONTINUE
 215           T = Tplus
 216           IF (.NOT.IsReject) THEN
 217               H = Hnew   ! Do not increase stepsize if previos step was rejected
 218           END IF
 219           IsReject = .false.
 220           Naccept = Naccept+1
 221         END IF
 222
 223
 224 C ======= End of the time loop ===============================
 225       if ( T .lt. Tnext ) go to 10
 226
 227
 228
 229 C ======= Output Information =================================
 230       Info(2) = Nfcn
 231       Info(3) = Njac
 232       Info(4) = Naccept
 233       Info(5) = Nreject
 234       Hstart  = H
 235
 236       RETURN
 237       END
 238
 239
 240
 241       SUBROUTINE FUNC_CHEM(N, T, Y, P)
 242       INCLUDE 'KPP_ROOT_params.h'
 243       INCLUDE 'KPP_ROOT_global.h'
 244       INTEGER N
 245       KPP_REAL   T, Told
 246       KPP_REAL   Y(NVAR), P(NVAR)
 247       Told = TIME
 248       TIME = T
 249       CALL Update_SUN()
 250       CALL Update_RCONST()
 251       CALL Fun( Y, FIX, RCONST, P )
 252       TIME = Told
 253       RETURN
 254       END
 255
 256
 257       SUBROUTINE JAC_CHEM(N, T, Y, J)
 258       INCLUDE 'KPP_ROOT_params.h'
 259       INCLUDE 'KPP_ROOT_global.h'
 260       INTEGER N
 261       KPP_REAL   Told, T
 262       KPP_REAL   Y(NVAR), J(LU_NONZERO)
 263       Told = TIME
 264       TIME = T
 265       CALL Update_SUN()
 266       CALL Update_RCONST()
 267       CALL Jac_SP( Y, FIX, RCONST, J )
 268       TIME = Told
 269       RETURN
 270       END
 271