source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/drv/general_soa.f90 @ 3961

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

Merge of branch palm4u into trunk

File size: 3.3 KB
Line 
1!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2!   Driver for the Second Order Adjoint (SOA) model
3!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4
5PROGRAM KPP_ROOT_SOA_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!~~~>  Number of second order adjoints
15!      i.e., number of vectors U_i s.t. SOA_i = Hess*U_i
16!      1 <= i <= NSOA
17
18      INTEGER, PARAMETER :: NSOA = 2
19
20      KPP_REAL :: Y_tlm(NVAR,NSOA)
21      KPP_REAL :: Y_adj(NVAR)
22      KPP_REAL :: Y_soa(NVAR,NSOA)
23 
24      STEPMIN = 0.0d0
25      STEPMAX = 0.0d0
26
27      DO i=1,NVAR
28        RTOL(i) = 1.0d-4
29        ATOL(i) = 1.0d-3
30      END DO
31     
32      CALL Initialize()
33!~~~>  Tangent linear variable values at the initial time
34      Y_tlm(1:NVAR,1:NSOA) = 0.0d0
35      Y_tlm(ind_1,1) = 1.0d0
36      Y_tlm(ind_2,2) = 1.0d0
37!~~~>  Adjoint values at the final time
38      Y_adj(1:NVAR) = 0.0d0
39      Y_adj(ind_1)  = 1.0d0
40!~~~>  2nd order adjoint values at the final time     
41      Y_soa(1:NVAR,1:NSOA) = 0.0d0
42      Y_soa(ind_1,1) = 1.0d0
43      Y_soa(ind_2,2) = 1.0d0
44
45!~~~~~~~~~~~ Time LOOP ~~~~~~~~~~~~~
46
47      CALL InitSaveData()
48
49      T = TSTART
50
51        CALL GetMass( C, DVAL )
52        WRITE(6,991) (T-TSTART)/(TEND-TSTART)*100, T,      &
53                    (TRIM(SPC_NAMES(MONITOR(i))),          &
54                      C(MONITOR(i))/CFACTOR, i=1,NMONITOR),&
55                    (TRIM(SMASS(i)),DVAL(i)/CFACTOR, i=1,NMASS)
56
57        CALL SaveData()
58
59        CALL Update_SUN() 
60        CALL Update_RCONST()
61
62        CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, T, TEND, &
63                            ATOL, RTOL)
64
65
66      CALL GetMass( C, DVAL )
67      WRITE(6,991) (T-TSTART)/(TEND-TSTART)*100, T,        &
68                    (TRIM(SPC_NAMES(MONITOR(i))),          &
69                      C(MONITOR(i))/CFACTOR, i=1,NMONITOR),&
70                    (TRIM(SMASS(i)),DVAL(i)/CFACTOR, i=1,NMASS)
71
72      CALL SaveData()
73
74!~~~~~~~~~~~ END Time LOOP ~~~~~~~~~~~~~
75
76      WRITE(6,*) '**************************************************'
77      WRITE(6,*) ' Results were written in the files'
78      WRITE(6,*) ' KPP_ROOT_[TLM|ADJ|SOA].m'
79      WRITE(6,*) '**************************************************'
80     
81      PRINT 995,TRIM(SPC_NAMES(ind_1)),       &
82            TRIM(SPC_NAMES(ind_1)),           &
83            Y_tlm(ind_1,1), Y_adj(ind_1)
84      PRINT 995,TRIM(SPC_NAMES(ind_1)),       &
85            TRIM(SPC_NAMES(ind_2)),           &
86            Y_tlm(ind_1,2), Y_adj(ind_2)
87
88      DO j=1,NSOA
89         PRINT 997, j,(TRIM(SPC_NAMES(i)),Y_soa(i,j),i=1,NVAR)
90      END DO
91     
92      OPEN(53,FILE='KPP_ROOT_TLM.m')
93      DO j=1, NSOA
94         WRITE(53,993) (Y_tlm(i,j),i=1,NVAR)
95      END DO
96      CLOSE(53)
97     
98      OPEN(54,FILE='KPP_ROOT_ADJ.m')
99      WRITE(54,993) (Y_adj(i),i=1,NVAR)
100      CLOSE(54)
101     
102      OPEN(55,FILE='KPP_ROOT_SOA.m')
103      DO j=1, NSOA
104         WRITE(55,993) (Y_soa(i,j),i=1,NVAR)
105      END DO
106      CLOSE(55)
107     
108 991  FORMAT(F6.1,'%. T=',E10.3,3X,20(A,'=',E10.4,';',1X))
109 993  FORMAT(1000(E24.16,2X))
110 995  FORMAT('d[',A,'](tf)/d[',A,'](t0). TLM=',E14.7,'  ADJ=',E14.7)
111 997  FORMAT('2nd ADJ (',I3,'): ',200(A,'=',E10.3,'; '))
112
113      CALL CloseSaveData()
114
115END PROGRAM KPP_ROOT_SOA_Driver
116
Note: See TracBrowser for help on using the repository browser.