chem/KPP/kpp/kpp-2.1/int/oldies/ros4.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 ROS4(NVAR,TIN,TOUT,STEPMIN,STEPMAX,
  18      +                   STEPMIN,VAR,ATOL,RTOL,
  19      +                   Info,FUNC_CHEM,JAC_CHEM)
  20
  21       RETURN
  22       END
  23
  24
  25
  26
  27       SUBROUTINE ROS4(N,T,Tnext,Hmin,Hmax,Hstart,
  28      +                   y,AbsTol,RelTol,
  29      +                   Info,FUNC_CHEM,JAC_CHEM)
  30       IMPLICIT NONE
  31       INCLUDE 'KPP_ROOT_params.h'
  32       INCLUDE 'KPP_ROOT_sparse.h'
  33 C
  34 C  Four Stages, Fourth Order L-stable Rosenbrock Method,
  35 C     with embedded L-stable, third order method for error control
  36 C     Simplified version of E. Hairer's atmros4; the coefficients are slightly
  37 C     different
  38 C
  39 C
  40 C  INPUT ARGUMENTS:
  41 C     y = Vector of (NVAR) concentrations, contains the
  42 C         initial values on input
  43 C     [T, Tnext] = the integration interval
  44 C     Hmin, Hmax = lower and upper bounds for the selected step-size.
  45 C          Note that for Step = Hmin the current computed
  46 C          solution is unconditionally accepted by the error
  47 C          control mechanism.
  48 C     AbsTol, RelTol = (NVAR) dimensional vectors of
  49 C          componentwise absolute and relative tolerances.
  50 C     FUNC_CHEM = name of routine of derivatives. KPP syntax.
  51 C          See the header below.
  52 C     JAC_CHEM = name of routine that computes the Jacobian, in
  53 C          sparse format. KPP syntax. See the header below.
  54 C     Info(1) = 1  for  Autonomous   system
  55 C             = 0  for nonAutonomous system
  56 C
  57 C  OUTPUT ARGUMENTS:
  58 C     y = the values of concentrations at Tend.
  59 C     T = equals TENDon output.
  60 C     Info(2) = # of FUNC_CHEM CALLs.
  61 C     Info(3) = # of JAC_CHEM CALLs.
  62 C     Info(4) = # of accepted steps.
  63 C     Info(5) = # of rejected steps.
  64 C     Hstart = The last accepted stepsize
  65 C
  66 C    Adrian Sandu, December 2001
  67 C
  68       KPP_REAL K1(NVAR), K2(NVAR), K3(NVAR), K4(NVAR)
  69       KPP_REAL F1(NVAR)
  70       KPP_REAL DFDT(NVAR)
  71       KPP_REAL JAC(LU_NONZERO)
  72       KPP_REAL Hmin,Hmax,Hstart
  73       KPP_REAL y(NVAR), ynew(NVAR)
  74       KPP_REAL AbsTol(NVAR), RelTol(NVAR)
  75       KPP_REAL T, Tnext, H, Hnew, Tplus
  76       KPP_REAL elo,ghinv,uround
  77       KPP_REAL ERR, factor, facmax
  78       KPP_REAL w, e, dround, tau
  79       KPP_REAL hgam1, hgam2, hgam3, hgam4
  80       KPP_REAL hc21, hc31, hc32, hc41, hc42, hc43
  81
  82       INTEGER    n,nfcn,njac,Naccept,Nreject,i,j,ier
  83       INTEGER    Info(5)
  84       LOGICAL    IsReject, Autonomous
  85       EXTERNAL   FUNC_CHEM, JAC_CHEM
  86
  87
  88 C   The method coefficients
  89       DOUBLE PRECISION gamma, gamma2, gamma3, gamma4
  90       PARAMETER ( gamma  =  0.5728200000000000D+00 )
  91       PARAMETER ( gamma2 = -0.1769193891319233D+01 )
  92       PARAMETER ( gamma3 =  0.7592633437920482D+00 )
  93       PARAMETER ( gamma4 = -0.1049021087100450D+00 )
  94       DOUBLE PRECISION a21, a31, a32, a41, a42, a43
  95       PARAMETER (  a21 = 0.2000000000000000D+01  )
  96       PARAMETER (  a31 = 0.1867943637803922D+01  )
  97       PARAMETER (  a32 = 0.2344449711399156D+00  )
  98       DOUBLE PRECISION alpha2, alpha3
  99       PARAMETER (  alpha2  = 0.1145640000000000D+01  )
 100       PARAMETER (  alpha3  = 0.6552168638155900D+00  )
 101       DOUBLE PRECISION c21, c31, c32, c41, c42, c43
 102       PARAMETER (  c21  = -0.7137615036412310D+01  )
 103       PARAMETER (  c31  =  0.2580708087951457D+01  )
 104       PARAMETER (  c32  =  0.6515950076447975D+00  )
 105       PARAMETER (  c41  = -0.2137148994382534D+01  )
 106       PARAMETER (  c42  = -0.3214669691237626D+00  )
 107       PARAMETER (  c43  = -0.6949742501781779D+00  )
 108       DOUBLE PRECISION b1, b2, b3, b4
 109       PARAMETER (  b1 = 0.2255570073418735D+01  )
 110       PARAMETER (  b2 = 0.2870493262186792D+00  )
 111       PARAMETER (  b3 = 0.4353179431840180D+00  )
 112       PARAMETER (  b4 = 0.1093502252409163D+01 )
 113       DOUBLE PRECISION d1, d2, d3, d4
 114       PARAMETER (  d1 = -0.2815431932141155D+00  )
 115       PARAMETER (  d2 = -0.7276199124938920D-01  )
 116       PARAMETER (  d3 = -0.1082196201495311D+00  )
 117       PARAMETER (  d4 = -0.1093502252409163D+01  )
 118
 119 c     Initialization of counters, etc.
 120       Autonomous = Info(1) .EQ. 1
 121       uround = 1.d-15
 122       dround = DSQRT(uround)
 123       IF (Hmax.le.0.D0) THEN
 124           Hmax = DABS(Tnext-T)
 125       END IF
 126       H = DMAX1(1.d-8, Hstart)
 127       Tplus = T
 128       IsReject = .false.
 129       Naccept  = 0
 130       Nreject  = 0
 131       Nfcn     = 0
 132       Njac     = 0
 133
 134 C === Starting the time loop ===
 135  10    CONTINUE
 136
 137        Tplus = T + H
 138        IF ( Tplus .gt. Tnext ) THEN
 139           H = Tnext - T
 140           Tplus = Tnext
 141        END IF
 142
 143 C              Initial Function and Jacobian values
 144        CALL FUNC_CHEM( T, y, F1 )
 145        CALL JAC_CHEM( T, y, JAC )
 146
 147 C              The time derivative for non-Autonomous case
 148        IF (.not. Autonomous) THEN
 149          tau = DSIGN(dround*DMAX1( 1.0d-6, DABS(T) ), T)
 150          CALL FUNC_CHEM( T+tau, y, K2 )
 151          nfcn=nfcn+1
 152          DO 20 j = 1,NVAR
 153            DFDT(j) = ( K2(j)-F1(j) )/tau
 154  20      CONTINUE
 155        END IF
 156
 157 C              Form the Prediction matrix and compute its LU factorization
 158        Njac = Njac+1
 159        ghinv = 1.0d0/(gamma*H)
 160        DO 30 j=1,LU_NONZERO
 161          JAC(j) = -JAC(j)
 162  30    CONTINUE
 163        DO 40 j=1,NVAR
 164          JAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) + ghinv
 165  40    CONTINUE
 166        CALL KppDecomp (JAC, ier)
 167 C
 168        IF (ier.ne.0) THEN
 169          IF ( H.gt.Hmin) THEN
 170             H = 5.0d-1*H
 171             GO TO 10
 172          ELSE
 173             PRINT *,'ROS4: Singular factorization at T=',T,'; H=',H
 174             STOP
 175          END IF
 176        END IF
 177
 178
 179 C ------------ STAGE 1-------------------------
 180        DO 50 j = 1,NVAR
 181          K1(j) =  F1(j)
 182  50    CONTINUE
 183        IF (.NOT. Autonomous) THEN
 184           hgam1 = H*gamma
 185           DO 60 j=1,NVAR
 186             K1(j) = K1(j) + hgam1*DFDT(j)
 187  60       CONTINUE
 188        END IF
 189        CALL KppSolve (JAC, K1)
 190
 191 C ----------- STAGE 2 -------------------------
 192        DO 70 j = 1,NVAR
 193          ynew(j) = y(j) + a21*K1(j)
 194  70    CONTINUE
 195        CALL FUNC_CHEM( T+alpha2*H, ynew, F1)
 196        nfcn=nfcn+1
 197        hc21 = c21/H
 198        DO 80 j = 1,NVAR
 199          K2(j) = F1(j) + hc21*K1(j)
 200  80    CONTINUE
 201        IF (.NOT. Autonomous) THEN
 202          hgam2 = H*gamma2
 203          DO 90 j=1,NVAR
 204             K2(j) = K2(j) + hgam2*DFDT(j)
 205  90      CONTINUE
 206        END IF
 207        CALL KppSolve (JAC, K2)
 208
 209
 210 C ------------ STAGE 3 -------------------------
 211        DO 100 j = 1,NVAR
 212          ynew(j) = y(j) + a31*K1(j) + a32*K2(j)
 213  100   CONTINUE
 214        CALL FUNC_CHEM( T+alpha3*H, ynew, F1)
 215        nfcn=nfcn+1
 216        hc31 = c31/H
 217        hc32 = c32/H
 218        DO 110 j = 1,NVAR
 219          K3(j) = F1(j) + hc31*K1(j) + hc32*K2(j)
 220  110   CONTINUE
 221        IF (.NOT. Autonomous) THEN
 222          hgam3 = H*gamma3
 223          DO 120 j=1,NVAR
 224             K3(j) = K3(j) + hgam3*DFDT(j)
 225  120     CONTINUE
 226        END IF
 227        CALL KppSolve (JAC, K3)
 228
 229 C ------------ STAGE 4 -------------------------
 230 C              Note: uses the same function value as stage 3
 231        hc41 = c41/H
 232        hc42 = c42/H
 233        hc43 = c43/H
 234        DO 140 j = 1,NVAR
 235          K4(j) = F1(j) + hc41*K1(j) + hc42*K2(j) + hc43*K3(j)
 236  140   CONTINUE
 237        IF (.NOT. Autonomous) THEN
 238           hgam4 = H*gamma4
 239           DO 150 j=1,NVAR
 240             K4(j) = K4(j) + hgam4*DFDT(j)
 241  150      CONTINUE
 242        END IF
 243        CALL KppSolve (JAC, K4)
 244
 245
 246
 247 C ---- The Solution ---
 248        DO 160 j = 1,NVAR
 249          ynew(j) = y(j) + b1*K1(j) + b2*K2(j) + b3*K3(j) + b4*K4(j)
 250  160   CONTINUE
 251
 252
 253 C ====== Error estimation ========
 254
 255         ERR=0.d0
 256         DO 170 j = 1,NVAR
 257            w = AbsTol(j) + RelTol(j)*DMAX1(DABS(y(j)),DABS(ynew(j)))
 258            e = d1*K1(j) + d2*K2(j) + d3*K3(j) + d4*K4(j)
 259            ERR = ERR + ( e/w )**2
 260  170    CONTINUE
 261         ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )
 262
 263 C ======= Choose the stepsize ===============================
 264
 265         elo    = 4.0D0 ! estimator local order
 266         factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
 267         Hnew   = DMIN1(Hmax,DMAX1(Hmin, H/factor))
 268
 269 C ======= Rejected/Accepted Step ============================
 270
 271         IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
 272           IsReject = .true.
 273           H = DMIN1(H/10,Hnew)
 274           Nreject  = Nreject+1
 275         ELSE
 276           DO 180 i=1,NVAR
 277              y(i)  = ynew(i)
 278  180      CONTINUE
 279           T = Tplus
 280           IF (.NOT.IsReject) THEN
 281               H = Hnew   ! Do not increase stepsize if previos step was rejected
 282           END IF
 283           IsReject = .false.
 284           Naccept = Naccept+1
 285         END IF
 286
 287 C ======= End of the time loop ===============================
 288       IF ( T .lt. Tnext ) GO TO 10
 289
 290 C ======= Output Information =================================
 291       Info(2) = Nfcn
 292       Info(3) = Njac
 293       Info(4) = Naccept
 294       Info(5) = Nreject
 295       Hstart  = H
 296
 297       RETURN
 298       END
 299
 300
 301       SUBROUTINE FUNC_CHEM( T, Y, P )
 302       INCLUDE 'KPP_ROOT_params.h'
 303       INCLUDE 'KPP_ROOT_global.h'
 304       KPP_REAL   T, Told
 305       KPP_REAL   Y(NVAR), P(NVAR)
 306       Told = TIME
 307       TIME = T
 308       CALL Update_SUN()
 309       CALL Update_RCONST()
 310       CALL Fun( Y,  FIX, RCONST, P )
 311       TIME = Told
 312       RETURN
 313       END
 314
 315
 316       SUBROUTINE JAC_CHEM( T, Y, J )
 317       INCLUDE 'KPP_ROOT_params.h'
 318       INCLUDE 'KPP_ROOT_global.h'
 319       KPP_REAL   Told, T
 320       KPP_REAL   Y(NVAR), J(LU_NONZERO)
 321       Told = TIME
 322       TIME = T
 323       CALL Update_SUN()
 324       CALL Update_RCONST()
 325       CALL Jac_SP( Y,  FIX, RCONST, J )
 326       TIME = Told
 327       RETURN
 328       END
 329
 330
 331
 332
 333
 334