Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |