chem/KPP/kpp/kpp-2.1/int/oldies/ros2.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 C--------------------------------
  16       INTEGER N_stepss, N_accepteds, N_rejecteds, N_jacs, ITOL, IERR
  17       SAVE N_stepss, N_accepteds, N_rejecteds, N_jacs
  18 C--------------------------------
  19
  20       INFO(1) = 1 ! Autonomous
  21       INFO(2) = 0
  22
  23       CALL ROS2(NVAR,TIN,TOUT,STEPMIN,STEPMAX,
  24      +                   STEPMIN,VAR,ATOL,RTOL,
  25      +                   Info,FUNC_CHEM,JAC_CHEM)
  26
  27 C--------------------------------
  28       N_stepss=N_stepss+Info(4)+Info(5)
  29       N_accepteds=N_accepteds+Info(4)
  30       N_rejecteds=N_rejecteds+Info(5)
  31       N_jacs=N_jacs+Info(3)
  32       PRINT*,'Step=',N_stepss,' Acc=',N_accepteds,' Rej=',N_rejecteds,
  33      &          ' Jac=',N_jacs
  34 C--------------------------------
  35
  36       RETURN
  37       END
  38
  39
  40
  41
  42       SUBROUTINE ROS2(N,T,Tnext,Hmin,Hmax,Hstart,
  43      +                   Y,AbsTol,RelTol,
  44      +                   Info,FUNC_CHEM,JAC_CHEM)
  45       IMPLICIT NONE
  46       INCLUDE 'KPP_ROOT_params.h'
  47       INCLUDE 'KPP_ROOT_sparse.h'
  48
  49 C  INPUT ARGUMENTS:
  50 C     Y = Vector of (NVAR) concentrations, contains the
  51 C         initial values on input
  52 C     [T, Tnext] = the integration interval
  53 C     Hmin, Hmax = lower and upper bounds for the selected step-size.
  54 C          Note that for Step = Hmin the current computed
  55 C          solution is unconditionally accepted by the error
  56 C          control mechanism.
  57 C     AbsTol, RelTol = (NVAR) dimensional vectors of
  58 C          componentwise absolute and relative tolerances.
  59 C     FUNC_CHEM = name of routine of derivatives. KPP syntax.
  60 C          See the header below.
  61 C     JAC_CHEM = name of routine that computes the Jacobian, in
  62 C          sparse format. KPP syntax. See the header below.
  63 C     Info(1) = 1  for  autonomous   system
  64 C             = 0  for nonautonomous system
  65 C     Info(2) = 1  for third order embedded formula
  66 C             = 0  for first order embedded formula
  67 C
  68 C   Note: Stage 3 used to build strongly A-stable order 3 formula for error control
  69 C            Embed3 = (Info(2).EQ.1)
  70 C         IF Embed3 = .true. THEN the third order embedded formula is used
  71 C                     .false. THEN a first order embedded  formula is used
  72 C
  73 C
  74 C  OUTPUT ARGUMENTS:
  75 C     Y = the values of concentrations at TEND.
  76 C     T = equals TEND on output.
  77 C     Info(2) = # of FUNC_CHEM CALLs.
  78 C     Info(3) = # of JAC_CHEM CALLs.
  79 C     Info(4) = # of accepted steps.
  80 C     Info(5) = # of rejected steps.
  81
  82       KPP_REAL  K1(NVAR), K2(NVAR), K3(NVAR)
  83       KPP_REAL  F1(NVAR), JAC(LU_NONZERO)
  84       KPP_REAL  DFDT(NVAR)
  85       KPP_REAL  Hmin,Hmax,Hnew,Hstart,ghinv,uround
  86       KPP_REAL  Y(NVAR), Ynew(NVAR)
  87       KPP_REAL  AbsTol(NVAR), RelTol(NVAR)
  88       KPP_REAL  T, Tnext, H, Hold, Tplus
  89       KPP_REAL  ERR, factor, facmax
  90       KPP_REAL  tau, beta, elo, dround, a21, c21, c31, c32
  91       KPP_REAL  gamma3, d1, d2, d3, gam
  92       INTEGER    n,nfcn,njac,Naccept,Nreject,i,j,ier
  93       INTEGER    Info(5)
  94       LOGICAL    IsReject, Autonomous, Embed3
  95       EXTERNAL    FUNC_CHEM, JAC_CHEM
  96
  97       KPP_REAL  gamma, m1, m2, alpha, beta1, beta2, delta, w, e
  98
  99 c     Initialization of counters, etc.
 100       Autonomous = Info(1) .EQ. 1
 101       Embed3  = Info(2) .EQ. 1
 102       uround  = 1.d-15
 103       dround  = dsqrt(uround)
 104       H = DMAX1(1.d-8, Hmin)
 105       Tplus = T
 106       IsReject = .false.
 107       Naccept  = 0
 108       Nreject  = 0
 109       Nfcn     = 0
 110       Njac     = 0
 111
 112 C     Method Parameters
 113       gamma = 1.d0 + 1.d0/sqrt(2.d0)
 114       a21   = - 1.d0/gamma
 115       m1 = -3.d0/(2.d0*gamma)
 116       m2 = -1.d0/(2.d0*gamma)
 117       c21 = -2.d0/gamma
 118       c31 = -1.0D0/gamma**2*(1.0D0-7.0D0*gamma+9.0D0*gamma**2)
 119      &         /(-1.0D0+2.0D0*gamma)
 120       c32 = -1.0D0/gamma**2*(1.0D0-6.0D0*gamma+6.0D0*gamma**2)
 121      &          /(-1.0D0+2.0D0*gamma)/2
 122       gamma3 = 0.5D0 - 2*gamma
 123       d1 = ((-9.0D0*gamma+8.0D0*gamma**2+2.0D0)/gamma**2/
 124      &          (-1.0D0+2*gamma))/6.0D0
 125       d2 = ((-1.0D0+3.0D0*gamma)/gamma**2/
 126      &          (-1.0D0+2.0D0*gamma))/6.0D0
 127       d3 = -1.0D0/(3.0D0*gamma)
 128
 129 C === Start the time loop ===
 130        DO WHILE (T .LT. Tnext)
 131
 132 10     CONTINUE
 133
 134        Tplus = T + H
 135        IF ( Tplus .gt. Tnext ) THEN
 136           H = Tnext - T
 137           Tplus = Tnext
 138        END IF
 139
 140        CALL JAC_CHEM( T, Y, JAC )
 141
 142        Njac = Njac+1
 143        ghinv = -1.0d0/(gamma*H)
 144        DO 20 j=1,NVAR
 145          JAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) + ghinv
 146  20    CONTINUE
 147        CALL KppDecomp (JAC, ier)
 148
 149        IF (ier.ne.0) THEN
 150          IF ( H.gt.Hmin) THEN
 151             H = 5.0d-1*H
 152             GO TO 10
 153          ELSE
 154             PRINT *,'IER <> 0, H=',H
 155             STOP
 156          END IF
 157        END IF
 158
 159        CALL FUNC_CHEM( T, Y, F1 )
 160
 161 C ====== NONAUTONOMOUS CASE ===============
 162        IF (.NOT. Autonomous) THEN
 163          tau = DSIGN(DROUND*DMAX1( 1.0d-6, DABS(T) ), T)
 164          CALL FUNC_CHEM( T+tau, Y, K2)
 165          nfcn=nfcn+1
 166          DO 30 j = 1,NVAR
 167            DFDT(j) = ( K2(j)-F1(j) )/tau
 168  30      CONTINUE
 169        END IF ! .NOT.Autonomous
 170
 171 C ----- STAGE 1 -----
 172        delta = gamma*H
 173        DO 40 j = 1,NVAR
 174           K1(j) =  F1(j)
 175  40    CONTINUE
 176        IF (.NOT.Autonomous) THEN
 177          DO 45 j = 1,NVAR
 178            K1(j) = K1(j) + delta*DFDT(j)
 179  45      CONTINUE
 180        END IF ! .NOT.Autonomous
 181        CALL KppSolve (JAC, K1)
 182
 183 C ----- STAGE 2  -----
 184        DO 50 j = 1,NVAR
 185          Ynew(j) = Y(j) + a21*K1(j)
 186  50    CONTINUE
 187        CALL FUNC_CHEM( T+H, Ynew, F1)
 188        nfcn=nfcn+1
 189        beta = -c21/H
 190        DO 55 j = 1,NVAR
 191          K2(j) = F1(j) + beta*K1(j)
 192  55    CONTINUE
 193        IF (.NOT.Autonomous) THEN
 194          delta = -gamma*H
 195          DO 56 j = 1,NVAR
 196            K2(j) = K2(j) + delta*DFDT(j)
 197  56      CONTINUE
 198        END IF ! .NOT.Autonomous
 199        CALL KppSolve (JAC, K2)
 200
 201 C ----- STAGE 3  -----
 202        IF (Embed3) THEN
 203        beta1 = -c31/H
 204        beta2 = -c32/H
 205        delta = gamma3*H
 206        DO 57 j = 1,NVAR
 207          K3(j) = F1(j) + beta1*K1(j) + beta2*K2(j)
 208  57    CONTINUE
 209        IF (.NOT.Autonomous) THEN
 210          DO 58 j = 1,NVAR
 211            K3(j) = K3(j) + delta*DFDT(j)
 212  58      CONTINUE
 213        END IF ! .NOT.Autonomous
 214        CALL KppSolve (JAC, K3)
 215        END IF ! Embed3
 216
 217
 218 C ---- The Solution ---
 219        DO 120 j = 1,NVAR
 220          Ynew(j) = Y(j) + m1*K1(j) + m2*K2(j)
 221  120   CONTINUE
 222
 223
 224 C ====== Error estimation ========
 225
 226         ERR=0.d0
 227         DO 130 i=1,NVAR
 228            w = AbsTol(i) + RelTol(i)*DMAX1(DABS(Y(i)),DABS(Ynew(i)))
 229            IF ( Embed3 ) THEN
 230              e = d1*K1(i) + d2*K2(i) + d3*K3(i)
 231            ELSE
 232              e = 1.d0/(2.d0*gamma)*(K1(i)+K2(i))
 233            END IF  ! Embed3
 234            ERR = ERR + ( e/w )**2
 235  130    CONTINUE
 236         ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )
 237
 238 C ======= Choose the stepsize ===============================
 239
 240         IF ( Embed3 ) THEN
 241             elo = 3.0D0 ! estimator local order
 242         ELSE
 243             elo = 2.0D0
 244         END IF
 245         factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
 246         Hnew   = DMIN1(Hmax,DMAX1(Hmin, H/factor))
 247
 248 C ======= Rejected/Accepted Step ============================
 249
 250         IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
 251           IsReject = .true.
 252           H = DMIN1(H/10,Hnew)
 253           Nreject  = Nreject+1
 254         ELSE
 255           DO 140 i=1,NVAR
 256              Y(i)  = Ynew(i)
 257  140      CONTINUE
 258           T = Tplus
 259           IF (.NOT.IsReject) THEN
 260               H = Hnew   ! Do not increase stepsize IF previous step was rejected
 261           END IF
 262           IsReject = .false.
 263           Naccept = Naccept+1
 264         END IF
 265
 266 C ======= END of the time loop ===============================
 267       END DO
 268
 269
 270 C ======= Output Information =================================
 271       Info(2) = Nfcn
 272       Info(3) = Njac
 273       Info(4) = Naccept
 274       Info(5) = Nreject
 275
 276       RETURN
 277       END
 278
 279
 280
 281       SUBROUTINE FUNC_CHEM( T, Y, P )
 282       INCLUDE 'KPP_ROOT_params.h'
 283       INCLUDE 'KPP_ROOT_global.h'
 284       KPP_REAL  T, Told
 285       KPP_REAL  Y(NVAR), P(NVAR)
 286       Told = TIME
 287       TIME = T
 288       CALL Update_SUN()
 289       CALL Update_RCONST()
 290       CALL Fun( Y,  FIX, RCONST, P )
 291       TIME = Told
 292       RETURN
 293       END
 294
 295
 296       SUBROUTINE JAC_CHEM( T, Y, J )
 297       INCLUDE 'KPP_ROOT_params.h'
 298       INCLUDE 'KPP_ROOT_global.h'
 299       KPP_REAL  Told, T
 300       KPP_REAL  Y(NVAR), J(LU_NONZERO)
 301       Told = TIME
 302       TIME = T
 303       CALL Update_SUN()
 304       CALL Update_RCONST()
 305       CALL Jac_SP( Y,  FIX, RCONST, J )
 306       TIME = Told
 307       RETURN
 308       END
 309
 310
 311
 312
 313
 314