chem/KPP/kpp/kpp-2.1/int/oldies/ros2.f90

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