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, 2 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
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//
13//Former revisions:
14//-----------------
15// alphabetic ordering of PUBLIC statments in templates/module_header     (28.03.2019, forkel)
16//
17// renamed get_mechanismname to get_mechanism_name                        (26.03.2019, forkel)
18//
19//
20// Deleted $Id since document_changes does not work for C and C++         (15.03.2019, forkel)
21//
22// OpenMP version    (15.03.2019, ketelsen)
23//
24// Added vector switch Kacc,Krej,IERRV, Commented add_line for istatf,    (05.03.2019, forkel)
25//      added ,pe after ierr_u,         
26//
27// Added create_set_cs and cs_mech and get_mechanismname in module_header (05.03.2019, forkel)
28//
29// exclude kco_compress from handling by global_variables2vector (30.10.2018, forkel)
30//
31// Added  automatic line with mechanism name (read mech_list) (25.09.2018, forkel)
32//
33// Added  vl_glo = size(tempi,1) (20.09.2018, forkel)
34//
35// Removed creation of fill_ Subroutine and creation of calls thereof (18.09.2018, ketelsen)
36//
37// Fix in order not to loose the values of qvap and fakt (12.09.2018, forkel)
38//
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)
41//
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)
44//
45//
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)
54//
55// Added phot                                                 (2017-09-14, forkel)
56// --> Change in module_header: Variables for photolyis added (2017-09-14, forkel)
57//
58// change of some output to lowercase with uppercase Fortran (2017, forkel)
59//
60// Intial version of KP4 adapted to PALM                      (Nov. 2016, ketelsen)
61//
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
69void 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;
81   cout << "De_indexing " << kpp_switches.de_indexing() <<endl;
82
83   create_fortran_files_and_read();
84
85// Generate first module lines
86
87     string first_line="MODULE " + module_name;
88   mz_kpp.add_line(first_line);
89   mz_kpp.add_line(" ");
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
140   cout << "##### Hier kpp_switches.is_vector          " <<endl;
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
146//   Change global variables to vector (except for kp4_compress, which has already the right form)
147     
148     for(it=kpp_subroutines.begin();it!=kpp_subroutines.end();it++) {
149       if(it->get_name() != "kco_compress" ) {
150         it->global_variables2vector (global_variable_list);
151       }
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       }
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       }
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++) {
184     if(it->get_name() == "Update_RCONST") {
185       it->edit_Update_RCONST(Vvar_list);
186     }
187
188     if(it->get_name() == "Initialize") {
189       it->edit_Initialize(Vvar_list);
190     }
191
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
204// Create subroutine to communicate mechanism name to other modules
205   create_set_cs();
206
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
230     buf = prefix + "interface            " + it->get_name() ;
231     mz_kpp.add_line(buf);
232     buf = prefix + "  module procedure   " + it->get_name();
233     mz_kpp.add_line(buf);
234     buf = prefix + "end interface        " + it->get_name();
235     mz_kpp.add_line(buf);
236     mz_kpp.add_line(" ");
237   }
238
239   mz_kpp.add_line(" ");
240
241// Declare variables THREADPRIVATE for OpenMP version
242
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
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
270   string last_line="end module " + module_name;
271   mz_kpp.add_line("");
272   mz_kpp.add_line(last_line);
273
274// Write the complete module to file: mz_kpp.f
275
276   write_module_file();
277
278   return;
279}
280
281void 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
352void 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
387void 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
400void 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
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
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()) {
445       mz_kpp.add_line("  logical,parameter          :: L_VECTOR = .TRUE.             ");
446   } else {
447       mz_kpp.add_line("  logical,parameter          :: L_VECTOR = .FALSE.            ");
448   }
449//  mz_pj_20070531+
450   sprintf(distr,"%i",kpp_switches.de_indexing());
451   diline = distr ;
452   mz_kpp.add_line("  integer,parameter          :: I_LU_DI = " + diline );
453//  mz_pj_20070531-
454
455   mz_kpp.add_line("  integer,parameter          :: VL_DIM = " 
456                 + kpp_switches.get_vector_length() ); 
457   mz_kpp.add_line("  integer                     :: vl                              "); 
458   mz_kpp.add_line("                                                                 "); 
459   mz_kpp.add_line("  integer                     :: VL_glo                          "); 
460   mz_kpp.add_line("  integer                     :: is,ie                           "); 
461   mz_kpp.add_line("                                                                 "); 
462   mz_kpp.add_line("                                                                 "); 
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   }
467   mz_kpp.add_line("  logical                     :: data_loaded = .false.             "); 
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   }
473   in.close();
474
475   return;
476}
477
478void 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
496void 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
505   kppi.set_name("get_mechanism_name");
506   kppi.add_line("SUBROUTINE get_mechanism_name                                       ");
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                                                            ");
528   kppi.add_line("END SUBROUTINE get_mechanism_name                                   ");
529   kppi.add_line("                                                                    ");
530   kpp_subroutines.push_back(kppi);
531
532   return;
533}
534
535
536void create_kpp_module::create_kpp_integrate() {
537   fortran_file          kppi;
538   vector<Vvar>::iterator               iv;
539   string                               xline;
540
541
542   kppi.set_name("chem_gasphase_integrate");
543
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 )  ");
545   kppi.add_line("                                                                    ");
546   kppi.add_line("  IMPLICIT NONE                                                     ");
547   kppi.add_line("                                                                    ");
548
549   kppi.add_line("  REAL(dp), INTENT(IN)                   :: time_step_len           ");
550   kppi.add_line("  REAL(dp),  DIMENSION(:,:),  INTENT(INOUT) :: conc                    ");
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         ");
563   kppi.add_line("                                                                    ");
564   kppi.add_line("  INTEGER                                 :: k   ! loop variable     ");
565   kppi.add_line("  REAL(dp)                                :: dt                      ");
566   kppi.add_line("  integer, dimension(20)                 :: istatus_u               ");
567   kppi.add_line("  integer                                :: ierr_u                  ");
568// kppi.add_line("  integer                                :: istatf                  ");
569   kppi.add_line("  integer                                :: vl_dim_lo               ");
570   kppi.add_line("                                                                    ");
571   kppi.add_line("                                                                    ");
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                      ");
575   kppi.add_line("                                                                    ");
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("                                                                    ");
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                                           ");
586   kppi.add_line("    is = k                                                          ");
587   kppi.add_line("    ie = min(k+vl_dim_lo-1,VL_glo)                                  ");
588   kppi.add_line("    vl = ie-is+1                                                    ");
589
590   kppi.add_line("                                                                    ");
591   if(kpp_switches.is_vector()) {
592     kppi.add_line("    C(1:vl,:) = Conc(is:ie,:)                                     ");
593   } else {
594     kppi.add_line("    C(:) = Conc(is,:)                                             ");
595   }
596
597   kppi.add_line("                                                                    ");
598   if(kpp_switches.is_vector()) {
599     kppi.add_line("    temp(1:vl) = tempi(is:ie)                                     ");
600   } else {
601     kppi.add_line("    temp = tempi(is)                                              ");
602   }
603   kppi.add_line("                                                                    ");
604   if(kpp_switches.is_vector()) {
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()) {
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()) {
635     kppi.add_line("    Conc(is:ie,:) = C(1:VL,:)                                     ");
636   } else {
637     kppi.add_line("   IF (PRESENT(l_debug) .AND. PRESENT(PE)) THEN                       ");
638     kppi.add_line("      IF (l_debug) CALL error_output(Conc(is,:),ierr_u,pe)            ");
639     kppi.add_line("   ENDIF                                                              ");
640     kppi.add_line("                                                                      ");
641     kppi.add_line("    Conc(is,:) = C(:)                                                 ");
642   }
643
644   kppi.add_line("                                                                    ");
645   kppi.add_line("    ! Return Diagnostic Information                                 ");
646   kppi.add_line("                                                                    ");
647   if(kpp_switches.is_vector()) {
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)               ");
651   } else {
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)                ");
655   }
656   kppi.add_line("                                                                    ");
657   kppi.add_line("    if (present (istatus) )  then                                   ");
658   if(kpp_switches.is_vector()) {
659     kppi.add_line("      istatus(4) =   istatus(4) + sum(Kacc(1:VL))                  ");
660     kppi.add_line("      istatus(5) =   istatus(5) + sum(Krej(1:VL))                  ");
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   }
667   kppi.add_line("    end if                                                          ");
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++) {
676//     kppi.add_line("  if (allocated("+ iv->name +"))   deallocate("+ iv->name +" )    ");
677   }
678
679   kppi.add_line("                                                                    ");
680   kppi.add_line("  data_loaded = .false.                                             ");
681   kppi.add_line("                                                                    ");
682   kppi.add_line("  return                                                            ");
683   kppi.add_line("END SUBROUTINE chem_gasphase_integrate                              ");
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.