Home

Context Navigation

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

Last change on this file since 3799 was 3799, checked in by forkel, 5 years ago
editing in kpp4palm: add statements for avoiding unused variables, remove $Id
File size: 25.5 KB

Rev	Line
[2696]	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	//
[3458]	12	//
[2696]	13	//Former revisions:
	14	//-----------------
[3799]	15	// Deleted $Id since document_changes does not work for C and C++ (15.03.2019, forkel)
	16	//
[3797]	17	// OpenMP version (15.03.2019, ketelsen)
	18	//
[3789]	19	// Added vector switch Kacc,Krej,IERRV, Commented add_line for istatf, (05.03.2019, forkel)
	20	// added ,pe after ierr_u,
	21	//
[3780]	22	// Added create_set_cs and cs_mech and get_mechanismname in module_header (05.03.2019, forkel)
	23	//
[3458]	24	// exclude kco_compress from handling by global_variables2vector (30.10.2018, forkel)
[2696]	25	//
[3458]	26	// Added automatic line with mechanism name (read mech_list) (25.09.2018, forkel)
[3298]	27	//
[3458]	28	// Added vl_glo = size(tempi,1) (20.09.2018, forkel)
[3298]	29	//
[3458]	30	// Removed creation of fill_ Subroutine and creation of calls thereof (18.09.2018, ketelsen)
[3298]	31	//
[3458]	32	// Fix in order not to loose the values of qvap and fakt (12.09.2018, forkel)
[3298]	33	//
[3458]	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)
[3298]	36	//
[3458]	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)
[2768]	39	//
[3298]	40	//
[3458]	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)
[2696]	49	//
[3458]	50	// Added phot (2017-09-14, forkel)
	51	// --> Change in module_header: Variables for photolyis added (2017-09-14, forkel)
[2696]	52	//
[3458]	53	// change of some output to lowercase with uppercase Fortran (2017, forkel)
[2696]	54	//
[3458]	55	// Intial version of KP4 adapted to PALM (Nov. 2016, ketelsen)
	56	//
[2696]	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;
[3298]	76	cout << "De_indexing " << kpp_switches.de_indexing() <<endl;
[2696]	77
	78	create_fortran_files_and_read();
	79
	80	// Generate first module lines
	81
[2768]	82	string first_line="MODULE " + module_name;
[3298]	83	mz_kpp.add_line(first_line);
	84	mz_kpp.add_line(" ");
[2696]	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
[3298]	135	cout << "##### Hier kpp_switches.is_vector " <<endl;
[2696]	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
[3458]	141	// Change global variables to vector (except for kp4_compress, which has already the right form)
[2696]	142
	143	for(it=kpp_subroutines.begin();it!=kpp_subroutines.end();it++) {
[3458]	144	if(it->get_name() != "kco_compress" ) {
	145	it->global_variables2vector (global_variable_list);
	146	}
[2696]	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	}
[3298]	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	}
[2696]	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++) {
[3298]	179	if(it->get_name() == "Update_RCONST") {
[2696]	180	it->edit_Update_RCONST(Vvar_list);
	181	}
[3298]	182
	183	if(it->get_name() == "Initialize") {
	184	it->edit_Initialize(Vvar_list);
	185	}
	186
[2696]	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
[3799]	199	// Create subroutine to communicate mechanism name to other modules
[3780]	200	create_set_cs();
	201
[2696]	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
[3298]	225	buf = prefix + "interface " + it->get_name() ;
[2696]	226	mz_kpp.add_line(buf);
[3298]	227	buf = prefix + " module procedure " + it->get_name();
[2696]	228	mz_kpp.add_line(buf);
[3298]	229	buf = prefix + "end interface " + it->get_name();
[2696]	230	mz_kpp.add_line(buf);
	231	mz_kpp.add_line(" ");
	232	}
	233
	234	mz_kpp.add_line(" ");
	235
