chem/KPP/kpp/kpp-2.1/drv/general_ddm_ic.f

   1 C --- Driver for computing sensitivity coefficients w.r.t. all initial values
   2       PROGRAM driver
   3
   4       INCLUDE 'KPP_ROOT_Parameters.h'
   5       INCLUDE 'KPP_ROOT_Global.h'
   6
   7 C ---  Number of sensitivity coefficients to compute
   8 C ---  Note: this value is set for sensitivities w.r.t. all initial values
   9 C ---       it may have to be changed for other applications
  10       INTEGER NSENSIT
  11       PARAMETER (NSENSIT = NVAR)
  12
  13       KPP_REAL  DVAL(NSPEC)
  14       KPP_REAL  Y(NVAR*(NSENSIT+1))
  15 C ---  Note: Y contains: (1:NVAR) concentrations, followed by
  16 C ---                   (1:NVAR) sensitivities w.r.t. first parameter, followed by
  17 C ---                   etc.,  followed by
  18 C ---                   (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
  19
  20       INTEGER i
  21
  22 C --- The type of sensitivity coefficients to compute
  23 C --- DDMTYPE = 0 : sensitivities w.r.t. initial values
  24 C --- DDMTYPE = 1 : sensitivities w.r.t. parameters
  25       DDMTYPE = 0
  26
  27 C ---- TIME VARIABLES ------------------
  28
  29       TSTART = 0
  30       TEND = TSTART + 600
  31       DT = 60.
  32       TEMP = 298
  33
  34       STEPMIN = 0.01
  35       STEPMAX = 900
  36
  37       PRINT*,'Please provide: RTOL = , ATOL = '
  38       READ*,RTOLS, ATOLS
  39       do i=1,NVAR
  40         RTOL(i) = RTOLS
  41         ATOL(i) = ATOLS
  42       end do
  43
  44 C ********** TIME LOOP *************************
  45
  46       CALL Initialize()
  47
  48 C -- Initialize Concentrations and Sensitivities
  49       DO i=1,NVAR
  50          Y(i) = VAR(i)
  51       END DO
  52
  53 C ---  Note: the initial values below are for sensitivities w.r.t. initial values;
  54 C ---       they may have to be changed for other applications
  55       DO j=1,NSENSIT
  56         DO i=1,NVAR
  57            Y(i+NVAR*j) = 0.0D0
  58         END DO
  59         Y(j+NVAR*j) = 1.0D0
  60       END DO
  61
  62       CALL InitSaveData()
  63
  64       WRITE(6,990) (SPC_NAMES(MONITOR(i)), i=1,NMONITOR),
  65      *             (SMASS(i), i=1,NMASS )
  66 990   FORMAT('done[%] Time[h] ',20(4X,A6))
  67
  68       TIME = TSTART
  69       DO WHILE (TIME .lt. TEND)
  70
  71         CALL GetMass( C, DVAL )
  72         WRITE(6,991) (TIME-TSTART)/(TEND-TSTART)*100, TIME/3600.,
  73      *               (C(MONITOR(i))/CFACTOR, i=1,NMONITOR),
  74      *               (DVAL(i)/CFACTOR, i=1,NMASS)
  75 991   FORMAT(F6.1,'% ',F7.2,3X,20(E10.4,2X))
  76
  77         CALL SaveData()
  78
  79         CALL Update_SUN()
  80         CALL Update_RCONST()
  81
  82         CALL INTEGRATE( NSENSIT, Y, TIME, TIME+DT )
  83
  84         DO i=1,NVAR
  85            VAR(i) = Y(i)
  86         END DO
  87
  88       END DO
  89
  90       CALL GetMass( C, DVAL )
  91       WRITE(6,991) (TIME-TSTART)/(TEND-TSTART)*100, TIME/3600.,
  92      *               (C(MONITOR(i))/CFACTOR, i=1,NMONITOR),
  93      *               (DVAL(i)/CFACTOR, i=1,NMASS)
  94
  95 C      DO i=1,NSENSIT
  96 C        WRITE(6,992) i, ( Y(NVAR*i+j), j=1,NVAR )
  97 C      END DO
  98 C 992  FORMAT('SEN(',I3,') = ',1000(E10.4,2X))
  99
 100       CALL SaveData()
 101
 102 C *********** END TIME LOOP ********
 103       OPEN(20, FILE='KPP_ROOT_results.m')
 104       WRITE(6,*) '**************************************************'
 105       WRITE(6,*) '* Concentrations and Sensitivities at final time *'
 106       WRITE(6,*) '*  were written in the file KPP_ROOT_results.m  *'
 107       WRITE(6,*) '**************************************************'
 108       DO i=0,NSENSIT
 109         WRITE(20,993) ( Y(NVAR*i+j), j=1,NVAR )
 110       END DO
 111  993  FORMAT(1000(E24.16,2X))
 112
 113       CALL CloseSaveData()
 114
 115       STOP
 116       END
 117