source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/drv/general_adj.f90 @ 4794

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

Merge of branch palm4u into trunk

File size: 4.0 KB
Line 
1!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2!   Driver for the Adjoint (ADJ) 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 = 1, ind_2 = 2
12      ! INTEGER :: ind_1 = ind_NO2, ind_2 = ind_O3
13
14!~~>  NADJ = Number of functionals for which sensitivities are computed
15!     Note: the setting below is for sensitivities of all final concentrations;
16!           the setting may have to be changed for other applications
17      INTEGER, PARAMETER :: NADJ = NVAR
18      KPP_REAL, DIMENSION(NVAR,NADJ) :: Y_adj, ATOL_adj, RTOL_adj
19
20!~~>  Control (in) and status (out) arguments for the integration
21      KPP_REAL, DIMENSION(20) :: RCNTRL, RSTATUS
22      INTEGER,       DIMENSION(20) :: ICNTRL, ISTATUS
23 
24      STEPMIN = 0.0d0
25      STEPMAX = 0.0d0
26
27!~~~> Tolerances for calculating concentrations       
28      DO i=1,NVAR
29        RTOL(i) = 1.0d-4
30        ATOL(i) = 1.0d-3
31      END DO
32     
33!~~~> Tolerances for calculating adjoints
34!     are used for controlling adjoint truncation error
35!     and for solving the linear adjoint equations by iterations 
36!     Note: Adjoints typically span many orders of magnitude
37!           and a careful tuning of ATOL_adj may be necessary     
38      DO j=1,NADJ
39        DO i=1,NVAR
40          RTOL_adj(i,j) = 1.0d-4
41          ATOL_adj(i,j) = 1.0d-10
42        END DO
43      END DO
44     
45      CALL Initialize()
46     
47!~~~>  The adjoint values at the final time
48      Y_adj(1:NVAR,1:NADJ) = 0.0d0
49      DO j=1,NADJ
50        Y_adj(j,j) = 1.0d0
51      END DO
52
53!~~~> Default control options
54      ICNTRL(1:20) = 0
55      RCNTRL(1:20) = 0.0d0       
56
57!~~~> Begin time loop
58
59      TIME = TSTART
60      CALL InitSaveData()
61
62      T = TSTART
63
64      CALL GetMass( C, DVAL )
65      WRITE(6,991) (T-TSTART)/(TEND-TSTART)*100, T,        &
66                  (TRIM(SPC_NAMES(MONITOR(i))),            &
67                    C(MONITOR(i))/CFACTOR, i=1,NMONITOR),  &
68                  (TRIM(SMASS(i)),DVAL(i)/CFACTOR, i=1,NMASS)
69
70      TIME = T
71      CALL SaveData()
72
73      CALL INTEGRATE_ADJ( NADJ, VAR, Y_adj, T, TEND,       &
74                          ATOL_adj, RTOL_adj,              &
75                          ICNTRL, RCNTRL, ISTATUS, RSTATUS )
76
77
78      CALL GetMass( C, DVAL )
79      WRITE(6,991) (TEND-TSTART)/(TEND-TSTART)*100, TEND,  &
80                  (TRIM(SPC_NAMES(MONITOR(i))),            &
81                    C(MONITOR(i))/CFACTOR, i=1,NMONITOR),  &
82                  (TRIM(SMASS(i)),DVAL(i)/CFACTOR, i=1,NMASS)
83
84      TIME = T
85      CALL SaveData()
86
87!~~~> End time loop ~~~~~~~~~~
88
89      OPEN(20, FILE='KPP_ROOT_ADJ_results.m')
90      WRITE(6,*) '**************************************************'
91      WRITE(6,*) ' Concentrations and Sensitivities at final time '
92      WRITE(6,*) ' were written in the file KPP_ROOT_ADJ_results.m'
93      WRITE(6,*) '**************************************************'
94      DO j=1,NADJ
95        WRITE(20,993) ( Y_adj(i,j), i=1,NVAR )         
96      END DO
97
98      WRITE(6,995) TRIM(SPC_NAMES(ind_1)),TRIM(SPC_NAMES(ind_1)), &
99                   Y_adj(ind_1,1)
100      WRITE(6,995) TRIM(SPC_NAMES(ind_2)),TRIM(SPC_NAMES(ind_2)), &
101                   Y_adj(ind_2,2)
102      WRITE(6,995) TRIM(SPC_NAMES(ind_2)),TRIM(SPC_NAMES(ind_1)), &
103                   Y_adj(ind_1,2)
104      WRITE(6,995) TRIM(SPC_NAMES(ind_1)),TRIM(SPC_NAMES(ind_2)), &
105                   Y_adj(ind_2,1)
106
107      CALL CloseSaveData()
108     
109 991  FORMAT(F6.1,'%. T=',E10.3,3X,20(A,'=',E10.4,';',1X))
110 993  FORMAT(1000(E24.16,2X))
111 995  FORMAT('ADJ: d[',A,'](tf)/d[',A,'](t0)=',E14.7)
112 
113      !~~~> The entire matrix of sensitivities
114      WRITE(6,996) ( 'd ',TRIM(SPC_NAMES(i)), i=1,NVAR ) 
115      DO j=1,NADJ
116        WRITE(6,997) TRIM(SPC_NAMES(j)),( Y_adj(j,i), i=1,NVAR )         
117      END DO
118 996  FORMAT(12X,100('  ',A2,A6,4X))
119 997  FORMAT('d/d',A6,' = ',100(E12.5,2X))
120
121END PROGRAM KPP_ROOT_ADJ_Driver
122
Note: See TracBrowser for help on using the repository browser.