Home

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/create_kpp_module.C @ 3458

Last change on this file since 3458 was 3458, checked in by kanani, 5 years ago
Reintegrated fixes/changes from branch chemistry
File size: 22.4 KB

Line
1
2	// ############################################################################
3	//
4	// create_kpp_module
5	//
6	// create scalar code from .f90 sources created by KPP
7	//
8	// COPYRIGHT Klaus Ketelsen and MPI-CH April 2007
9	//
10	// ############################################################################
11	//
12	//Current revisions:
13	//------------------
14	//
15	//
16	//Former revisions:
17	//-----------------
18	//$Id: create_kpp_module.C 3453 2018-10-30 13:21:51Z forkel $
19	// exclude kco_compress from handling by global_variables2vector (30.10.2018, forkel)
20	//
21	// Added automatic line with mechanism name (read mech_list) (25.09.2018, forkel)
22	//
23	// Added vl_glo = size(tempi,1) (20.09.2018, forkel)
24	//
25	// Removed creation of fill_ Subroutine and creation of calls thereof (18.09.2018, ketelsen)
26	//
27	// Fix in order not to loose the values of qvap and fakt (12.09.2018, forkel)
28	//
29	// Bug fixes: moved kppi.add_line(" CALL initialize after fakt = fakti(is)
30	// Deleted definition of qvap,fakt in create_kpp_integrate again (03.09.2018, forkel)
31	//
32	// Changes for vector mode (edit_WAXPY, edit_FunTemplate, edit_JacTemplate,
33	// some cleanup of comments, various changes in create_kpp_integrate) (July 2018, ketelsen)
34	//
35	//
36	// Added qvap and fakt (June 2018, forkel)
37	// --> Change in module_header: qvap, fakt added (June 2018, forkel)
38	//
39	// re-established original uppercase/lowercase (June 2018, forkel)
40	// --> Change in module_header: reset case in Initialize, Integrate, and
41	// Update_rconst (June 2018, forkel)
42	//
43	// Removed preprocessor directive __chem again (2017-09-14, forkel)
44	//
45	// Added phot (2017-09-14, forkel)
46	// --> Change in module_header: Variables for photolyis added (2017-09-14, forkel)
47	//
48	// change of some output to lowercase with uppercase Fortran (2017, forkel)
49	//
50	// Intial version of KP4 adapted to PALM (Nov. 2016, ketelsen)
51	//
52	#include <stdio.h>
53	// stdlib is necessary to define getenv:
54	#include <stdlib.h>
55
56	#include "create_kpp_module.h"
57	#include "utils.h"
58
59	void create_kpp_module::do_work (string s) {
60	vector<fortran_file>::iterator it;
61	vector<string>::iterator ic;
62	vector<Vvar>::iterator iv;
63
64	expand_decomp exp_de;
65
66	prefix = s;
67	module_name = prefix;
68
69	cout << "Create " << module_name << " from kpp Fortran sources" <<endl;
70	cout << "Vector mode " << kpp_switches.is_vector() <<endl;
71	cout << "De_indexing " << kpp_switches.de_indexing() <<endl;
72
73	create_fortran_files_and_read();
74
75	// Generate first module lines
76
77	string first_line="MODULE " + module_name;
78	mz_kpp.add_line(first_line);
79	mz_kpp.add_line(" ");
80
81	// string e5_line = first_line +"_e5";
82	// e5_kpp.add_line(e5_line);
83	// e5_line = " USE " + module_name;
84	// e5_kpp.add_line(e5_line);
85	// e5_kpp.add_line(" ");
86
87	// edit include files
88
89	for(it=kpp_includes.begin();it!=kpp_includes.end();it++) {
90	it->edit_inc(header_variables);
91
92	// Create variable Species list and vector variable list
93
94	if(it->get_name() == module_name + "_Parameters") {
95	it->create_species_list(species_list);
96	}
97	if(it->get_name() == module_name + "_Global") {
98	it->vector_variable_list(Vvar_list);
99	}
100	}
101
102	// Prepare expansion of decomposition subroutine
103
104	if(kpp_switches.de_indexing () > 0 ) {
105	exp_de.create_sparse_info (kpp_includes, module_name);
106	}
107
108	// edit FORTRAN files
109
110	for(it=kpp_files.begin();it!=kpp_files.end();it++) {
111	it->edit_fortran ();
112	}
113
114	// Generate a list of single subroutines from kpp-files
115	// kpp files are modules containing several subroutines
116
117	copy_files_to_subroutines ();
118
119	// All header_variables to include list
120
121	kpp_includes.push_back(header_variables);
122
123	// Create decomposition subroutine
124	if(kpp_switches.de_indexing () > 0 ) {
125	exp_de.create_routine (kpp_subroutines);
126	}
127
128	if(kpp_switches.is_vector()) {
129
130	cout << "##### Hier kpp_switches.is_vector " <<endl;
131	// Change header section
132	for(it=kpp_includes.begin();it!=kpp_includes.end();it++) {
133	it->edit_inc_vec(global_variable_list);
134	}
135
136	// Change global variables to vector (except for kp4_compress, which has already the right form)
137
138	for(it=kpp_subroutines.begin();it!=kpp_subroutines.end();it++) {
139	if(it->get_name() != "kco_compress" ) {
140	it->global_variables2vector (global_variable_list);
141	}
142	}
143
144	// Edit individual subroutines
145
146	for(it=kpp_subroutines.begin();it!=kpp_subroutines.end();it++) {
147	if(it->get_name() == "KppDecomp") {
148	it->edit_KppDecomp();
149	}
150	if(it->get_name() == "KppSolve") {
151	it->edit_KppSolve();
152	}
153	if(it->get_name() == "Jac_SP" ) {
154	it->edit_Jac_SP();
155	}
156	if(it->get_name() == "Fun" ) {
157	it->edit_Fun();
158	}
159	if(it->get_name() == "WAXPY" ) {
160	it->edit_WAXPY();
161	}
162	if(it->get_name() == "FunTemplate" ) {
163	it->edit_FunTemplate();
164	}
165	if(it->get_name() == "JacTemplate" ) {
166	it->edit_JacTemplate();
167	}
168	}
169	}
170
171	// Update_RCONST has to be changed also in scalar mode
172
173	for(it=kpp_subroutines.begin();it!=kpp_subroutines.end();it++) {
174	if(it->get_name() == "Update_RCONST") {
175	it->edit_Update_RCONST(Vvar_list);
176	}
177
178	if(it->get_name() == "Initialize") {
179	it->edit_Initialize(Vvar_list);
180	}
181
182	}
183
184	// Add Solver template to subroutine list
185	if(kpp_switches.is_vector()) {
186	add_solver_to_subroutine_list ();
187	}
188
189	// The module header will be taken from ../templates/module_header.
190	// Please edit if header has to be changed.
191
192	generate_module_header();
193
194	// Create kpp_integrate subroutine (chem_gasphase_integrate) for skalar and vector mode
195
196	create_kpp_integrate();
197	// Copy include files
198
199	for(it=kpp_includes.begin();it!=kpp_includes.end();it++) {
200	it->copy_to_MZ_KPP(mz_kpp);
201	}
202
203	mz_kpp.add_line(" ");
204	mz_kpp.add_line("! Interface Block ");
205	mz_kpp.add_line(" ");
206	for(it=kpp_subroutines.begin();it!=kpp_subroutines.end();it++) {
207	string buf;
208
209	string prefix = " ";
210	for(ic=interface_ignore.begin();ic!=interface_ignore.end();ic++) {
211	if(it->get_name() == *ic) {
212	prefix = "!interface not working ";
213	break;
214	}
215	}
216
217	buf = prefix + "interface " + it->get_name() ;
218	mz_kpp.add_line(buf);
219	buf = prefix + " module procedure " + it->get_name();
220	mz_kpp.add_line(buf);
221	buf = prefix + "end interface " + it->get_name();
222	mz_kpp.add_line(buf);
223	mz_kpp.add_line(" ");
224	}
225
226	mz_kpp.add_line(" ");
227
228
229	// Copy FORTRAN subroutines to mz_kpp
230
231	mz_kpp.add_line(" CONTAINS");
232
233	for(it=kpp_subroutines.begin();it!=kpp_subroutines.end();it++) {
234	mz_kpp.add_line(" ");
235	it->copy_to_MZ_KPP(mz_kpp);
236	}
237
238	// Finish module
239
240	string last_line="end module " + module_name;
241	mz_kpp.add_line("");
242	mz_kpp.add_line(last_line);
243
244	// Write the complete module to file: mz_kpp.f
245
246	write_module_file();
247
248	return;
249	}
250
251	void create_kpp_module::create_fortran_files_and_read() {
252
253	string name;
254	ifstream in,in_c,in_b,in_i;
255	fortran_file f_file;
256	vector<fortran_file>::iterator it;
257
258	// Open file with list of FORTRAN routines
259
260	in.open("file_list");
261	if( !in ) {
262	cout << "cannot open " << endl; my_abort("file_list");
263	}
264
265	// Create kpp_fortran routines
266	while ( 1 ) {
267	in >> name;
268	if( in.eof() ) break;
269	if( in.bad() ) my_abort("ERROR_READ_1");
270	f_file.set_name(name);
271	kpp_files.push_back(f_file);
272	}
273	in.close();
274
275	// Read FORTRAN code
276
277	for(it=kpp_files.begin();it!=kpp_files.end();it++) {
278	it->read();
279	}
280
281	// Open file with list of include files
282
283	in_c.open("include_list");
284	if( !in_c ) {
285	cout << "cannot open " << endl; my_abort("include_list");
286	}
287
288	// Create kpp_includes vector
289	while ( 1 ) {
290	in_c >> name;
291	if( in_c.eof() ) break;
292	if( in_c.bad() ) my_abort("ERROR_READ_3");
293	f_file.set_name(name);
294	kpp_includes.push_back(f_file);
295	}
296	in_c.close();
297
298	// Read include files
299
300	for(it=kpp_includes.begin();it!=kpp_includes.end();it++) {
301	it->read();
302	}
303
304	// Read Ignore list
305
306	in_i.open("interface_ignore_list");
307	if( !in_i ) {
308	cout << "cannot open " << endl; my_abort("include_list");
309	}
310
311	// Create kpp_includes vector
312	while ( 1 ) {
313	in_i >> name;
314	if( in_i.eof() ) break;
315	if( in_i.bad() ) my_abort("ERROR_READ_4");
316	interface_ignore.push_back(name);
317	}
318	in_c.close();
319
320	}
321
322	void create_kpp_module::copy_files_to_subroutines () {
323	string name;
324	ifstream in;
325	fortran_file s_file;
326	vector<fortran_file>::iterator it;
327
328	// Open file with list of FORTRAN routines
329
330	in.open("subroutine_list");
331	if( !in ) {
332	cout << "cannot open " << endl; my_abort("subroutine_list");
333	}
334
335	// Create vector kpp_subroutines
336
337	while ( 1 ) {
338	in >> name;
339	if( in.eof() ) break;
340	if( in.bad() ) my_abort("ERROR_READ_S1");
341	s_file.set_name(name);
342	kpp_subroutines.push_back(s_file);
343	}
344	in.close();
345
346	header_variables.add_line(" ");
347	header_variables.add_line("! variable definations from individual module headers ");
348	header_variables.add_line(" ");
349
350	// Loop over all FORTRAN Files
351
352	for(it=kpp_files.begin();it!=kpp_files.end();it++) {
353	it->copy_to_subroutine_vector(kpp_subroutines, header_variables);
354	}
355	}
356
357	void create_kpp_module::add_solver_to_subroutine_list () {
358	fortran_file s_file;
359
360	string solver_name = getenv("KPP_SOLVER");
361	cout << "KPP_SOLVER " <<solver_name <<endl;
362
363	s_file.set_name(solver_name);
364	s_file.read();
365	kpp_subroutines.push_back(s_file);
366
367	return;
368	}
369
370	void create_kpp_module::generate_module_header() {
371
372	string buf;
373	ifstream in;
374	ifstream in_e5;
375	program_line line;
376	vector<fortran_file>::iterator it;
377	char distr[2];
378	string diline;
379
380	// Read mechanism from mech_list
381
382	in.open("mech_list");
383	if( !in ) {
384	cout << "cannot open " << endl; my_abort("mech_list");
385	}
386
387	while ( 1 ) {
388	getline (in, buf);
389	if( in.eof() ) break;
390	if( in.bad() ) my_abort("ERROR_READ_4");
391	line.set_line(buf);
392	mz_kpp.add_line(line);
393	}
394	in.close();
395
396
397	// Read Modul Header from file $MZ_KPP_HOME/templates/module_header
398
399	in.open("module_header");
400	if( !in ) {
401	cout << "cannot open " << endl; my_abort("module_header");
402	}
403
404	while ( 1 ) {
405	getline (in, buf);
406	if( in.eof() ) break;
407	if( in.bad() ) my_abort("ERROR_READ_4");
408	line.set_line(buf);
409	mz_kpp.add_line(line);
410	}
411	mz_kpp.add_line(" ");
412	mz_kpp.add_line("! Variables used for vector mode ");
413	mz_kpp.add_line(" ");
414	if(kpp_switches.is_vector()) {
415	mz_kpp.add_line(" logical,parameter :: L_VECTOR = .TRUE. ");
416	} else {
417	mz_kpp.add_line(" logical,parameter :: L_VECTOR = .FALSE. ");
418	}
419	// mz_pj_20070531+
420	sprintf(distr,"%i",kpp_switches.de_indexing());
421	diline = distr ;
422	mz_kpp.add_line(" integer,parameter :: I_LU_DI = " + diline );
423	// mz_pj_20070531-
424
425	mz_kpp.add_line(" integer,parameter :: VL_DIM = "
426	+ kpp_switches.get_vector_length() );
427	mz_kpp.add_line(" integer :: vl ");
428	mz_kpp.add_line(" ");
429	mz_kpp.add_line(" integer :: VL_glo ");
430	mz_kpp.add_line(" integer :: is,ie ");
431	mz_kpp.add_line(" ");
432	mz_kpp.add_line(" ");
433	mz_kpp.add_line(" integer, dimension(VL_dim) :: Kacc,Krej ");
434	mz_kpp.add_line(" integer, dimension(VL_dim) :: IERRV ");
435	mz_kpp.add_line(" logical :: data_loaded = .false. ");
436
437	in.close();
438
439	return;
440	}
441
442	void create_kpp_module::write_module_file() {
443	ofstream out;
444	ofstream out_e5;
445
446	string out_file = "kk_kpp.f90";
447	out.open(out_file.c_str(), ios::out);
448	if( !out ) {
449	cout << "cannot open " << endl; my_abort(out_file);
450	}
451
452	mz_kpp.write_file (out);
453
454	out.close();
455
456
457	return;
458	}
459
460	void create_kpp_module::create_kpp_integrate() {
461	fortran_file kppi;
462	vector<Vvar>::iterator iv;
463	string line;
464
465	kppi.set_name("chem_gasphase_integrate");
466
467	kppi.add_line("SUBROUTINE chem_gasphase_integrate (time_step_len, conc, tempi, qvapi, fakti, photo, ierrf, xnacc, xnrej, istatus, l_debug, pe, icntrl_i, rcntrl_i ) ");
468	kppi.add_line(" ");
469	kppi.add_line(" IMPLICIT NONE ");
470	kppi.add_line(" ");
471
472	kppi.add_line(" REAL(dp), INTENT(IN) :: time_step_len ");
473	kppi.add_line(" REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: conc ");
474	kppi.add_line(" REAL(dp), DIMENSION(:,:), INTENT(IN) :: photo ");
475	kppi.add_line(" REAL(dp), DIMENSION(:), INTENT(IN) :: tempi ");
476	kppi.add_line(" REAL(dp), DIMENSION(:), INTENT(IN) :: qvapi ");
477	kppi.add_line(" REAL(dp), DIMENSION(:), INTENT(IN) :: fakti ");
478	kppi.add_line(" INTEGER, INTENT(OUT), OPTIONAL :: ierrf(:) ");
479	kppi.add_line(" INTEGER, INTENT(OUT), OPTIONAL :: xNacc(:) ");
480	kppi.add_line(" INTEGER, INTENT(OUT), OPTIONAL :: xNrej(:) ");
481	kppi.add_line(" INTEGER, INTENT(INOUT), OPTIONAL :: istatus(:) ");
482	kppi.add_line(" INTEGER, INTENT(IN), OPTIONAL :: PE ");
483	kppi.add_line(" LOGICAL, INTENT(IN), OPTIONAL :: l_debug ");
484	kppi.add_line(" INTEGER, DIMENSION(nkppctrl),INTENT(IN), OPTIONAL :: icntrl_i ");
485	kppi.add_line(" REAL(dp), DIMENSION(nkppctrl),INTENT(IN), OPTIONAL :: rcntrl_i ");
486	kppi.add_line(" ");
487	kppi.add_line(" INTEGER :: k ! loop variable ");
488	kppi.add_line(" REAL(dp) :: dt ");
489	kppi.add_line(" integer, dimension(20) :: istatus_u ");
490	kppi.add_line(" integer :: ierr_u ");
491	kppi.add_line(" integer :: istatf ");
492	kppi.add_line(" integer :: vl_dim_lo ");
493	kppi.add_line(" ");
494	kppi.add_line(" ");
495	kppi.add_line(" if (present (istatus) ) istatus = 0 ");
496	kppi.add_line(" if (present (icntrl_i) ) icntrl = icntrl_i ");
497	kppi.add_line(" if (present (rcntrl_i) ) rcntrl = rcntrl_i ");
498	kppi.add_line(" ");
499	kppi.add_line(" vl_glo = size(tempi,1) ");
500	kppi.add_line(" ");
501	kppi.add_line(" vl_dim_lo = VL_DIM ");
502	kppi.add_line(" DO k=1,VL_glo,vl_dim_lo ");
503	kppi.add_line(" is = k ");
504	kppi.add_line(" ie = min(k+vl_dim_lo-1,VL_glo) ");
505	kppi.add_line(" vl = ie-is+1 ");
506
507	kppi.add_line(" ");
508	if(kpp_switches.is_vector()) {
509	kppi.add_line(" C(1:vl,:) = Conc(is:ie,:) ");
510	} else {
511	kppi.add_line(" C(:) = Conc(is,:) ");
512	}
513
514	kppi.add_line(" ");
515	if(kpp_switches.is_vector()) {
516	kppi.add_line(" temp(1:vl) = tempi(is:ie) ");
517	} else {
518	kppi.add_line(" temp = tempi(is) ");
519	}
520	kppi.add_line(" ");
521	if(kpp_switches.is_vector()) {
522	kppi.add_line(" qvap(1:vl) = qvapi(is:ie) ");
523	} else {
524	kppi.add_line(" qvap = qvapi(is) ");
525	}
526	kppi.add_line(" ");
527	if(kpp_switches.is_vector()) {
528	kppi.add_line(" fakt(1:vl) = fakti(is:ie) ");
529	} else {
530	kppi.add_line(" fakt = fakti(is) ");
531	}
532
533	kppi.add_line(" ");
534	kppi.add_line(" CALL initialize ");
535
536	kppi.add_line(" ");
537	if(kpp_switches.is_vector()) {
538	kppi.add_line(" phot(1:vl,:) = photo(is:ie,:) ");
539	} else {
540	kppi.add_line(" phot(:) = photo(is,:) ");
541	}
542	kppi.add_line(" ");
543	kppi.add_line(" CALL update_rconst ");
544	kppi.add_line(" ");
545	kppi.add_line(" dt = time_step_len ");
546	kppi.add_line(" ");
547	kppi.add_line(" ! integrate from t=0 to t=dt ");
548	kppi.add_line(" CALL integrate(0._dp, dt, icntrl, rcntrl, istatus_u = istatus_u, ierr_u=ierr_u)");
549	kppi.add_line(" ");
550	kppi.add_line(" ");
551	if(kpp_switches.is_vector()) {
552	kppi.add_line(" Conc(is:ie,:) = C(1:VL,:) ");
553	} else {
554	kppi.add_line(" IF (PRESENT(l_debug) .AND. PRESENT(PE)) THEN ");
555	kppi.add_line(" IF (l_debug) CALL error_output(Conc(is,:),ierr_u, PE) ");
556	kppi.add_line(" ENDIF ");
557	kppi.add_line(" ");
558	kppi.add_line(" Conc(is,:) = C(:) ");
559	}
560
561	kppi.add_line(" ");
562	kppi.add_line(" ! Return Diagnostic Information ");
563	kppi.add_line(" ");
564	if(kpp_switches.is_vector()) {
565	kppi.add_line(" if(Present(ierrf)) ierrf(is:ie) = IERRV(1:VL) ");
566	kppi.add_line(" if(Present(xNacc)) xNacc(is:ie) = Kacc(1:VL) ");
567	kppi.add_line(" if(Present(xNrej)) xNrej(is:ie) = Krej(1:VL) ");
568	} else {
569	kppi.add_line(" if(Present(ierrf)) ierrf(is) = IERR_U ");
570	kppi.add_line(" if(Present(xNacc)) xNacc(is) = istatus_u(4) ");
571	kppi.add_line(" if(Present(xNrej)) xNrej(is) = istatus_u(5) ");
572	}
573	kppi.add_line(" ");
574	kppi.add_line(" if (present (istatus) ) then ");
575	if(kpp_switches.is_vector()) {
576	kppi.add_line(" istatus(4) = istatus(4) + sum(Kacc(1:VL)) ");
577	kppi.add_line(" istatus(5) = istatus(5) + sum(Krej(1:VL)) ");
578	kppi.add_line(" istatus(3) = istatus(4) + istatus(5) ");
579	kppi.add_line(" istatus(6) = istatus(6) + istatus_u(6) ");
580	kppi.add_line(" istatus(7) = istatus(7) + istatus_u(7) ");
581	} else {
582	kppi.add_line(" istatus(1:8) = istatus(1:8) + istatus_u(1:8) ");
583	}
584	kppi.add_line(" end if ");
585	kppi.add_line(" ");
586	kppi.add_line(" END DO ");
587	kppi.add_line(" ");
588
589	kppi.add_line(" ");
590	kppi.add_line("! Deallocate input arrays ");
591	kppi.add_line(" ");
592	for(iv=Vvar_list.begin();iv!=Vvar_list.end();iv++) {
593	// kppi.add_line(" if (allocated("+ iv->name +")) deallocate("+ iv->name +" ) ");
594	}
595
596	kppi.add_line(" ");
597	kppi.add_line(" data_loaded = .false. ");
598	kppi.add_line(" ");
599	kppi.add_line(" return ");
600	kppi.add_line("END SUBROUTINE chem_gasphase_integrate ");
601
602	// e5_subroutines.push_back(kppi);
603	kpp_subroutines.push_back(kppi);
604
605	return;
606	}
607

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |