Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |