forked from KineticPreProcessor/KPP-Standalone
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathkpp_standalone.F90
233 lines (190 loc) · 7.51 KB
/
kpp_standalone.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
program main
! The KPP Standalone for GEOS-Chem Mechanism Analysis
!
!
! Program Description:
!
! This program runs the GEOS-Chem KPP Standalone for a given set of initial conditions.
! It reads an input file generated by the KPP Standalone Interface that generates model
! output of the full chemical state of grid cells in 3D GEOS-Chem, GCHP, and GEOS-CF runs.
! The full mechanism is run to replicate the chemistry of the specified grid cell.
! Obin Sturm ([email protected]), Michael S Long, Christoph Keller
! The KPP Standalone is adapted from a box model used in the following publication:
! Lin, H., Long, M. S., Sander, R., Sandu, A., Yantosca, R. M., Estrada, L. A., et al. (2023).
! An adaptive auto-reduction solver for speeding up integration of chemical kinetics in atmospheric chemistry models:
! Implementation and evaluation in the Kinetic Pre-Processor (KPP) version 3.0.0.
! Journal of Advances in Modeling Earth Systems, 15, e2022MS003293. https://doi.org/10.1029/2022MS003293
! Updates:
! - 2024/05/06, Obin Sturm: Simplification of the code for the GEOS-Chem KPP Standalone,
! removed all autoreduce and convergence criteria testing.
! The more general tool just runs one operator timestep
! and then prints out the results to an output file.
USE GCKPP_GLOBAL
USE GCKPP_JACOBIANSP
USE GCKPP_PARAMETERS
USE GCKPP_MONITOR
USE GCKPP_MODEL
USE KPP_STANDALONE_INIT, ONLY: read_input
IMPLICIT NONE
! Local variables
INTEGER :: ICNTRL(20), IERR, I
INTEGER :: ISTATUS(20)
INTEGER :: fileTotSteps
INTEGER :: level
REAL(dp) :: OperatorTimestep
REAL(dp) :: RCNTRL(20)
REAL(dp) :: Hstart
REAL(dp) :: Hexit
REAL(dp) :: cosSZA
REAL(dp) :: RSTATE(20)
REAL(dp) :: T, TIN, TOUT
REAL(dp) :: start, end
REAL(dp) :: Vloc(NVAR), Cinit(NSPEC), R(NREACT)
REAL :: full_sumtime, full_avg
LOGICAL :: OUTPUT
! Vars for reading files
character(len=256) :: inputfile
character(len=256) :: outputfile
! Check if an argument was provided
if (command_argument_count() .ge. 1) then
! Get the first argument
call get_command_argument(1, inputfile)
print*, 'Processing sample: ', trim(inputfile)
else
print*, 'No sample provided. Exiting.'
stop
endif
! If a second argument is provided, use it as the output file
if (command_argument_count() .ge. 2) then
! Get the second argument
call get_command_argument(2, outputfile)
print*, 'Output file: ', trim(outputfile)
endif
! Initialize
OUTPUT = .false.
RCONST = 0.0_dp
C = 0.0_dp
! Read the input file
call read_input( inputfile, R, Cinit, SPC_NAMES, &
Hstart, Hexit, cosSZA, level, &
fileTotSteps, OperatorTimestep, ICNTRL, RCNTRL, &
ATOL )
! TODO: Pass RTOL from the commmand line
! Run the full mechanism
call fullmech( RTOL_VALUE = 0.5e-2_dp )
! Write the output file
if (command_argument_count() .ge. 2) then
call write_output(inputfile,outputfile)
endif
CONTAINS
! TODO: Pass RTOL as vector
subroutine fullmech(RTOL_VALUE)
USE GCKPP_INTEGRATOR
USE GCKPP_RATES
USE GCKPP_INITIALIZE
USE GCKPP_GLOBAL
IMPLICIT NONE
! Inputs
REAL(dp), INTENT(IN) :: RTOL_VALUE ! Relative tolerance
! Initializze
IERR = 0
ISTATUS = 0
RSTATE = 0.0_dp
! For most integrators, RCNTRL(3) is the starting value of the
! integration step size, so override the initial setting with this.
RCNTRL(3) = Hstart
! Absolute tolerance (ATOL):
! Set to a default value if not defined in the input file.
WHERE( ATOL < 0.0_dp )
ATOL = 1.0e-2_dp
ENDWHERE
! Relative tolerance (RTOL)
RTOL = RTOL_VALUE
! Set ENV
T = 0.0_dp
TIN = T
TOUT = T + OperatorTimestep
! Set initial concentrations (C) and reacton rates (RCONST)
! to values read from the input file
C = Cinit
RCONST = R
! Initialize timings
full_avg = 0.0
full_sumtime = 0.0
start = 0.0_dp
end = 0.0_dp
! Integrate the mechanism for an operator timestep
CALL Integrate( TIN, TOUT, ICNTRL, RCNTRL, ISTATUS, RSTATE, IERR )
! Write results
write( 6, 10 ) fileTotSteps
10 format( " Number of internal timesteps (from 3D run): ", i5 )
write( 6, 11 ) ISTATUS(3)
11 format( " Number of internal timesteps ( standalone): ", i5 )
write( 6, 12 ) Hexit
12 format( " Hexit (from 3D run): ", f10.2 )
write( 6, 13 ) RSTATE(2)
13 format( " Hexit ( standalone): ", f10.2 )
! Check if 3D results are consistent with standalone
if ( fileTotSteps /= ISTATUS(3) ) then
write( 6, 14 )
14 format( "Warning: Number of internal steps do not match 3D grid cell" )
endif
if ( abs( Hexit - RSTATE(2) ) / Hexit > 0.001_dp ) then
write( 6, 15 )
15 format( "Warning: final timestep does not match 3D grid cell within 0.1%" )
endif
end subroutine fullmech
subroutine write_output(inputfile, outputfile)
character(len=256) :: outputfile
character(len=256) :: inputfile
character(len=256) :: header(30)
integer :: i
integer :: ierr
character(len=256) :: line
! Write meteo data lines of the input to the output file
open(20, file=outputfile)
open(10, file=inputfile, status='old')
! Write the header lines to the output file
write(20, '(A)') "43"
write(20, '(A)') REPEAT( "=", 79 )
write(20, '(A)') ""
write(20, '(A)') "KPP Standalone Output"
write(20, '(A)') "This file contains the concentrations of all the chemical species"
write(20, '(A)') "in a single grid cell of a GEOS-Chem 3D run as replicated by the "
write(20, '(A)') "KPP Standalone. Concentrations before and after the operator timestep"
write(20, '(A)') "are in CSV format, below."
write(20, '(A)') ""
write(20, '(A)') "Generated by the GEOS-Chem KPP Standalone:"
write(20, '(A)') "https://github.com/KineticPreProcessor/KPP-Standalone"
write(20, '(A)') ""
write(20, '(A)') "Input file used: " // trim(inputfile)
write(20, '(A)')
! Skip until we get to the metadata section
do
read( 10, '(A)', iostat=ierr ) line
if ( ierr /= 0 ) exit
if ( index( line, 'Meteorological and general' ) > 0 ) then
write(20, '(A)') trim(line)
exit
endif
enddo
! Copy the metadata from the input file to the output file
do
read( 10, '(A)', iostat=ierr ) line
if ( ierr /= 0 ) exit
if ( index( line, 'CSV data of full chemical state' ) > 0 ) exit
write( 20, '(A)' ) trim(line)
enddo
! Close input file
close(10)
! Write initial & final data to output file
write(20, '(A)') REPEAT( "=", 79 )
write(20, '(A)') "Species Name,Initial Concentration (molec/cm3),Final Concentration (molec/cm3)"
! write the species names, initial and final concentrations
do i=1,NSPEC
write(20, '(A, ES25.16E3, A, ES25.16E3)') &
trim(SPC_NAMES(i)) // ",", Cinit(i), ",", C(i)
enddo
close(20)
end subroutine write_output
end program main