Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |