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

Last change on this file since 3458 was 3458, checked in by kanani, 3 years ago

Reintegrated fixes/changes from branch chemistry

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