chem/KPP/kpp/kpp-2.1/int/oldies/ros1.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 ROS1(NVAR,TIN,TOUT,STEPMIN,VAR,
  18      +                   Info,FUNC_CHEM,JAC_CHEM)
  19
  20
  21       RETURN
  22       END
  23
  24
  25
  26
  27       SUBROUTINE ROS1(N,T,Tnext,Hstart,
  28      +                   y,Info,FUNC_CHEM,JAC_CHEM)
  29
  30       INCLUDE 'KPP_ROOT_params.h'
  31       INCLUDE 'KPP_ROOT_sparse.h'
  32 C
  33 C  Linearly Implicit Euler
  34 C  A method of theoretical interest but of no practical value
  35 C
  36 C  INPUT ARGUMENTS:
  37 C     y = Vector of (NVAR) concentrations, contains the
  38 C         initial values on input
  39 C     [T, Tnext] = the integration interval
  40 C     Hmin, Hmax = lower and upper bounds for the selected step-size.
  41 C          Note that for Step = Hmin the current computed
  42 C          solution is unconditionally accepted by the error
  43 C          control mechanism.
  44 C     AbsTol, RelTol = (NVAR) dimensional vectors of
  45 C          componentwise absolute and relative tolerances.
  46 C     FUNC_CHEM = name of routine of derivatives. KPP syntax.
  47 C          See the header below.
  48 C     JAC_CHEM = name of routine that computes the Jacobian, in
  49 C          sparse format. KPP syntax. See the header below.
  50 C     Info(1) = 1  for  Autonomous   system
  51 C             = 0  for nonAutonomous system
  52 C
  53 C  OUTPUT ARGUMENTS:
  54 C     y = the values of concentrations at Tend.
  55 C     T = equals TENDon output.
  56 C     Info(2) = # of FUNC_CHEM CALLs.
  57 C     Info(3) = # of JAC_CHEM CALLs.
  58 C     Info(4) = # of accepted steps.
  59 C     Info(5) = # of rejected steps.
  60 C     Hstart = The last accepted stepsize
  61 C
  62 C    Adrian Sandu, December 2001
  63 C
  64       KPP_REAL Fv(NVAR)
  65       KPP_REAL JAC(LU_NONZERO)
  66       KPP_REAL H, Hstart
  67       KPP_REAL y(NVAR)
  68       KPP_REAL T, Tnext, Tplus
  69       KPP_REAL elo,ghinv,uround
  70
  71       INTEGER    n,nfcn,njac,Naccept,Nreject,i,j
  72       INTEGER    Info(5)
  73       LOGICAL    IsReject, Autonomous
  74       EXTERNAL    FUNC_CHEM, JAC_CHEM
  75
  76
  77       H     = Hstart
  78       Tplus = T
  79       Nfcn  = 0
  80       Njac  = 0
  81
  82 C === Starting the time loop ===
  83  10    CONTINUE
  84
  85        Tplus = T + H
  86        IF ( Tplus .gt. Tnext ) THEN
  87           H = Tnext - T
  88           Tplus = Tnext
  89        END IF
  90
  91 C              Initial Function and Jacobian values
  92        CALL FUNC_CHEM(NVAR, T, y, Fv)
  93        Nfcn = Nfcn+1
  94        CALL JAC_CHEM(NVAR, T, y, JAC)
  95        Njac = Njac+1
  96
  97 C              Form the Prediction matrix and compute its LU factorization
  98        DO 40 j=1,NVAR
  99          JAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) - 1.0d0/H
 100  40    CONTINUE
 101        CALL KppDecomp (JAC, ier)
 102 C
 103        IF (ier.ne.0) THEN
 104           PRINT *,'ROS1: Singular factorization at T=',T,'; H=',H
 105           STOP
 106        END IF
 107
 108 C ------------ STAGE 1-------------------------
 109        CALL KppSolve (JAC, Fv)
 110
 111 C ---- The Solution ---
 112        DO 160 j = 1,NVAR
 113          y(j) = y(j) - Fv(j)
 114  160   CONTINUE
 115        T = T + H
 116
 117 C ======= End of the time loop ===============================
 118       IF ( T .lt. Tnext ) GO TO 10
 119
 120 C ======= Output Information =================================
 121       Info(2) = Nfcn
 122       Info(3) = Njac
 123       Info(4) = Njac
 124       Info(5) = 0
 125       Hstart  = H
 126
 127       RETURN
 128       END
 129
 130
 131       SUBROUTINE FUNC_CHEM(N, T, Y, P)
 132       INCLUDE 'KPP_ROOT_params.h'
 133       INCLUDE 'KPP_ROOT_global.h'
 134       INTEGER N
 135       KPP_REAL   T, Told
 136       KPP_REAL   Y(NVAR), P(NVAR)
 137       Told = TIME
 138       TIME = T
 139       CALL Update_SUN()
 140       CALL Update_RCONST()
 141       CALL Fun( Y,  FIX, RCONST, P )
 142       TIME = Told
 143       RETURN
 144       END
 145
 146
 147       SUBROUTINE JAC_CHEM(N, T, Y, J)
 148       INCLUDE 'KPP_ROOT_params.h'
 149       INCLUDE 'KPP_ROOT_global.h'
 150       INTEGER N
 151       KPP_REAL   Told, T
 152       KPP_REAL   Y(NVAR), J(LU_NONZERO)
 153       Told = TIME
 154       TIME = T
 155       CALL Update_SUN()
 156       CALL Update_RCONST()
 157       CALL Jac_SP( Y,  FIX, RCONST, J )
 158       TIME = Told
 159       RETURN
 160       END
 161
 162
 163
 164
 165
 166