source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/drv/general_complete.f90 @ 3225

Last change on this file since 3225 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

File size: 3.7 KB
Line 
1!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2!   Driver for the tangent linear model
3!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4
5PROGRAM KPP_ROOT_ADJ_Driver
6
7  USE KPP_ROOT_Model
8  USE KPP_ROOT_Initialize, ONLY: Initialize
9
10      KPP_REAL :: T, DVAL(NSPEC)
11      INTEGER :: i, j, ind_1 = ind_NO2, ind_2 = ind_O3
12
13! ---  Number of functional for which sensitivities are computed
14! ---  Note: this value is set for sensitivities w.r.t. all initial values
15! ---       it may have to be changed for other applications
16      INTEGER NADJ
17      PARAMETER (NADJ = 2)
18      KPP_REAL  Y_ADJ(NVAR,NADJ)
19      REAL(kind=dble_p)  R1(NVAR), R2(NVAR), V1, V2
20 
21! ---- TIME VARIABLES ------------------     
22
23      STEPMIN = 0.0d0
24      STEPMAX = 0.0d0
25
26      CALL SRAND(89)
27      RTOLS = 1.0d-3
28      DO i=1,NVAR
29        RTOL(i) = RTOLS
30        ATOL(i) = 1.0d-2
31        R1(i) = 10*(RAND()-0.5d0)
32        R2(i) = 10*(RAND()-0.5d0)
33      END DO
34     
35      CALL Initialize()
36! ---  Note: the initial values below are adjoint values at the final time
37     Y_ADJ(1:NVAR,1) = R1(1:NVAR)
38     Y_ADJ(1:NVAR,2) = R2(1:NVAR)
39
40! ********** T LOOP *************************
41
42      CALL InitSaveData()
43
44      WRITE(6,990) (SPC_NAMES(MONITOR(i)), i=1,NMONITOR)
45990   FORMAT('DOne[%] Time[h] ',20(4X,A12))
46
47      T = TSTART
48
49        CALL GetMass( C, DVAL )
50        WRITE(6,991) (T-TSTART)/(TEND-TSTART)*100, T/3600.,   &
51                    (C(MONITOR(i))/CFACTOR, i=1,NMONITOR),  &
52                   (DVAL(i)/CFACTOR, i=1,NMASS)
53991     FORMAT(F6.1,'% ',F7.2,3X,20(E10.4,2X))
54
55        CALL SaveData()
56
57        CALL Update_SUN() 
58        CALL Update_RCONST()
59
60        CALL INTEGRATE_ADJ( NADJ, Y_ADJ, T, TEND )
61
62      V1 = 0.0d0
63      V2 = 0.0d0
64      DO i=1,NVAR
65        V1 = V1 + Y_ADJ(i,1)*R2(i)
66        V2 = V2 + Y_ADJ(i,2)*R1(i)
67      END DO
68
69      PRINT*,'**************************************************'
70      WRITE(6,887) V1
71      WRITE(6,888) V2
72      WRITE(6,889) ABS(V1-V2)/MAX(ABS(V1),ABS(V2))
73887   FORMAT('u.M''*M''.v = ',E24.14 )
74888   FORMAT('v.M''*M''.u = ',E24.14 )
75889   FORMAT('RelativeErr=',E10.3 )
76      PRINT*,'**************************************************'
77
78      CALL GetMass( C, DVAL )
79      WRITE(6,991) (T-TSTART)/(TEND-TSTART)*100, T/3600.,    &
80                    (C(MONITOR(i))/CFACTOR, i=1,NMONITOR), &
81                    (DVAL(i)/CFACTOR, i=1,NMASS)
82
83      CALL SaveData()
84
85! *********** END TIME LOOP ********
86      OPEN(20, FILE='KPP_ROOT_ADJ_results.m')
87      WRITE(6,*) '**************************************************'
88      WRITE(6,*) ' Concentrations and Sensitivities at final time '
89      WRITE(6,*) ' were written in the file KPP_ROOT_ADJ_results.m'
90      WRITE(6,*) '**************************************************'
91      DO j=1,NADJ
92        WRITE(20,993) ( Y_ADJ(i,j), i=1,NVAR )         
93      END DO
94 993  FORMAT(1000(E24.16,2X))
95
96      PRINT*,'ADJ: d[',TRIM(SPC_NAMES(ind_1)),'](tf) / d[', &
97            TRIM(SPC_NAMES(ind_1)),'](t0)=', Y_ADJ(ind_1,1)
98      PRINT*,'ADJ: d[',TRIM(SPC_NAMES(ind_2)),'](tf) / d[', &
99            TRIM(SPC_NAMES(ind_2)),'](t0)=', Y_ADJ(ind_2,2)
100      PRINT*,'ADJ: d[',TRIM(SPC_NAMES(ind_1)),'](tf) / d[', &
101            TRIM(SPC_NAMES(ind_2)),'](t0)=', Y_ADJ(ind_1,2)
102      PRINT*,'ADJ: d[',TRIM(SPC_NAMES(ind_2)),'](tf) / d[', &
103            TRIM(SPC_NAMES(ind_1)),'](t0)=', Y_ADJ(ind_2,1)
104
105      PRINT*,'TLM: d[NO2](tf) / d[NO2](t0)=  1.714961808143527E-002'
106      PRINT*,'TLM: d[ O3](tf) / d[ O3](t0)= -4.447774183920545E-003'
107      PRINT*,'TLM: d[NO2](tf) / d[ O3](t0)=  0.897512294491540'
108      PRINT*,'TLM: d[ O3](tf) / d[NO2](t0)= -5.543729901774693E-005'
109 
110
111      CALL CloseSaveData()
112
113END PROGRAM KPP_ROOT_ADJ_Driver
114
Note: See TracBrowser for help on using the repository browser.