Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |