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

Last change on this file since 3820 was 3820, checked in by forkel, 3 years ago

renaming of get_mechanismname, do_emiss and do_depo, sorting in chem_modules

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