[3797]	236	// Declare variables THREADPRIVATE for OpenMP version
[2696]	237
[3797]	238	mz_kpp.add_line(" ! OpenMP directives generated by kp4 ");
	239	mz_kpp.add_line(" ");
	240	mz_kpp.add_line(" !$OMP THREADPRIVATE (vl,vl_glo,is,ie,data_loaded)");
	241	mz_kpp.add_line(" !$OMP THREADPRIVATE (c,var,fix,rconst,time,temp,stepmin,cfactor)");
	242	mz_kpp.add_line(" !$OMP THREADPRIVATE (qvap,fakt,cs_mech,a,icntrl,rcntrl)");
	243	mz_kpp.add_line(" ");
	244	if(kpp_switches.is_vector()) {
	245	mz_kpp.add_line(" ! Vector mode Only ");
	246	mz_kpp.add_line(" ");
	247	mz_kpp.add_line(" !$OMP THREADPRIVATE (kacc,krej,ierrv)");
	248	mz_kpp.add_line(" !$OMP THREADPRIVATE (kpoints,kpoints_SAVE,index_org,done_check,index_step,cell_done)");
	249	mz_kpp.add_line(" !$OMP THREADPRIVATE (f_done,kacc_done,krej_done,ierr_done,compress_done)");
	250	mz_kpp.add_line(" ");
	251	}
	252
	253
[2696]	254	// Copy FORTRAN subroutines to mz_kpp
	255
	256	mz_kpp.add_line(" CONTAINS");
	257
	258	for(it=kpp_subroutines.begin();it!=kpp_subroutines.end();it++) {
	259	mz_kpp.add_line(" ");
	260	it->copy_to_MZ_KPP(mz_kpp);
	261	}
	262
	263	// Finish module
	264
[3298]	265	string last_line="end module " + module_name;
	266	mz_kpp.add_line("");
	267	mz_kpp.add_line(last_line);
