- Timestamp:
- Oct 2, 2018 12:21:11 PM (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/create_kpp_module.C
r2768 r3298 17 17 //----------------- 18 18 //$Id: create_kpp_module.C 2470 2017-09-14 13:56:42Z forkel $ 19 // forkel 25.09.2018: added automatic line with meechanism name (read mech_list) 20 // forkel 20.09.2018: added vl_glo = size(tempi,1) 21 // ketelsen 18.09.2018: Removed creation of fill_ Subroutine and creation of calls thereof 22 // 23 // forkel 12.09.2018: Fix in order not to loose the values of qvap and fakt 24 // 25 // forkel 03.09.2018: Bug fix: moved kppi.add_line(" CALL initialize after fakt = fakti(is) 26 // forkel June 2018: added qvap and fakt, re-established original uppercase/lowercase 27 // Deleted definition of qvap,fakt in create_kpp_integrate again 28 // 29 // ketelsen July 2018: Changes for vector mode (edit_WAXPY, edit_FunTemplate, edit_JacTemplate, 30 // some cleanup of comments, various changes in create_kpp_integrate) 31 // 32 // forkel June 2018: Added qvap and fakt, re-established original uppercase/lowercase 33 // 34 // 35 / 2017-09-14 13:56:42Z forkel $ 19 36 // Removed preprocessor directive __chem again 20 37 // 21 / 2017-09-14 13:56:42Z forkel $38 // 2017-09-14 13:56:42Z forkel $ 22 39 // 23 40 // … … 25 42 //change of some output to lowercase with uppercase Fortran 26 43 // 27 // Nov 2016: Intial version (Klaus Ketelsen)44 // Nov 2016: Intial version of KP4 adapted to PALM (Klaus Ketelsen) 28 45 // 29 46 … … 49 66 cout << "Create " << module_name << " from kpp Fortran sources" <<endl; 50 67 cout << "Vector mode " << kpp_switches.is_vector() <<endl; 68 cout << "De_indexing " << kpp_switches.de_indexing() <<endl; 51 69 52 70 create_fortran_files_and_read(); 53 cout << "## after create_fortran_files_and_read " <<endl;54 71 55 72 // Generate first module lines 56 73 57 74 string first_line="MODULE " + module_name; 58 mz_kpp.add_line(first_line); 59 mz_kpp.add_line("!"); 60 // mz_kpp.add_line("#if defined( __chem )"); 61 // mz_kpp.add_line(" "); 75 mz_kpp.add_line(first_line); 76 mz_kpp.add_line(" "); 62 77 63 78 // string e5_line = first_line +"_e5"; … … 93 108 it->edit_fortran (); 94 109 } 95 cout << "## after edit FORTRAN files " <<endl;96 110 97 111 // Generate a list of single subroutines from kpp-files … … 109 123 } 110 124 111 112 if(kpp_switches.is_vector()) { 113 125 if(kpp_switches.is_vector()) { 126 127 cout << "##### Hier kpp_switches.is_vector " <<endl; 114 128 // Change header section 115 129 for(it=kpp_includes.begin();it!=kpp_includes.end();it++) { … … 138 152 it->edit_Fun(); 139 153 } 140 } 141 cout << "## after Edit individual subroutines " <<endl; 154 if(it->get_name() == "WAXPY" ) { 155 it->edit_WAXPY(); 156 } 157 if(it->get_name() == "FunTemplate" ) { 158 it->edit_FunTemplate(); 159 } 160 if(it->get_name() == "JacTemplate" ) { 161 it->edit_JacTemplate(); 162 } 163 } 164 // cout << "## after Edit individual subroutines " <<endl; 142 165 143 166 } … … 146 169 147 170 for(it=kpp_subroutines.begin();it!=kpp_subroutines.end();it++) { 148 if(it->get_name() == " update_rconst") {171 if(it->get_name() == "Update_RCONST") { 149 172 it->edit_Update_RCONST(Vvar_list); 150 173 } 174 175 if(it->get_name() == "Initialize") { 176 it->edit_Initialize(Vvar_list); 177 } 178 151 179 } 152 180 … … 184 212 } 185 213 186 buf = prefix + " INTERFACE" + it->get_name() ;214 buf = prefix + "interface " + it->get_name() ; 187 215 mz_kpp.add_line(buf); 188 buf = prefix + " MODULE PROCEDURE" + it->get_name();216 buf = prefix + " module procedure " + it->get_name(); 189 217 mz_kpp.add_line(buf); 190 buf = prefix + " END INTERFACE" + it->get_name();218 buf = prefix + "end interface " + it->get_name(); 191 219 mz_kpp.add_line(buf); 192 220 mz_kpp.add_line(" "); 193 221 } 194 for(iv=Vvar_list.begin();iv!=Vvar_list.end();iv++) {195 string buf;196 197 string sub_name = "fill_" + iv->name;198 buf = " INTERFACE " + sub_name;199 mz_kpp.add_line(buf);200 buf = " MODULE PROCEDURE " + sub_name;201 mz_kpp.add_line(buf);202 buf = " END INTERFACE " + sub_name;203 mz_kpp.add_line(buf);204 buf = " PUBLIC " + sub_name;205 mz_kpp.add_line(buf);206 mz_kpp.add_line(" ");207 }208 222 209 223 mz_kpp.add_line(" "); 210 224 211 for(iv=Vvar_list.begin();iv!=Vvar_list.end();iv++) {212 create_fill_routine(kpp_subroutines, *iv);213 }214 225 215 226 // Copy FORTRAN subroutines to mz_kpp … … 224 235 // Finish module 225 236 226 // mz_kpp.add_line("#endif"); 227 string last_line="END MODULE " + module_name; 228 mz_kpp.add_line(""); 229 mz_kpp.add_line(last_line); 230 231 // Handle e5 module 232 233 // for(it=e5_subroutines.begin();it!=e5_subroutines.end();it++) { 234 // e5_kpp.add_line(" "); 235 // it->copy_to_MZ_KPP(e5_kpp); 236 // } 237 238 // last_line = last_line + "_e5"; 239 // e5_kpp.add_line(" "); 240 // e5_kpp.add_line(last_line); 237 string last_line="end module " + module_name; 238 mz_kpp.add_line(""); 239 mz_kpp.add_line(last_line); 241 240 242 241 // Write the complete module to file: mz_kpp.f … … 376 375 string diline; 377 376 377 // Read mechanism from mech_list 378 379 in.open("mech_list"); 380 if( !in ) { 381 cout << "cannot open " << endl; my_abort("mech_list"); 382 } 383 384 while ( 1 ) { 385 getline (in, buf); 386 if( in.eof() ) break; 387 if( in.bad() ) my_abort("ERROR_READ_4"); 388 line.set_line(buf); 389 mz_kpp.add_line(line); 390 } 391 in.close(); 392 393 378 394 // Read Modul Header from file $MZ_KPP_HOME/templates/module_header 379 395 … … 394 410 mz_kpp.add_line(" "); 395 411 if(kpp_switches.is_vector()) { 396 mz_kpp.add_line(" LOGICAL, PARAMETER :: l_vector= .TRUE. ");397 } else { 398 mz_kpp.add_line(" LOGICAL, PARAMETER :: l_vector= .FALSE. ");412 mz_kpp.add_line(" logical,parameter :: L_VECTOR = .TRUE. "); 413 } else { 414 mz_kpp.add_line(" logical,parameter :: L_VECTOR = .FALSE. "); 399 415 } 400 416 // mz_pj_20070531+ 401 417 sprintf(distr,"%i",kpp_switches.de_indexing()); 402 418 diline = distr ; 403 mz_kpp.add_line(" INTEGER, PARAMETER :: i_lu_di= " + diline );419 mz_kpp.add_line(" integer,parameter :: I_LU_DI = " + diline ); 404 420 // mz_pj_20070531- 405 421 406 mz_kpp.add_line(" INTEGER, PARAMETER :: vl_dim= "407 408 mz_kpp.add_line(" INTEGER:: vl ");422 mz_kpp.add_line(" integer,parameter :: VL_DIM = " 423 + kpp_switches.get_vector_length() ); 424 mz_kpp.add_line(" integer :: vl "); 409 425 mz_kpp.add_line(" "); 410 mz_kpp.add_line(" INTEGER :: vl_glo ");411 mz_kpp.add_line(" INTEGER:: is,ie ");426 mz_kpp.add_line(" integer :: VL_glo "); 427 mz_kpp.add_line(" integer :: is,ie "); 412 428 mz_kpp.add_line(" "); 413 mz_kpp.add_line(" INTEGER, DIMENSION(vl_dim) :: kacc, krej "); 414 mz_kpp.add_line(" INTEGER, DIMENSION(vl_dim) :: ierrv "); 415 mz_kpp.add_line(" LOGICAL :: data_loaded = .false. "); 416 429 mz_kpp.add_line(" "); 430 mz_kpp.add_line(" integer, dimension(VL_dim) :: Kacc,Krej "); 431 mz_kpp.add_line(" integer, dimension(VL_dim) :: IERRV "); 432 mz_kpp.add_line(" logical :: data_loaded = .false. "); 433 417 434 in.close(); 418 419 // in_e5.open("module_header_e5");420 // if( !in_e5 ) {421 // cout << "cannot open " << endl; my_abort("module_header_e5");422 // }423 424 // while ( 1 ) {425 // getline (in_e5, buf);426 // if( in_e5.eof() ) break;427 // if( in_e5.bad() ) my_abort("ERROR_READ_4");428 // line.set_line(buf);429 // e5_kpp.add_line(line);430 // }431 // in_e5.close();432 435 433 436 return; … … 448 451 out.close(); 449 452 450 // out_file = "kk_mecca_kpp_e5.f90";451 // out_e5.open(out_file.c_str(), ios::out);452 // if( !out_e5 ) {453 // cout << "cannot open " << endl; my_abort(out_file);454 // }455 456 // e5_kpp.write_file (out_e5);457 458 // out_e5.close();459 453 460 454 return; … … 468 462 kppi.set_name("chem_gasphase_integrate"); 469 463 470 kppi.add_line("SUBROUTINE chem_gasphase_integrate (time_step_len, conc, temp k, photo, ierrf, xnacc, xnrej, istatus, l_debug, pe) ");464 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 ) "); 471 465 kppi.add_line(" "); 472 466 kppi.add_line(" IMPLICIT NONE "); 473 467 kppi.add_line(" "); 474 468 475 kppi.add_line(" REAL(dp), INTENT(IN):: time_step_len ");469 kppi.add_line(" REAL(dp), INTENT(IN) :: time_step_len "); 476 470 kppi.add_line(" REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: conc "); 477 kppi.add_line(" REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: photo "); 478 kppi.add_line(" REAL(dp), DIMENSION(:), INTENT(IN) :: tempk "); 479 kppi.add_line(" INTEGER, INTENT(OUT), OPTIONAL :: ierrf(:) "); 480 kppi.add_line(" INTEGER, INTENT(OUT), OPTIONAL :: xnacc(:) "); 481 kppi.add_line(" INTEGER, INTENT(OUT), OPTIONAL :: xnrej(:) "); 482 kppi.add_line(" INTEGER, INTENT(INOUT), OPTIONAL :: istatus(:) "); 483 kppi.add_line(" INTEGER, INTENT(IN), OPTIONAL :: pe "); 484 kppi.add_line(" LOGICAL, INTENT(IN), OPTIONAL :: l_debug "); 471 kppi.add_line(" REAL(dp), DIMENSION(:,:), INTENT(IN) :: photo "); 472 kppi.add_line(" REAL(dp), DIMENSION(:), INTENT(IN) :: tempi "); 473 kppi.add_line(" REAL(dp), DIMENSION(:), INTENT(IN) :: qvapi "); 474 kppi.add_line(" REAL(dp), DIMENSION(:), INTENT(IN) :: fakti "); 475 kppi.add_line(" INTEGER, INTENT(OUT), OPTIONAL :: ierrf(:) "); 476 kppi.add_line(" INTEGER, INTENT(OUT), OPTIONAL :: xNacc(:) "); 477 kppi.add_line(" INTEGER, INTENT(OUT), OPTIONAL :: xNrej(:) "); 478 kppi.add_line(" INTEGER, INTENT(INOUT), OPTIONAL :: istatus(:) "); 479 kppi.add_line(" INTEGER, INTENT(IN), OPTIONAL :: PE "); 480 kppi.add_line(" LOGICAL, INTENT(IN), OPTIONAL :: l_debug "); 481 kppi.add_line(" INTEGER, DIMENSION(nkppctrl),INTENT(IN), OPTIONAL :: icntrl_i "); 482 kppi.add_line(" REAL(dp), DIMENSION(nkppctrl),INTENT(IN), OPTIONAL :: rcntrl_i "); 485 483 kppi.add_line(" "); 486 484 kppi.add_line(" INTEGER :: k ! loop variable "); 487 485 kppi.add_line(" REAL(dp) :: dt "); 488 kppi.add_line(" INTEGER, DIMENSION(20) :: istatus_u "); 489 kppi.add_line(" INTEGER :: ierr_u "); 490 kppi.add_line(" "); 491 kppi.add_line(" "); 492 kppi.add_line(" if (present (istatus) ) istatus = 0 "); 493 kppi.add_line(" "); 494 // kppi.add_line(" vk_glo = size(tempk,1) "); 495 // kppi.add_line(" "); 496 497 kppi.add_line(" DO k=1,vl_glo,vl_dim "); 486 kppi.add_line(" integer, dimension(20) :: istatus_u "); 487 kppi.add_line(" integer :: ierr_u "); 488 kppi.add_line(" integer :: istatf "); 489 kppi.add_line(" integer :: vl_dim_lo "); 490 kppi.add_line(" "); 491 kppi.add_line(" "); 492 kppi.add_line(" if (present (istatus) ) istatus = 0 "); 493 kppi.add_line(" if (present (icntrl_i) ) icntrl = icntrl_i "); 494 kppi.add_line(" if (present (rcntrl_i) ) rcntrl = rcntrl_i "); 495 kppi.add_line(" "); 496 kppi.add_line(" vl_glo = size(tempi,1) "); 497 kppi.add_line(" "); 498 kppi.add_line(" vl_dim_lo = VL_DIM "); 499 kppi.add_line(" DO k=1,VL_glo,vl_dim_lo "); 498 500 kppi.add_line(" is = k "); 499 kppi.add_line(" ie = min(k+vl_dim -1,vl_glo)");501 kppi.add_line(" ie = min(k+vl_dim_lo-1,VL_glo) "); 500 502 kppi.add_line(" vl = ie-is+1 "); 501 503 502 504 kppi.add_line(" "); 503 505 if(kpp_switches.is_vector()) { 504 kppi.add_line(" c(1:vl,:) = conc(is:ie,:) "); 505 } else { 506 kppi.add_line(" c(:) = conc(is,:) "); 507 } 508 kppi.add_line(" "); 509 if(kpp_switches.is_vector()) { 510 kppi.add_line(" temp(1:vl) = tempk(is:ie) "); 511 } else { 512 kppi.add_line(" temp = tempk(is) "); 513 } 506 kppi.add_line(" C(1:vl,:) = Conc(is:ie,:) "); 507 } else { 508 kppi.add_line(" C(:) = Conc(is,:) "); 509 } 510 511 kppi.add_line(" "); 512 if(kpp_switches.is_vector()) { 513 kppi.add_line(" temp(1:vl) = tempi(is:ie) "); 514 } else { 515 kppi.add_line(" temp = tempi(is) "); 516 } 517 kppi.add_line(" "); 518 if(kpp_switches.is_vector()) { 519 kppi.add_line(" qvap(1:vl) = qvapi(is:ie) "); 520 } else { 521 kppi.add_line(" qvap = qvapi(is) "); 522 } 523 kppi.add_line(" "); 524 if(kpp_switches.is_vector()) { 525 kppi.add_line(" fakt(1:vl) = fakti(is:ie) "); 526 } else { 527 kppi.add_line(" fakt = fakti(is) "); 528 } 529 530 kppi.add_line(" "); 531 kppi.add_line(" CALL initialize "); 532 514 533 kppi.add_line(" "); 515 534 if(kpp_switches.is_vector()) { … … 528 547 kppi.add_line(" "); 529 548 if(kpp_switches.is_vector()) { 530 kppi.add_line(" conc(is:ie,:) = c(1:VL,:) "); 531 } else { 532 kppi.add_line(" IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN "); 533 kppi.add_line(" IF (l_debug) CALL error_output(conc(is,:),ierr_u, pe) "); 534 kppi.add_line(" ENDIF "); 535 kppi.add_line(" conc(is,:) = c(:) "); 549 kppi.add_line(" Conc(is:ie,:) = C(1:VL,:) "); 550 } else { 551 kppi.add_line(" IF (PRESENT(l_debug) .AND. PRESENT(PE)) THEN "); 552 kppi.add_line(" IF (l_debug) CALL error_output(Conc(is,:),ierr_u, PE) "); 553 kppi.add_line(" ENDIF "); 554 kppi.add_line(" "); 555 kppi.add_line(" Conc(is,:) = C(:) "); 536 556 } 537 557 … … 540 560 kppi.add_line(" "); 541 561 if(kpp_switches.is_vector()) { 542 kppi.add_line(" if(P RESENT(ierrf)) ierrf(is:ie) = ierrv(1:vl) ");543 kppi.add_line(" if(P RESENT(xnacc)) xnacc(is:ie) = kacc(1:vl) ");544 kppi.add_line(" if(P RESENT(xnrej)) xnrej(is:ie) = krej(1:vl) ");545 } else { 546 kppi.add_line(" if(P RESENT(ierrf)) ierrf(is) = ierr_u");547 kppi.add_line(" if(P RESENT(xnacc)) xnacc(is) = istatus_u(4) ");548 kppi.add_line(" if(P RESENT(xnrej)) xnrej(is) = istatus_u(5) ");549 } 550 kppi.add_line(" "); 551 kppi.add_line(" if ( PRESENT(istatus) ) then ");552 if(kpp_switches.is_vector()) { 553 kppi.add_line(" istatus(4) = istatus(4) + sum( kacc(1:vl)) ");554 kppi.add_line(" istatus(5) = istatus(5) + sum( krej(1:vl)) ");562 kppi.add_line(" if(Present(ierrf)) ierrf(is:ie) = IERRV(1:VL) "); 563 kppi.add_line(" if(Present(xNacc)) xNacc(is:ie) = Kacc(1:VL) "); 564 kppi.add_line(" if(Present(xNrej)) xNrej(is:ie) = Krej(1:VL) "); 565 } else { 566 kppi.add_line(" if(Present(ierrf)) ierrf(is) = IERR_U "); 567 kppi.add_line(" if(Present(xNacc)) xNacc(is) = istatus_u(4) "); 568 kppi.add_line(" if(Present(xNrej)) xNrej(is) = istatus_u(5) "); 569 } 570 kppi.add_line(" "); 571 kppi.add_line(" if (present (istatus) ) then "); 572 if(kpp_switches.is_vector()) { 573 kppi.add_line(" istatus(4) = istatus(4) + sum(Kacc(1:VL)) "); 574 kppi.add_line(" istatus(5) = istatus(5) + sum(Krej(1:VL)) "); 555 575 kppi.add_line(" istatus(3) = istatus(4) + istatus(5) "); 556 576 kppi.add_line(" istatus(6) = istatus(6) + istatus_u(6) "); … … 559 579 kppi.add_line(" istatus(1:8) = istatus(1:8) + istatus_u(1:8) "); 560 580 } 561 kppi.add_line(" ENDIF");581 kppi.add_line(" end if "); 562 582 kppi.add_line(" "); 563 583 kppi.add_line(" END DO "); … … 568 588 kppi.add_line(" "); 569 589 for(iv=Vvar_list.begin();iv!=Vvar_list.end();iv++) { 570 kppi.add_line(" if (ALLOCATED("+ iv->name +")) DEALLOCATE("+ iv->name +" ) ");590 // kppi.add_line(" if (allocated("+ iv->name +")) deallocate("+ iv->name +" ) "); 571 591 } 572 592 … … 574 594 kppi.add_line(" data_loaded = .false. "); 575 595 kppi.add_line(" "); 576 kppi.add_line(" RETURN");596 kppi.add_line(" return "); 577 597 kppi.add_line("END SUBROUTINE chem_gasphase_integrate "); 578 598 … … 583 603 } 584 604 585 void create_kpp_module::create_fill_routine(vector<fortran_file> &fi_list, Vvar & var) {586 fortran_file fi;587 vector<string>::iterator is;588 string line;589 590 cout << "Generate fill subroutine for " << var.name << endl;591 592 fi.set_name(var.name);593 line = " SUBROUTINE fill_" + var.name;594 fi.add_line(line + "(status, array) ");595 fi.add_line(" ");596 fi.add_line(" INTEGER, INTENT(OUT) :: status ");597 if(var.nr_dim() == 0) {598 fi.add_line(" REAL(dp), INTENT(IN), DIMENSION(:) :: array ");599 fi.add_line(" ");600 fi.add_line(" status = 0");601 fi.add_line(" IF (.not. ALLOCATED("+var.name+")) & ");602 fi.add_line(" ALLOCATE("+var.name+"(size(array))) ");603 } else if(var.nr_dim() == 1) {604 fi.add_line(" REAL (dp), INTENT(IN), DIMENSION(:,:) :: array ");605 fi.add_line(" ");606 fi.add_line(" status = 0 ");607 fi.add_line(" if (.not. ALLOCATED("+var.name+")) & ");608 fi.add_line(" ALLOCATE("+var.name+"(size(array,1),"+var.dim_var[0]+")) ");609 } else if(var.nr_dim() == 2) {610 fi.add_line(" ");611 fi.add_line(" REAL (dp), INTENT(IN), DIMENSION(:,:,:) :: array ");612 fi.add_line(" ");613 fi.add_line(" status = 0 ");614 fi.add_line(" if (.not. ALLOCATED("+var.name+")) & ");615 fi.add_line(" ALLOCATE("+var.name+"(size(array,1),"+var.dim_var[0]+var.dim_var[1]+")) ");616 } else {617 fi.add_line(" REAL (dp), INTENT(IN), DIMENSION(:,:,:,:) :: array ");618 fi.add_line(" ");619 fi.add_line(" status = 0 ");620 fi.add_line(" IF (.not. ALLOCATED("+var.name+")) & ");621 fi.add_line(" ALLOCATE("+var.name+"(size(array,1),"+var.dim_var[0]622 +var.dim_var[1]+var.dim_var[3]+")) ");623 }624 625 fi.add_line(" ");626 fi.add_line(" IF (data_loaded .AND. (vl_glo /= size(array,1)) ) THEN ");627 fi.add_line(" status = 1 ");628 fi.add_line(" RETURN ");629 fi.add_line(" END IF ");630 fi.add_line(" ");631 fi.add_line(" vl_glo = size(array,1) ");632 fi.add_line(" "+var.name+ " = array ");633 fi.add_line(" data_loaded = .TRUE. ");634 fi.add_line(" ");635 fi.add_line(" RETURN");636 fi.add_line(" ");637 fi.add_line(" END " + line);638 639 fi_list.push_back(fi);640 641 return;642 }
Note: See TracChangeset
for help on using the changeset viewer.
