wrfv2_fire/chem/KPP/kpp/kpp-2.1/int/oldies/ros3_ddm.f

   1       SUBROUTINE INTEGRATE( NSENSIT, Y, 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 C Y - Concentrations and Sensitivities
  11       KPP_REAL Y(NVAR*(NSENSIT+1))
  12 C ---  Note: Y contains: (1:NVAR) concentrations, followed by
  13 C ---                   (1:NVAR) sensitivities w.r.t. first parameter, followed by
  14 C ---                   etc.,  followed by
  15 C ---                   (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
  16
  17       INTEGER    INFO(5)
  18
  19       EXTERNAL   FUNC_CHEM, JAC_CHEM
  20
  21       INFO(1) = Autonomous
  22
  23         CALL ROS3_DDM(NVAR,NSENSIT,TIN,TOUT,STEPMIN,STEPMAX,
  24      +                   STEPMIN,Y,ATOL,RTOL,
  25      +                   Info,FUNC_CHEM,JAC_CHEM)
  26
  27       RETURN
  28       END
  29
  30
  31       SUBROUTINE ROS3_DDM(N,NSENSIT,T,Tnext,Hmin,Hmax,Hstart,
  32      +                   y,AbsTol,RelTol,
  33      +                   Info,FUNC_CHEM,JAC_CHEM)
  34       IMPLICIT NONE
  35       INCLUDE 'KPP_ROOT_params.h'
  36       INCLUDE 'KPP_ROOT_global.h'
  37       INCLUDE 'KPP_ROOT_sparse.h'
  38
  39 C       L-stable Rosenbrock 3(2), with
  40 C     strongly A-stable embedded formula for error control.
  41 C
  42 C Direct decoupled computation of sensitivities.
  43 C The global variable DDMTYPE distinguishes between:
  44 C      DDMTYPE = 0 : sensitivities w.r.t. initial values
  45 C      DDMTYPE = 1 : sensitivities w.r.t. parameters
  46 C
  47 C     All the arguments aggree with the KPP syntax.
  48 C
  49 C  INPUT ARGUMENTS:
  50 C     y = Vector of:   (1:NVAR) concentrations, followed by
  51 C                      (1:NVAR) sensitivities w.r.t. first parameter, followed by
  52 C                       etc.,  followed by
  53 C                      (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
  54 C         (y contains initial values at input, final values at output)
  55 C     [T, Tnext] = the integration interval
  56 C     Hmin, Hmax = lower and upper bounds for the selected step-size.
  57 C          Note that for Step = Hmin the current computed
  58 C          solution is unconditionally accepted by the error
  59 C          control mechanism.
  60 C     AbsTol, RelTol = (NVAR) dimensional vectors of
  61 C          componentwise absolute and relative tolerances.
  62 C     FUNC_CHEM = name of routine of derivatives. KPP syntax.
  63 C          See the header below.
  64 C     JAC_CHEM = name of routine that computes the Jacobian, in
  65 C          sparse format. KPP syntax. See the header below.
  66 C     Info(1) = 1  for  Autonomous   system
  67 C             = 0  for nonAutonomous system
  68 C
  69 C  OUTPUT ARGUMENTS:
  70 C     y = the values of concentrations at TEND.
  71 C     T = equals TEND on output.
  72 C     Info(2) = # of FUNC_CHEM calls.
  73 C     Info(3) = # of JAC_CHEM calls.
  74 C     Info(4) = # of accepted steps.
  75 C     Info(5) = # of rejected steps.
  76 C
  77 C     Adrian Sandu, April 1996
  78 C     The Center for Global and Regional Environmental Research
  79       INTEGER     NSENSIT
  80       KPP_REAL y(NVAR*(NSENSIT+1)), ynew(NVAR*(NSENSIT+1))
  81       KPP_REAL K1(NVAR*(NSENSIT+1))
  82       KPP_REAL K2(NVAR*(NSENSIT+1))
  83       KPP_REAL K3(NVAR*(NSENSIT+1))
  84       KPP_REAL DFDT(NVAR*(NSENSIT+1))
  85       KPP_REAL DFDP(NVAR*NSENSIT), DFDPDT(NVAR*NSENSIT)
  86       KPP_REAL DJDP(NVAR*NSENSIT)
  87       KPP_REAL JAC(LU_NONZERO), AJAC(LU_NONZERO)
  88       KPP_REAL DJDT(LU_NONZERO)
  89       KPP_REAL Fv(NVAR), Hv(NVAR)
  90       KPP_REAL HESS(NHESS)
  91       KPP_REAL Hmin,Hmax,Hstart,ghinv,uround
  92       KPP_REAL AbsTol(NVAR), RelTol(NVAR)
  93       KPP_REAL T, Tnext, Tplus, H, Hnew, elo
  94       KPP_REAL ERR, factor, facmax
  95       INTEGER    n,nfcn,njac,Naccept,Nreject,i,j,ier
  96       INTEGER    Info(5)
  97       LOGICAL    IsReject,Autonomous
  98       EXTERNAL   FUNC_CHEM, JAC_CHEM, HESS_CHEM
  99
 100       KPP_REAL gamma, c21, c31,c32,b1,b2,b3,d1,d2,d3,a21,a31,a32
 101       KPP_REAL alpha2, alpha3, g1, g2, g3, x1, x2, x3, ytol
 102       KPP_REAL dround, tau
 103
 104       gamma= .43586652150845899941601945119356d+00
 105       c21=  -.10156171083877702091975600115545d+01
 106       c31=   .40759956452537699824805835358067d+01
 107       c32=   .92076794298330791242156818474003d+01
 108        b1=   .10000000000000000000000000000000d+01
 109        b2=   .61697947043828245592553615689730d+01
 110        b3=  -.42772256543218573326238373806514d+00
 111        d1=   .50000000000000000000000000000000d+00
 112        d2=  -.29079558716805469821718236208017d+01
 113        d3=   .22354069897811569627360909276199d+00
 114        a21 = 1.d0
 115        a31 = 1.d0
 116        a32 = 0.d0
 117        alpha2 = gamma
 118        g1=   .43586652150845899941601945119356d+00
 119        g2=   .24291996454816804366592249683314d+00
 120        g3=   .21851380027664058511513169485832d+01
 121
 122
 123 c     Initialization of counters, etc.
 124       Autonomous = Info(1) .EQ. 1
 125       uround = 1.d-15
 126       dround = DSQRT(uround)
 127       IF (Hmax.le.0.D0) THEN
 128           Hmax = DABS(Tnext-T)
 129       END IF
 130       H = DMAX1(1.d-8, Hstart)
 131       Tplus = T
 132       IsReject = .false.
 133       Naccept  = 0
 134       Nreject  = 0
 135       Nfcn     = 0
 136       Njac     = 0
 137
 138
 139 C === Starting the time loop ===
 140  10    CONTINUE
 141
 142 C ====== Initial Function, Jacobian, and Hessian values ===============
 143        CALL FUNC_CHEM(NVAR, T, y, Fv)
 144        Nfcn = Nfcn + 1
 145        CALL JAC_CHEM(NVAR, T, y, JAC)
 146        Njac = Njac + 1
 147        CALL HESS_CHEM( NVAR, T, y, HESS )
 148        IF (DDMTYPE .EQ. 1) THEN
 149           CALL DFUNDPAR(NVAR, NSENSIT, T, y, DFDP)
 150        END IF
 151
 152 C ====== Time derivatives for NONAutonomousous CASE ===============
 153        IF (.not. Autonomous) THEN
 154          tau = DSIGN(dround*DMAX1( 1.0d-6, DABS(T) ), T)
 155          CALL FUNC_CHEM(NVAR, T+tau, y, K2)
 156          nfcn=nfcn+1
 157          DO 20 j = 1,NVAR
 158            DFDT(j) = ( K2(j)-Fv(j) )/tau
 159  20      CONTINUE
 160          CALL JAC_CHEM(NVAR, T+tau, y, AJAC)
 161          DO 30 j = 1,LU_NONZERO
 162            DJDT(j) = ( AJAC(j)-JAC(j) )/tau
 163  30      CONTINUE
 164          DO 40 i=1,NSENSIT
 165             CALL Jac_SP_Vec (DJDT,y(i*NVAR+1),DFDT(i*NVAR+1))
 166  40      CONTINUE
 167        END IF
 168
 169
 170        Tplus = T + H
 171        IF ( Tplus .gt. Tnext ) then
 172           H = Tnext - T
 173           Tplus = Tnext
 174        END IF
 175
 176        gHinv = 1.0d0/(gamma*H)
 177        DO 50 j=1,LU_NONZERO
 178          AJAC(j) = -JAC(j)
 179  50    CONTINUE
 180        DO 60 j=1,NVAR
 181          AJAC(LU_DIAG(j)) = AJAC(LU_DIAG(j)) + gHinv
 182  60    CONTINUE
 183        CALL KppDecomp (AJAC, ier)
 184
 185        IF (ier.NE.0) THEN
 186          IF ( H.GT.Hmin) THEN
 187             H = 5.0d-1*H
 188             GO TO 10
 189          ELSE
 190             PRINT *,'IER <> 0, H=',H
 191             STOP
 192          END IF
 193        END IF
 194
 195        Autonomous = .true.
 196
 197 C ------------------------------- STAGE 1 --------------------------------------
 198        DO 70 j = 1,NVAR
 199            K1(j) =  Fv(j)
 200  70    CONTINUE
 201        IF (.NOT.Autonomous) THEN
 202          x1 = gamma*H
 203          DO 80 j = 1,NVAR
 204             K1(j) = K1(j) + x1*DFDT(j)
 205  80     CONTINUE
 206        END IF
 207        CALL KppSolve (AJAC, K1)
 208 C               --- If  derivative w.r.t. parameters
 209          IF (DDMTYPE .EQ. 1) THEN
 210             CALL DJACDPAR(NVAR, NSENSIT, T, y, K1(1), DJDP)
 211          END IF
 212 C               --- End of derivative w.r.t. parameters
 213
 214        DO 110 i=1,NSENSIT
 215           CALL Jac_SP_Vec (JAC,y(i*NVAR+1),K1(i*NVAR+1))
 216           CALL Hess_Vec ( HESS, y(i*NVAR+1), K1(1), Hv )
 217           DO 90 j=1,NVAR
 218             K1(i*NVAR+j) = K1(i*NVAR+j) + Hv(j)
 219  90       CONTINUE
 220           IF (.NOT. Autonomous) THEN
 221             DO 100 j=1,NVAR
 222               K1(i*NVAR+j) = K1(i*NVAR+j) + x1*DFDT(i*NVAR+j)
 223  100        CONTINUE
 224           END IF
 225 C               --- If  derivative w.r.t. parameters
 226           IF (DDMTYPE .EQ. 1) THEN
 227             DO 44 j = 1,NVAR
 228                K1(i*NVAR+j) = K1(i*NVAR+j) + DFDP((i-1)*NVAR+j)
 229      &                           + DJDP((i-1)*NVAR+j)
 230  44        CONTINUE
 231           END IF
 232 C               --- End of derivative w.r.t. parameters
 233           CALL KppSolve (AJAC, K1(i*NVAR+1))
 234  110   CONTINUE
 235
 236 C ------------------------------- STAGE 2 --------------------------------------
 237        DO 120 j = 1,NVAR*(NSENSIT+1)
 238          ynew(j) = y(j) + a21*K1(j)
 239  120   CONTINUE
 240        CALL FUNC_CHEM(NVAR, T + alpha2*H, ynew, Fv)
 241        IF (DDMTYPE .EQ. 1) THEN
 242          CALL DFUNDPAR(NVAR, NSENSIT, T+alpha3*H, ynew, DFDP)
 243        END IF
 244        nfcn=nfcn+1
 245        x1 = c21/H
 246        DO 130 j = 1,NVAR
 247            K2(j) = Fv(j) + x1*K1(j)
 248  130   CONTINUE
 249        IF (.NOT.Autonomous) THEN
 250          x2 = g2*H
 251          DO 140 j = 1,NVAR
 252             K2(j) = K2(j) + x2*DFDT(j)
 253  140     CONTINUE
 254        END IF
 255        CALL KppSolve (AJAC, K2)
 256 C               --- If  derivative w.r.t. parameters
 257        IF (DDMTYPE .EQ. 1) THEN
 258           CALL DJACDPAR(NVAR, NSENSIT, T, y, K2(1), DJDP)
 259        END IF
 260 C               --- End of derivative w.r.t. parameters
 261
 262        CALL JAC_CHEM(NVAR, T+alpha2*H, ynew, JAC)
 263        njac=njac+1
 264        DO 170 i=1,NSENSIT
 265           CALL Jac_SP_Vec (JAC,ynew(i*NVAR+1),K2(i*NVAR+1))
 266           CALL Hess_Vec ( HESS, y(i*NVAR+1), K2(1), Hv )
 267           DO 150 j = 1,NVAR
 268              K2(i*NVAR+j) = K2(i*NVAR+j) + x1*K1(i*NVAR+j)
 269      &                          + Hv(j)
 270  150      CONTINUE
 271           IF (.NOT. Autonomous) THEN
 272              DO 160 j=1,NVAR
 273                 K2(i*NVAR+j) = K2(i*NVAR+j) + x2*DFDT(i*NVAR+j)
 274  160         CONTINUE
 275           END IF
 276 C               --- If  derivative w.r.t. parameters
 277           IF (DDMTYPE .EQ. 1) THEN
 278             DO 165 j = 1,NVAR
 279                K2(i*NVAR+j) = K2(i*NVAR+j) + DFDP((i-1)*NVAR+j)
 280      &                           + DJDP((i-1)*NVAR+j)
 281  165        CONTINUE
 282           END IF
 283 C               --- End of derivative w.r.t. parameters
 284           CALL KppSolve (AJAC, K2(i*NVAR+1))
 285  170   CONTINUE
 286
 287 C ------------------------------- STAGE 3  --------------------------------------
 288        x1 = c31/H
 289        x2 = c32/H
 290        DO 180 j = 1,NVAR
 291          K3(j) = Fv(j) + x1*K1(j) + x2*K2(j)
 292  180   CONTINUE
 293        IF (.NOT.Autonomous) THEN
 294          x3 = g3*H
 295          DO 190 j = 1,NVAR
 296             K3(j) = K3(j) + x3*DFDT(j)
 297  190     CONTINUE
 298        END IF
 299        CALL KppSolve (AJAC, K3)
 300 C               --- If  derivative w.r.t. parameters
 301          IF (DDMTYPE .EQ. 1) THEN
 302             CALL DJACDPAR(NVAR, NSENSIT, T, y, K3(1), DJDP)
 303          END IF
 304 C               --- End of derivative w.r.t. parameters
 305
 306        DO 220 i=1,NSENSIT
 307           CALL Jac_SP_Vec (JAC,ynew(i*NVAR+1),K3(i*NVAR+1))
 308           CALL Hess_Vec ( HESS, y(i*NVAR+1), K3(1), Hv )
 309           DO 200 j = 1,NVAR
 310              K3(i*NVAR+j) = K3(i*NVAR+j) +x1*K1(i*NVAR+j)
 311      &                       + x2*K2(i*NVAR+j) + Hv(j)
 312  200      CONTINUE
 313           IF (.NOT. Autonomous) THEN
 314              DO 210 j=1,NVAR
 315                 K3(i*NVAR+j) = K3(i*NVAR+j) + x3*DFDT(i*NVAR+j)
 316  210         CONTINUE
 317           END IF
 318 C               --- If  derivative w.r.t. parameters
 319           IF (DDMTYPE .EQ. 1) THEN
 320              DO 215 j = 1,NVAR
 321                K3(i*NVAR+j) = K3(i*NVAR+j) + DFDP((i-1)*NVAR+j)
 322      &                           + DJDP((i-1)*NVAR+j)
 323  215        CONTINUE
 324           END IF
 325 C               --- End of derivative w.r.t. parameters
 326           CALL KppSolve (AJAC, K3(i*NVAR+1))
 327  220   CONTINUE
 328
 329 C ------------------------------ The Solution ---
 330
 331        DO 230 j = 1,NVAR*(NSENSIT+1)
 332          ynew(j) = y(j) + b1*K1(j) + b2*K2(j) + b3*K3(j)
 333  230   CONTINUE
 334
 335
 336 C ====== Error estimation ========
 337
 338         ERR=0.d0
 339         DO 240 i=1,NVAR
 340            ytol = AbsTol(i) + RelTol(i)*DABS(ynew(i))
 341            ERR=ERR+((d1*K1(i)+d2*K2(i)+d3*K3(i))/ytol)**2
 342  240    CONTINUE
 343         ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )
 344
 345 C ======= Choose the stepsize ===============================
 346
 347         elo    = 3.0D0 ! estimator local order
 348         factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
 349         Hnew   = DMIN1(Hmax,DMAX1(Hmin, H/factor))
 350
 351 C ======= Rejected/Accepted Step ============================
 352
 353         IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
 354           IsReject = .true.
 355           H = DMIN1(H/10,Hnew)
 356           Nreject  = Nreject+1
 357           GO TO 10
 358         ELSE
 359           DO 250 j=1,NVAR*(NSENSIT+1)
 360              y(j)  = ynew(j)
 361  250      CONTINUE
 362           T = Tplus
 363           IF (.NOT.IsReject) THEN
 364               H = Hnew   ! Do not increase stepsize IF previos step was rejected
 365           END IF
 366           IsReject = .false.
 367           Naccept = Naccept+1
 368         END IF
 369
 370
 371 C ======= End of the time loop ===============================
 372       IF ( T .lt. Tnext ) GO TO 10
 373
 374
 375 C ======= Output Information =================================
 376       Info(2) = Nfcn
 377       Info(3) = Njac
 378       Info(4) = Naccept
 379       Info(5) = Nreject
 380       Hstart  = H
 381
 382       RETURN
 383       END
 384
 385
 386
 387       SUBROUTINE FUNC_CHEM(N, T, Y, P)
 388       INCLUDE 'KPP_ROOT_params.h'
 389       INCLUDE 'KPP_ROOT_global.h'
 390       INTEGER N
 391       KPP_REAL   T, Told
 392       KPP_REAL   Y(NVAR), P(NVAR)
 393       Told = TIME
 394       TIME = T
 395       CALL Update_SUN()
 396       CALL Update_RCONST()
 397       CALL Fun( Y,  FIX, RCONST, P )
 398       TIME = Told
 399       RETURN
 400       END
 401
 402
 403       SUBROUTINE DFUNDPAR(N, NSENSIT, T, Y, P)
 404 C ---  Computes the partial derivatives of FUNC_CHEM w.r.t. parameters
 405       INCLUDE 'KPP_ROOT_params.h'
 406       INCLUDE 'KPP_ROOT_global.h'
 407 C ---  NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients
 408       INTEGER N
 409       INTEGER NCOEFF, JCOEFF(NREACT)
 410       COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
 411
 412       KPP_REAL   T, Told
 413       KPP_REAL   Y(NVAR), P(NVAR*NSENSIT)
 414       Told = TIME
 415       TIME = T
 416       CALL Update_SUN()
 417       CALL Update_RCONST()
 418 C
 419       IF (DDMTYPE .EQ. 0) THEN
 420 C ---  Note: the values below are for sensitivities w.r.t. initial values;
 421 C ---       they may have to be changed for other applications
 422         DO j=1,NSENSIT
 423           DO i=1,NVAR
 424             P(i+NVAR*(j-1)) = 0.0D0
 425           END DO
 426         END DO
 427       ELSE
 428 C ---  Example: the call below is for sensitivities w.r.t. rate coefficients;
 429 C ---       JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
 430 C ---       w.r.t. which one differentiates
 431         CALL dFun_dRcoeff( Y,  FIX, NCOEFF, JCOEFF, P )
 432       END IF
 433       TIME = Told
 434       RETURN
 435       END
 436
 437       SUBROUTINE JAC_CHEM(N, T, Y, J)
 438       INCLUDE 'KPP_ROOT_params.h'
 439       INCLUDE 'KPP_ROOT_global.h'
 440       INTEGER N
 441       KPP_REAL   Told, T
 442       KPP_REAL   Y(NVAR), J(LU_NONZERO)
 443       Told = TIME
 444       TIME = T
 445       CALL Update_SUN()
 446       CALL Update_RCONST()
 447       CALL Jac_SP( Y,  FIX, RCONST, J )
 448       TIME = Told
 449       RETURN
 450       END
 451
 452
 453       SUBROUTINE DJACDPAR(N, NSENSIT, T, Y, U, P)
 454 C ---  Computes the partial derivatives of JAC w.r.t. parameters times user vector U
 455       INCLUDE 'KPP_ROOT_params.h'
 456       INCLUDE 'KPP_ROOT_global.h'
 457 C ---  NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients
 458       INTEGER N
 459       INTEGER NCOEFF, JCOEFF(NREACT)
 460       COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
 461
 462       KPP_REAL   T, Told
 463       KPP_REAL   Y(NVAR), U(NVAR)
 464       KPP_REAL   P(NVAR*NSENSIT)
 465       Told = TIME
 466       TIME = T
 467       CALL Update_SUN()
 468       CALL Update_RCONST()
 469 C
 470       IF (DDMTYPE .EQ. 0) THEN
 471 C ---  Note: the values below are for sensitivities w.r.t. initial values;
 472 C ---       they may have to be changed for other applications
 473         DO j=1,NSENSIT
 474           DO i=1,NVAR
 475             P(i+NVAR*(j-1)) = 0.0D0
 476           END DO
 477         END DO
 478       ELSE
 479 C ---  Example: the call below is for sensitivities w.r.t. rate coefficients;
 480 C ---       JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
 481 C ---       w.r.t. which one differentiates
 482         CALL dJac_dRcoeff( Y,  FIX, U, NCOEFF, JCOEFF, P )
 483       END IF
 484       TIME = Told
 485       RETURN
 486       END
 487
 488
 489       SUBROUTINE HESS_CHEM(N, T, Y, HESS)
 490       INCLUDE 'KPP_ROOT_params.h'
 491       INCLUDE 'KPP_ROOT_global.h'
 492       INTEGER N
 493       KPP_REAL   Told, T
 494       KPP_REAL   Y(NVAR), HESS(NHESS)
 495       Told = TIME
 496       TIME = T
 497       CALL Update_SUN()
 498       CALL Update_RCONST()
 499       CALL Hessian( Y,  FIX, RCONST, HESS )
 500       TIME = Told
 501       RETURN
 502       END
 503