[2696]	268
	269	// Write the complete module to file: mz_kpp.f
	270
	271	write_module_file();
	272
	273	return;
	274	}
	275
	276	void create_kpp_module::create_fortran_files_and_read() {
	277
	278	string name;
	279	ifstream in,in_c,in_b,in_i;
	280	fortran_file f_file;
	281	vector<fortran_file>::iterator it;
	282
	283	// Open file with list of FORTRAN routines
	284
	285	in.open("file_list");
	286	if( !in ) {
	287	cout << "cannot open " << endl; my_abort("file_list");
	288	}
	289
	290	// Create kpp_fortran routines
	291	while ( 1 ) {
	292	in >> name;
	293	if( in.eof() ) break;
	294	if( in.bad() ) my_abort("ERROR_READ_1");
	295	f_file.set_name(name);
	296	kpp_files.push_back(f_file);
	297	}
	298	in.close();
	299
	300	// Read FORTRAN code
	301
	302	for(it=kpp_files.begin();it!=kpp_files.end();it++) {
	303	it->read();
	304	}
	305
	306	// Open file with list of include files
	307
	308	in_c.open("include_list");
	309	if( !in_c ) {
	310	cout << "cannot open " << endl; my_abort("include_list");
	311	}
	312
	313	// Create kpp_includes vector
	314	while ( 1 ) {
	315	in_c >> name;
	316	if( in_c.eof() ) break;
	317	if( in_c.bad() ) my_abort("ERROR_READ_3");
	318	f_file.set_name(name);
	319	kpp_includes.push_back(f_file);
	320	}
	321	in_c.close();
	322
	323	// Read include files
	324
	325	for(it=kpp_includes.begin();it!=kpp_includes.end();it++) {
	326	it->read();
	327	}
	328
	329	// Read Ignore list
	330
	331	in_i.open("interface_ignore_list");
	332	if( !in_i ) {
	333	cout << "cannot open " << endl; my_abort("include_list");
	334	}
	335
	336	// Create kpp_includes vector
	337	while ( 1 ) {
	338	in_i >> name;
	339	if( in_i.eof() ) break;
	340	if( in_i.bad() ) my_abort("ERROR_READ_4");
	341	interface_ignore.push_back(name);
	342	}
	343	in_c.close();
	344
	345	}
	346
	347	void create_kpp_module::copy_files_to_subroutines () {
	348	string name;
	349	ifstream in;
	350	fortran_file s_file;
	351	vector<fortran_file>::iterator it;
	352
	353	// Open file with list of FORTRAN routines
	354
	355	in.open("subroutine_list");
	356	if( !in ) {
	357	cout << "cannot open " << endl; my_abort("subroutine_list");
	358	}
	359
	360	// Create vector kpp_subroutines
	361
	362	while ( 1 ) {
	363	in >> name;
	364	if( in.eof() ) break;
	365	if( in.bad() ) my_abort("ERROR_READ_S1");
	366	s_file.set_name(name);
	367	kpp_subroutines.push_back(s_file);
	368	}
	369	in.close();
	370
	371	header_variables.add_line(" ");
	372	header_variables.add_line("! variable definations from individual module headers ");
	373	header_variables.add_line(" ");
	374
	375	// Loop over all FORTRAN Files
	376
	377	for(it=kpp_files.begin();it!=kpp_files.end();it++) {
	378	it->copy_to_subroutine_vector(kpp_subroutines, header_variables);
	379	}
	380	}
	381
	382	void create_kpp_module::add_solver_to_subroutine_list () {
	383	fortran_file s_file;
	384
	385	string solver_name = getenv("KPP_SOLVER");
	386	cout << "KPP_SOLVER " <<solver_name <<endl;
	387
	388	s_file.set_name(solver_name);
	389	s_file.read();
	390	kpp_subroutines.push_back(s_file);
	391
	392	return;
	393	}
	394
	395	void create_kpp_module::generate_module_header() {
	396
	397	string buf;
	398	ifstream in;
	399	ifstream in_e5;
	400	program_line line;
	401	vector<fortran_file>::iterator it;
	402	char distr[2];
	403	string diline;
	404
[3298]	405	// Read mechanism from mech_list
	406
	407	in.open("mech_list");
	408	if( !in ) {
	409	cout << "cannot open " << endl; my_abort("mech_list");
	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	in.close();
	420
	421
[2696]	422	// Read Modul Header from file $MZ_KPP_HOME/templates/module_header
	423
	424	in.open("module_header");
	425	if( !in ) {
	426	cout << "cannot open " << endl; my_abort("module_header");
	427	}
	428
	429	while ( 1 ) {
	430	getline (in, buf);
	431	if( in.eof() ) break;
	432	if( in.bad() ) my_abort("ERROR_READ_4");
	433	line.set_line(buf);
	434	mz_kpp.add_line(line);
	435	}
	436	mz_kpp.add_line(" ");
	437	mz_kpp.add_line("! Variables used for vector mode ");
	438	mz_kpp.add_line(" ");
	439	if(kpp_switches.is_vector()) {
[3298]	440	mz_kpp.add_line(" logical,parameter :: L_VECTOR = .TRUE. ");
[2696]	441	} else {
[3298]	442	mz_kpp.add_line(" logical,parameter :: L_VECTOR = .FALSE. ");
[2696]	443	}
	444	// mz_pj_20070531+
	445	sprintf(distr,"%i",kpp_switches.de_indexing());
	446	diline = distr ;
[3298]	447	mz_kpp.add_line(" integer,parameter :: I_LU_DI = " + diline );
[2696]	448	// mz_pj_20070531-
	449
[3298]	450	mz_kpp.add_line(" integer,parameter :: VL_DIM = "
	451	+ kpp_switches.get_vector_length() );
	452	mz_kpp.add_line(" integer :: vl ");
[2696]	453	mz_kpp.add_line(" ");
[3298]	454	mz_kpp.add_line(" integer :: VL_glo ");
	455	mz_kpp.add_line(" integer :: is,ie ");
[2696]	456	mz_kpp.add_line(" ");
[3298]	457	mz_kpp.add_line(" ");
[3789]	458	if(kpp_switches.is_vector()) {
	459	mz_kpp.add_line(" integer, dimension(VL_dim) :: Kacc,Krej ");
	460	mz_kpp.add_line(" integer, dimension(VL_dim) :: IERRV ");
	461	}
[3298]	462	mz_kpp.add_line(" logical :: data_loaded = .false. ");
[3797]	463	if(kpp_switches.is_vector()) {
	464	mz_kpp.add_line(" REAL(dp),POINTER,DIMENSION(:,:),CONTIGUOUS :: var ");
	465	} else {
	466	mz_kpp.add_line(" REAL(dp),POINTER,DIMENSION(:),CONTIGUOUS :: var ");
	467	}
[2696]	468	in.close();
	469
	470	return;
	471	}
	472
	473	void create_kpp_module::write_module_file() {
	474	ofstream out;
	475	ofstream out_e5;
	476
	477	string out_file = "kk_kpp.f90";
	478	out.open(out_file.c_str(), ios::out);
	479	if( !out ) {
	480	cout << "cannot open " << endl; my_abort(out_file);
	481	}
	482
	483	mz_kpp.write_file (out);
	484
	485	out.close();
	486
	487
	488	return;
	489	}
	490
[3780]	491	void create_kpp_module::create_set_cs() {
	492	fortran_file kppi;
	493	vector<Vvar>::iterator iv;
	494	string xline;
	495
	496	string buf;
	497	ifstream in;
	498	program_line line;
	499
	500	kppi.set_name("get_mechanismname");
	501	kppi.add_line("SUBROUTINE get_mechanismname ");
	502	kppi.add_line(" ");
	503	kppi.add_line(" IMPLICIT NONE ");
	504	// Read mechanism from set_cm
	505	// Tis got an own own subroutine to aviod being called at each timestep
	506
	507	in.open("set_cm");
	508	if( !in ) {
	509	cout << "cannot open " << endl; my_abort("set_cm");
	510	}
	511
	512	while ( 1 ) {
	513	getline (in, buf);
	514	if( in.eof() ) break;
	515	if( in.bad() ) my_abort("ERROR_READ_4");
	516	line.set_line(buf);
	517	kppi.add_line(line);
	518	}
	519	in.close();
	520
	521	kppi.add_line(" ");
	522	kppi.add_line(" return ");
	523	kppi.add_line("END SUBROUTINE get_mechanismname ");
	524	kppi.add_line(" ");
	525	kpp_subroutines.push_back(kppi);
	526
	527	return;
	528	}
	529
	530
[2696]	531	void create_kpp_module::create_kpp_integrate() {
	532	fortran_file kppi;
	533	vector<Vvar>::iterator iv;
[3780]	534	string xline;
[2696]	535
[3780]	536
[2696]	537	kppi.set_name("chem_gasphase_integrate");
	538
[3298]	539	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 ) ");
[2696]	540	kppi.add_line(" ");
	541	kppi.add_line(" IMPLICIT NONE ");
	542	kppi.add_line(" ");
	543
[3298]	544	kppi.add_line(" REAL(dp), INTENT(IN) :: time_step_len ");
[2696]	545	kppi.add_line(" REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: conc ");
[3298]	546	kppi.add_line(" REAL(dp), DIMENSION(:,:), INTENT(IN) :: photo ");
	547	kppi.add_line(" REAL(dp), DIMENSION(:), INTENT(IN) :: tempi ");
	548	kppi.add_line(" REAL(dp), DIMENSION(:), INTENT(IN) :: qvapi ");
	549	kppi.add_line(" REAL(dp), DIMENSION(:), INTENT(IN) :: fakti ");
	550	kppi.add_line(" INTEGER, INTENT(OUT), OPTIONAL :: ierrf(:) ");
	551	kppi.add_line(" INTEGER, INTENT(OUT), OPTIONAL :: xNacc(:) ");
	552	kppi.add_line(" INTEGER, INTENT(OUT), OPTIONAL :: xNrej(:) ");
	553	kppi.add_line(" INTEGER, INTENT(INOUT), OPTIONAL :: istatus(:) ");
	554	kppi.add_line(" INTEGER, INTENT(IN), OPTIONAL :: PE ");
	555	kppi.add_line(" LOGICAL, INTENT(IN), OPTIONAL :: l_debug ");
	556	kppi.add_line(" INTEGER, DIMENSION(nkppctrl),INTENT(IN), OPTIONAL :: icntrl_i ");
	557	kppi.add_line(" REAL(dp), DIMENSION(nkppctrl),INTENT(IN), OPTIONAL :: rcntrl_i ");
[2696]	558	kppi.add_line(" ");
	559	kppi.add_line(" INTEGER :: k ! loop variable ");
	560	kppi.add_line(" REAL(dp) :: dt ");
[3298]	561	kppi.add_line(" integer, dimension(20) :: istatus_u ");
	562	kppi.add_line(" integer :: ierr_u ");
[3789]	563	// kppi.add_line(" integer :: istatf ");
[3298]	564	kppi.add_line(" integer :: vl_dim_lo ");
[2696]	565	kppi.add_line(" ");
	566	kppi.add_line(" ");
[3298]	567	kppi.add_line(" if (present (istatus) ) istatus = 0 ");
	568	kppi.add_line(" if (present (icntrl_i) ) icntrl = icntrl_i ");
	569	kppi.add_line(" if (present (rcntrl_i) ) rcntrl = rcntrl_i ");
[2696]	570	kppi.add_line(" ");
[3797]	571	if(kpp_switches.is_vector()) {
	572	kppi.add_line(" var => c(:,1:nvar) ");
	573	} else {
	574	kppi.add_line(" var => c(1:nvar) ");
	575	}
	576	kppi.add_line(" ");
[3298]	577	kppi.add_line(" vl_glo = size(tempi,1) ");
	578	kppi.add_line(" ");
	579	kppi.add_line(" vl_dim_lo = VL_DIM ");
	580	kppi.add_line(" DO k=1,VL_glo,vl_dim_lo ");
[2696]	581	kppi.add_line(" is = k ");
[3298]	582	kppi.add_line(" ie = min(k+vl_dim_lo-1,VL_glo) ");
[2696]	583	kppi.add_line(" vl = ie-is+1 ");
	584
	585	kppi.add_line(" ");
	586	if(kpp_switches.is_vector()) {
[3298]	587	kppi.add_line(" C(1:vl,:) = Conc(is:ie,:) ");
[2696]	588	} else {
[3298]	589	kppi.add_line(" C(:) = Conc(is,:) ");
[2696]	590	}
[3298]	591
[2696]	592	kppi.add_line(" ");
	593	if(kpp_switches.is_vector()) {
[3298]	594	kppi.add_line(" temp(1:vl) = tempi(is:ie) ");
[2696]	595	} else {
[3298]	596	kppi.add_line(" temp = tempi(is) ");
[2696]	597	}
	598	kppi.add_line(" ");
	599	if(kpp_switches.is_vector()) {
[3298]	600	kppi.add_line(" qvap(1:vl) = qvapi(is:ie) ");
	601	} else {
	602	kppi.add_line(" qvap = qvapi(is) ");
	603	}
	604	kppi.add_line(" ");
	605	if(kpp_switches.is_vector()) {
	606	kppi.add_line(" fakt(1:vl) = fakti(is:ie) ");
	607	} else {
	608	kppi.add_line(" fakt = fakti(is) ");
	609	}
	610
	611	kppi.add_line(" ");
	612	kppi.add_line(" CALL initialize ");
	613
	614	kppi.add_line(" ");
	615	if(kpp_switches.is_vector()) {
[2696]	616	kppi.add_line(" phot(1:vl,:) = photo(is:ie,:) ");
	617	} else {
	618	kppi.add_line(" phot(:) = photo(is,:) ");
	619	}
	620	kppi.add_line(" ");
	621	kppi.add_line(" CALL update_rconst ");
	622	kppi.add_line(" ");
	623	kppi.add_line(" dt = time_step_len ");
	624	kppi.add_line(" ");
	625	kppi.add_line(" ! integrate from t=0 to t=dt ");
	626	kppi.add_line(" CALL integrate(0._dp, dt, icntrl, rcntrl, istatus_u = istatus_u, ierr_u=ierr_u)");
	627	kppi.add_line(" ");
	628	kppi.add_line(" ");
	629	if(kpp_switches.is_vector()) {
[3298]	630	kppi.add_line(" Conc(is:ie,:) = C(1:VL,:) ");
[2696]	631	} else {
[3298]	632	kppi.add_line(" IF (PRESENT(l_debug) .AND. PRESENT(PE)) THEN ");
[3789]	633	kppi.add_line(" IF (l_debug) CALL error_output(Conc(is,:),ierr_u,pe) ");
[3298]	634	kppi.add_line(" ENDIF ");
	635	kppi.add_line(" ");
	636	kppi.add_line(" Conc(is,:) = C(:) ");
[2696]	637	}
	638
	639	kppi.add_line(" ");
	640	kppi.add_line(" ! Return Diagnostic Information ");
	641	kppi.add_line(" ");
	642	if(kpp_switches.is_vector()) {
[3298]	643	kppi.add_line(" if(Present(ierrf)) ierrf(is:ie) = IERRV(1:VL) ");
	644	kppi.add_line(" if(Present(xNacc)) xNacc(is:ie) = Kacc(1:VL) ");
	645	kppi.add_line(" if(Present(xNrej)) xNrej(is:ie) = Krej(1:VL) ");
[2696]	646	} else {
[3298]	647	kppi.add_line(" if(Present(ierrf)) ierrf(is) = IERR_U ");
	648	kppi.add_line(" if(Present(xNacc)) xNacc(is) = istatus_u(4) ");
	649	kppi.add_line(" if(Present(xNrej)) xNrej(is) = istatus_u(5) ");
[2696]	650	}
	651	kppi.add_line(" ");
[3298]	652	kppi.add_line(" if (present (istatus) ) then ");
[2696]	653	if(kpp_switches.is_vector()) {
[3298]	654	kppi.add_line(" istatus(4) = istatus(4) + sum(Kacc(1:VL)) ");
	655	kppi.add_line(" istatus(5) = istatus(5) + sum(Krej(1:VL)) ");
[2696]	656	kppi.add_line(" istatus(3) = istatus(4) + istatus(5) ");
	657	kppi.add_line(" istatus(6) = istatus(6) + istatus_u(6) ");
	658	kppi.add_line(" istatus(7) = istatus(7) + istatus_u(7) ");
	659	} else {
	660	kppi.add_line(" istatus(1:8) = istatus(1:8) + istatus_u(1:8) ");
	661	}
[3298]	662	kppi.add_line(" end if ");
[2696]	663	kppi.add_line(" ");
	664	kppi.add_line(" END DO ");
	665	kppi.add_line(" ");
	666
	667	kppi.add_line(" ");
	668	kppi.add_line("! Deallocate input arrays ");
	669	kppi.add_line(" ");
	670	for(iv=Vvar_list.begin();iv!=Vvar_list.end();iv++) {
[3298]	671	// kppi.add_line(" if (allocated("+ iv->name +")) deallocate("+ iv->name +" ) ");
[2696]	672	}
	673
	674	kppi.add_line(" ");
	675	kppi.add_line(" data_loaded = .false. ");
	676	kppi.add_line(" ");
[3298]	677	kppi.add_line(" return ");
[3780]	678	kppi.add_line("END SUBROUTINE chem_gasphase_integrate ");
[2696]	679
	680	// e5_subroutines.push_back(kppi);
	681	kpp_subroutines.push_back(kppi);
	682
	683	return;
	684	}
	685

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |