Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |