Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |