[2696] | 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 | |
---|