source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_phstat/chem_gasphase_mod.f90 @ 4901

Last change on this file since 4901 was 4841, checked in by forkel, 19 months ago

updated copyright statements for chem_gasphase_mod.f90. This must be done in UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header and NOT in SOURCE/chem_gasphase_mod.f90.

File size: 81.2 KB
Line 
1MODULE chem_gasphase_mod
2 
3!   Mechanism: phstat
4!
5!------------------------------------------------------------------------------!
6!
7! ******Module chem_gasphase_mod is automatically generated by kpp4palm ******
8!
9!   *********Please do NOT change this Code,it will be ovewritten *********
10!
11!------------------------------------------------------------------------------!
12! This file was created by KPP (http://people.cs.vt.edu/asandu/Software/Kpp/)
13! and kpp4palm (created by Klaus Ketelsen). kpp4palm is an adapted version
14! of KP4 (Jöckel,P.,Kerkweg,A.,Pozzer,A.,Sander,R.,Tost,H.,Riede,
15! H.,Baumgaertner,A.,Gromov,S.,and Kern,B.,2010: Development cycle 2 of
16! the Modular Earth Submodel System (MESSy2),Geosci. Model Dev.,3,717-752,
17! https://doi.org/10.5194/gmd-3-717-2010). KP4 is part of the Modular Earth
18! Submodel System (MESSy),which is is available under the  GNU General Public
19! License (GPL).
20!
21! KPP is free software; you can redistribute it and/or modify it under the terms
22! of the General Public Licence as published by the Free Software Foundation;
23! either version 2 of the License,or (at your option) any later version.
24! KPP is distributed in the hope that it will be useful,but WITHOUT ANY WARRANTY;
25! without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
26! PURPOSE. See the GNU General Public Licence for more details.
27!
28!------------------------------------------------------------------------------!
29! This file is part of the PALM model system.
30!
31! PALM is free software: you can redistribute it and/or modify it under the
32! terms of the GNU General Public License as published by the Free Software
33! Foundation,either version 3 of the License,or (at your option) any later
34! version.
35!
36! PALM is distributed in the hope that it will be useful,but WITHOUT ANY
37! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
38! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
39!
40! You should have received a copy of the GNU General Public License along with
41! PALM. If not,see <http://www.gnu.org/licenses/>.
42!
43! Copyright 2017-2021 Leibniz Universitaet Hannover
44! Copyright 2016-2021 Klaus Ketelsen
45! Copyright 2017-2021 Karlsruhe Institute of Technology
46!--------------------------------------------------------------------------------!
47!
48!
49! MODULE HEADER TEMPLATE
50!
51!  Initial version (Nov. 2016,ketelsen),for later modifications of module_header
52!  see comments in kpp4palm/src/create_kpp_module.C
53
54! Set kpp Double Precision to PALM Default Precision
55
56  USE kinds,           ONLY: dp=>wp
57
58  USE pegrid,          ONLY: myid, threads_per_task
59
60  IMPLICIT        NONE
61  PRIVATE
62  !SAVE  ! note: occurs again in automatically generated code ...
63
64! Public variables
65  PUBLIC :: atol
66  PUBLIC :: cs_mech
67  PUBLIC :: eqn_names
68  PUBLIC :: fakt
69  PUBLIC :: nmaxfixsteps
70  PUBLIC :: nphot
71  PUBLIC :: nreact
72  PUBLIC :: nspec
73  PUBLIC :: nvar
74  PUBLIC :: qvap
75  PUBLIC :: phot
76  PUBLIC :: phot_names
77  PUBLIC :: rconst
78  PUBLIC :: rtol
79  PUBLIC :: spc_names
80  PUBLIC :: temp
81  PUBLIC :: vl_dim                     !< PUBLIC to enable other MODULEs to distiguish between scalar and vec
82 
83! Public routines
84  PUBLIC :: chem_gasphase_integrate
85  PUBLIC :: get_mechanism_name
86  PUBLIC :: initialize
87  PUBLIC :: initialize_kpp_ctrl
88  PUBLIC :: integrate
89  PUBLIC :: update_rconst
90
91! END OF MODULE HEADER TEMPLATE
92                                                                 
93! Variables used for vector mode                                 
94                                                                 
95  LOGICAL, PARAMETER          :: l_vector = .FALSE.           
96  INTEGER, PARAMETER          :: i_lu_di = 2
97  INTEGER, PARAMETER          :: vl_dim = 1
98  INTEGER                     :: vl                             
99                                                                 
100  INTEGER                     :: vl_glo                         
101  INTEGER                     :: is, ie                           
102                                                                 
103                                                                 
104  LOGICAL                     :: data_loaded = .FALSE.             
105  REAL(dp), POINTER, DIMENSION(:), CONTIGUOUS    :: var             
106! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
107!
108! Parameter Module File
109!
110! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
111!       (http://www.cs.vt.edu/~asandu/Software/KPP)
112! KPP is distributed under GPL,the general public licence
113!       (http://www.gnu.org/copyleft/gpl.html)
114! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
115! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
116!     With important contributions from:
117!        M. Damian,Villanova University,USA
118!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
119!
120! File                 : chem_gasphase_mod_Parameters.f90
121! Time                 : Tue Jan 12 18:41:28 2021
122! Working directory    : /home/forkel-r/palmstuff/work/trunk20210112/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
123! Equation file        : chem_gasphase_mod.kpp
124! Output root filename : chem_gasphase_mod
125!
126! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
127
128
129
130
131
132
133! NSPEC - Number of chemical species
134  INTEGER, PARAMETER :: nspec = 3 
135! NVAR - Number of Variable species
136  INTEGER, PARAMETER :: nvar = 3 
137! NVARACT - Number of Active species
138  INTEGER, PARAMETER :: nvaract = 3 
139! NFIX - Number of Fixed species
140  INTEGER, PARAMETER :: nfix = 1 
141! NREACT - Number of reactions
142  INTEGER, PARAMETER :: nreact = 2 
143! NVARST - Starting of variables in conc. vect.
144  INTEGER, PARAMETER :: nvarst = 1 
145! NFIXST - Starting of fixed in conc. vect.
146  INTEGER, PARAMETER :: nfixst = 4 
147! NONZERO - Number of nonzero entries in Jacobian
148  INTEGER, PARAMETER :: nonzero = 9 
149! LU_NONZERO - Number of nonzero entries in LU factoriz. of Jacobian
150  INTEGER, PARAMETER :: lu_nonzero = 9 
151! CNVAR - (NVAR+1) Number of elements in compressed row format
152  INTEGER, PARAMETER :: cnvar = 4 
153! CNEQN - (NREACT+1) Number stoicm elements in compressed col format
154  INTEGER, PARAMETER :: cneqn = 3 
155! NHESS - Length of Sparse Hessian
156  INTEGER, PARAMETER :: nhess = 3 
157! NMASS - Number of atoms to check mass balance
158  INTEGER, PARAMETER :: nmass = 1 
159
160! Index declaration for variable species in C and VAR
161!   VAR(ind_spc) = C(ind_spc)
162
163  INTEGER, PARAMETER, PUBLIC :: ind_o3 = 1 
164  INTEGER, PARAMETER, PUBLIC :: ind_no = 2 
165  INTEGER, PARAMETER, PUBLIC :: ind_no2 = 3 
166
167! Index declaration for fixed species in C
168!   C(ind_spc)
169
170
171! Index declaration for fixed species in FIX
172!    FIX(indf_spc) = C(ind_spc) = C(NVAR+indf_spc)
173
174
175! NJVRP - Length of sparse Jacobian JVRP
176  INTEGER, PARAMETER :: njvrp = 3 
177
178! NSTOICM - Length of Sparse Stoichiometric Matrix
179  INTEGER, PARAMETER :: nstoicm = 6 
180
181
182! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
183!
184! Global Data Module File
185!
186! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
187!       (http://www.cs.vt.edu/~asandu/Software/KPP)
188! KPP is distributed under GPL,the general public licence
189!       (http://www.gnu.org/copyleft/gpl.html)
190! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
191! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
192!     With important contributions from:
193!        M. Damian,Villanova University,USA
194!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
195!
196! File                 : chem_gasphase_mod_Global.f90
197! Time                 : Tue Jan 12 18:41:28 2021
198! Working directory    : /home/forkel-r/palmstuff/work/trunk20210112/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
199! Equation file        : chem_gasphase_mod.kpp
200! Output root filename : chem_gasphase_mod
201!
202! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
203
204
205
206
207
208
209! Declaration of global variables
210
211! C - Concentration of all species
212  REAL(kind=dp), TARGET    :: c(nspec)
213! VAR - Concentrations of variable species (global)
214! REAL(kind=dp):: var(nvar)  var is now POINTER
215! FIX - Concentrations of fixed species (global)
216  REAL(kind=dp):: fix(nfix)
217! VAR,FIX are chunks of array C
218! RCONST - Rate constants (global)
219  REAL(kind=dp):: rconst(nreact)
220! TIME - Current integration time
221  REAL(kind=dp):: time
222! TEMP - Temperature
223  REAL(kind=dp):: temp
224! ATOL - Absolute tolerance
225  REAL(kind=dp):: atol(nvar)
226! RTOL - Relative tolerance
227  REAL(kind=dp):: rtol(nvar)
228! STEPMIN - Lower bound for integration step
229  REAL(kind=dp):: stepmin
230! CFACTOR - Conversion factor for concentration units
231  REAL(kind=dp):: cfactor
232
233! INLINED global variable declarations
234
235! QVAP - Water vapor
236  REAL(kind=dp):: qvap
237! FAKT - Conversion factor
238  REAL(kind=dp):: fakt
239
240! CS_MECH for check of mechanism name with namelist
241  CHARACTER(len=30):: cs_mech
242
243! INLINED global variable declarations
244
245
246
247! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
248!
249! Sparse Jacobian Data Structures File
250!
251! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
252!       (http://www.cs.vt.edu/~asandu/Software/KPP)
253! KPP is distributed under GPL,the general public licence
254!       (http://www.gnu.org/copyleft/gpl.html)
255! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
256! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
257!     With important contributions from:
258!        M. Damian,Villanova University,USA
259!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
260!
261! File                 : chem_gasphase_mod_JacobianSP.f90
262! Time                 : Tue Jan 12 18:41:28 2021
263! Working directory    : /home/forkel-r/palmstuff/work/trunk20210112/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
264! Equation file        : chem_gasphase_mod.kpp
265! Output root filename : chem_gasphase_mod
266!
267! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
268
269
270
271
272
273
274! Sparse Jacobian Data
275
276
277  INTEGER, PARAMETER, DIMENSION(9):: lu_irow =  (/ &
278       1, 1, 1, 2, 2, 2, 3, 3, 3 /) 
279
280  INTEGER, PARAMETER, DIMENSION(9):: lu_icol =  (/ &
281       1, 2, 3, 1, 2, 3, 1, 2, 3 /) 
282
283  INTEGER, PARAMETER, DIMENSION(4):: lu_crow =  (/ &
284       1, 4, 7, 10 /) 
285
286  INTEGER, PARAMETER, DIMENSION(4):: lu_diag =  (/ &
287       1, 5, 9, 10 /) 
288
289
290
291! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
292!
293! Utility Data Module File
294!
295! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
296!       (http://www.cs.vt.edu/~asandu/Software/KPP)
297! KPP is distributed under GPL,the general public licence
298!       (http://www.gnu.org/copyleft/gpl.html)
299! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
300! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
301!     With important contributions from:
302!        M. Damian,Villanova University,USA
303!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
304!
305! File                 : chem_gasphase_mod_Monitor.f90
306! Time                 : Tue Jan 12 18:41:28 2021
307! Working directory    : /home/forkel-r/palmstuff/work/trunk20210112/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
308! Equation file        : chem_gasphase_mod.kpp
309! Output root filename : chem_gasphase_mod
310!
311! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
312
313
314
315
316
317  CHARACTER(len=15), PARAMETER, DIMENSION(3):: spc_names =  (/ &
318     'O3             ','NO             ','NO2            ' /)
319
320  CHARACTER(len=100), PARAMETER, DIMENSION(2):: eqn_names =  (/ &
321     '    NO2 --> O3 + NO                                                                                 ',&
322     'O3 + NO --> NO2                                                                                     ' /)
323
324! INLINED global variables
325
326  !   inline f90_data: declaration of global variables for photolysis
327  !   REAL(kind=dp):: phot(nphot)must eventually be moved to global later for
328  INTEGER, PARAMETER :: nphot = 1
329  !   phot photolysis frequencies
330  REAL(kind=dp):: phot(nphot)
331
332  INTEGER, PARAMETER, PUBLIC :: j_no2 = 1
333
334  CHARACTER(len=15), PARAMETER, DIMENSION(nphot):: phot_names =   (/ &
335     'J_NO2          '/)
336
337! End INLINED global variables
338
339
340! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
341 
342! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
343 
344! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
345 
346! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
347 
348 
349!  variable definations from  individual module headers
350 
351! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
352!
353! Initialization File
354!
355! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
356!       (http://www.cs.vt.edu/~asandu/Software/KPP)
357! KPP is distributed under GPL,the general public licence
358!       (http://www.gnu.org/copyleft/gpl.html)
359! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
360! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
361!     With important contributions from:
362!        M. Damian,Villanova University,USA
363!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
364!
365! File                 : chem_gasphase_mod_Initialize.f90
366! Time                 : Tue Jan 12 18:41:28 2021
367! Working directory    : /home/forkel-r/palmstuff/work/trunk20210112/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
368! Equation file        : chem_gasphase_mod.kpp
369! Output root filename : chem_gasphase_mod
370!
371! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
372
373
374
375
376
377! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
378!
379! Numerical Integrator (Time-Stepping) File
380!
381! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
382!       (http://www.cs.vt.edu/~asandu/Software/KPP)
383! KPP is distributed under GPL,the general public licence
384!       (http://www.gnu.org/copyleft/gpl.html)
385! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
386! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
387!     With important contributions from:
388!        M. Damian,Villanova University,USA
389!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
390!
391! File                 : chem_gasphase_mod_Integrator.f90
392! Time                 : Tue Jan 12 18:41:28 2021
393! Working directory    : /home/forkel-r/palmstuff/work/trunk20210112/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
394! Equation file        : chem_gasphase_mod.kpp
395! Output root filename : chem_gasphase_mod
396!
397! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
398
399
400
401! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
402!
403! INTEGRATE - Integrator routine
404!   Arguments :
405!      TIN       - Start Time for Integration
406!      TOUT      - End Time for Integration
407!
408! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
409
410!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
411!  Rosenbrock - Implementation of several Rosenbrock methods:             !
412!               *Ros2                                                    !
413!               *Ros3                                                    !
414!               *Ros4                                                    !
415!               *Rodas3                                                  !
416!               *Rodas4                                                  !
417!  By default the code employs the KPP sparse linear algebra routines     !
418!  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
419!                                                                         !
420!    (C)  Adrian Sandu,August 2004                                       !
421!    Virginia Polytechnic Institute and State University                  !
422!    Contact: sandu@cs.vt.edu                                             !
423!    Revised by Philipp Miehe and Adrian Sandu,May 2006                  !                               !
424!    This implementation is part of KPP - the Kinetic PreProcessor        !
425!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
426
427
428  SAVE
429 
430!~~~>  statistics on the work performed by the rosenbrock method
431  INTEGER, PARAMETER :: nfun=1, njac=2, nstp=3, nacc=4, &
432                        nrej=5, ndec=6, nsol=7, nsng=8, &
433                        ntexit=1, nhexit=2, nhnew = 3
434
435! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
436!
437! Linear Algebra Data and Routines File
438!
439! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
440!       (http://www.cs.vt.edu/~asandu/Software/KPP)
441! KPP is distributed under GPL,the general public licence
442!       (http://www.gnu.org/copyleft/gpl.html)
443! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
444! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
445!     With important contributions from:
446!        M. Damian,Villanova University,USA
447!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
448!
449! File                 : chem_gasphase_mod_LinearAlgebra.f90
450! Time                 : Tue Jan 12 18:41:28 2021
451! Working directory    : /home/forkel-r/palmstuff/work/trunk20210112/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
452! Equation file        : chem_gasphase_mod.kpp
453! Output root filename : chem_gasphase_mod
454!
455! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
456
457
458
459
460
461
462! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
463!
464! The ODE Jacobian of Chemical Model File
465!
466! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
467!       (http://www.cs.vt.edu/~asandu/Software/KPP)
468! KPP is distributed under GPL,the general public licence
469!       (http://www.gnu.org/copyleft/gpl.html)
470! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
471! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
472!     With important contributions from:
473!        M. Damian,Villanova University,USA
474!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
475!
476! File                 : chem_gasphase_mod_Jacobian.f90
477! Time                 : Tue Jan 12 18:41:28 2021
478! Working directory    : /home/forkel-r/palmstuff/work/trunk20210112/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
479! Equation file        : chem_gasphase_mod.kpp
480! Output root filename : chem_gasphase_mod
481!
482! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
483
484
485
486
487
488
489! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
490!
491! The ODE Function of Chemical Model File
492!
493! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
494!       (http://www.cs.vt.edu/~asandu/Software/KPP)
495! KPP is distributed under GPL,the general public licence
496!       (http://www.gnu.org/copyleft/gpl.html)
497! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
498! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
499!     With important contributions from:
500!        M. Damian,Villanova University,USA
501!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
502!
503! File                 : chem_gasphase_mod_Function.f90
504! Time                 : Tue Jan 12 18:41:28 2021
505! Working directory    : /home/forkel-r/palmstuff/work/trunk20210112/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
506! Equation file        : chem_gasphase_mod.kpp
507! Output root filename : chem_gasphase_mod
508!
509! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
510
511
512
513
514
515! A - Rate for each equation
516  REAL(kind=dp):: a(nreact)
517
518! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
519!
520! The Reaction Rates File
521!
522! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
523!       (http://www.cs.vt.edu/~asandu/Software/KPP)
524! KPP is distributed under GPL,the general public licence
525!       (http://www.gnu.org/copyleft/gpl.html)
526! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
527! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
528!     With important contributions from:
529!        M. Damian,Villanova University,USA
530!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
531!
532! File                 : chem_gasphase_mod_Rates.f90
533! Time                 : Tue Jan 12 18:41:28 2021
534! Working directory    : /home/forkel-r/palmstuff/work/trunk20210112/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
535! Equation file        : chem_gasphase_mod.kpp
536! Output root filename : chem_gasphase_mod
537!
538! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
539
540
541
542
543
544! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
545!
546! Auxiliary Routines File
547!
548! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
549!       (http://www.cs.vt.edu/~asandu/Software/KPP)
550! KPP is distributed under GPL,the general public licence
551!       (http://www.gnu.org/copyleft/gpl.html)
552! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
553! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
554!     With important contributions from:
555!        M. Damian,Villanova University,USA
556!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
557!
558! File                 : chem_gasphase_mod_Util.f90
559! Time                 : Tue Jan 12 18:41:28 2021
560! Working directory    : /home/forkel-r/palmstuff/work/trunk20210112/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
561! Equation file        : chem_gasphase_mod.kpp
562! Output root filename : chem_gasphase_mod
563!
564! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
565
566
567
568
569
570
571  ! header MODULE initialize_kpp_ctrl_template
572
573  ! notes:
574  ! - l_vector is automatically defined by kp4
575  ! - vl_dim is automatically defined by kp4
576  ! - i_lu_di is automatically defined by kp4
577  ! - wanted is automatically defined by xmecca
578  ! - icntrl rcntrl are automatically defined by kpp
579  ! - "USE messy_main_tools" is in MODULE_header of messy_mecca_kpp.f90
580  ! - SAVE will be automatically added by kp4
581
582  !SAVE
583
584  ! for fixed time step control
585  ! ... max. number of fixed time steps (sum must be 1)
586  INTEGER, PARAMETER         :: nmaxfixsteps = 50
587  ! ... switch for fixed time stepping
588  LOGICAL, PUBLIC            :: l_fixed_step = .FALSE.
589  INTEGER, PUBLIC            :: nfsteps = 1
590  ! ... number of kpp control PARAMETERs
591  INTEGER, PARAMETER, PUBLIC :: nkppctrl = 20
592  !
593  ! steering PARAMETERs for chemistry solver (see kpp domumentation)
594  INTEGER, DIMENSION(nkppctrl), PUBLIC      :: icntrl = 0
595  REAL(dp), DIMENSION(nkppctrl), PUBLIC     :: rcntrl = 0.0_dp
596  ! t_steps: fixed time steps in vector mode
597  REAL(dp), DIMENSION(nmaxfixsteps), PUBLIC :: t_steps = 0.0_dp
598
599  ! END header MODULE initialize_kpp_ctrl_template
600
601 
602! Interface Block
603 
604  INTERFACE            initialize
605    MODULE PROCEDURE   initialize
606  END INTERFACE        initialize
607 
608  INTERFACE            integrate
609    MODULE PROCEDURE   integrate
610  END INTERFACE        integrate
611 
612  INTERFACE            fun
613    MODULE PROCEDURE   fun
614  END INTERFACE        fun
615 
616  INTERFACE            kppsolve
617    MODULE PROCEDURE   kppsolve
618  END INTERFACE        kppsolve
619 
620  INTERFACE            jac_sp
621    MODULE PROCEDURE   jac_sp
622  END INTERFACE        jac_sp
623 
624  INTERFACE            k_arr
625    MODULE PROCEDURE   k_arr
626  END INTERFACE        k_arr
627 
628  INTERFACE            update_rconst
629    MODULE PROCEDURE   update_rconst
630  END INTERFACE        update_rconst
631 
632  INTERFACE            arr2
633    MODULE PROCEDURE   arr2
634  END INTERFACE        arr2
635 
636  INTERFACE            initialize_kpp_ctrl
637    MODULE PROCEDURE   initialize_kpp_ctrl
638  END INTERFACE        initialize_kpp_ctrl
639 
640  INTERFACE            error_output
641    MODULE PROCEDURE   error_output
642  END INTERFACE        error_output
643 
644  INTERFACE            wscal
645    MODULE PROCEDURE   wscal
646  END INTERFACE        wscal
647 
648!INTERFACE not working  INTERFACE            waxpy
649!INTERFACE not working    MODULE PROCEDURE   waxpy
650!INTERFACE not working  END INTERFACE        waxpy
651 
652  INTERFACE            rosenbrock
653    MODULE PROCEDURE   rosenbrock
654  END INTERFACE        rosenbrock
655 
656  INTERFACE            funtemplate
657    MODULE PROCEDURE   funtemplate
658  END INTERFACE        funtemplate
659 
660  INTERFACE            jactemplate
661    MODULE PROCEDURE   jactemplate
662  END INTERFACE        jactemplate
663 
664  INTERFACE            kppdecomp
665    MODULE PROCEDURE   kppdecomp
666  END INTERFACE        kppdecomp
667 
668  INTERFACE            get_mechanism_name
669    MODULE PROCEDURE   get_mechanism_name
670  END INTERFACE        get_mechanism_name
671 
672  INTERFACE            chem_gasphase_integrate
673    MODULE PROCEDURE   chem_gasphase_integrate
674  END INTERFACE        chem_gasphase_integrate
675 
676 
677  ! openmp directives generated by kp4
678 
679  !$OMP THREADPRIVATE (vl, vl_glo, is, ie, data_loaded)
680  !$OMP THREADPRIVATE (c, var, fix, rconst, time, temp, stepmin, cfactor)
681  !$OMP THREADPRIVATE (qvap, fakt, cs_mech, a, icntrl, rcntrl)
682 
683 CONTAINS
684 
685SUBROUTINE initialize()
686
687
688 INTEGER         :: k
689
690  INTEGER :: i
691  REAL(kind=dp):: x
692  k = is
693  cfactor = 1.000000e+00_dp
694
695! Following line is just to avoid compiler message about unused variables
696   IF ( lu_crow(1) == 1  .OR.  lu_icol(1) == 1  .OR.  lu_irow(1) == 1 )  CONTINUE 
697
698
699  x = (0.) * cfactor
700  DO i = 1 , nvar
701  ENDDO
702
703  x = (0.) * cfactor
704  DO i = 1 , nfix
705    fix(i) = x
706  ENDDO
707
708! constant rate coefficients
709! END constant rate coefficients
710
711! INLINED initializations
712
713! End INLINED initializations
714
715     
716END SUBROUTINE initialize
717 
718SUBROUTINE integrate( tin, tout, &
719  icntrl_u, rcntrl_u, istatus_u, rstatus_u, ierr_u)
720
721
722   REAL(kind=dp), INTENT(IN):: tin  ! start time
723   REAL(kind=dp), INTENT(IN):: tout ! END time
724   ! OPTIONAL input PARAMETERs and statistics
725   INTEGER,      INTENT(IN), OPTIONAL :: icntrl_u(20)
726   REAL(kind=dp), INTENT(IN), OPTIONAL :: rcntrl_u(20)
727   INTEGER,      INTENT(OUT), OPTIONAL :: istatus_u(20)
728   REAL(kind=dp), INTENT(OUT), OPTIONAL :: rstatus_u(20)
729   INTEGER,      INTENT(OUT), OPTIONAL :: ierr_u
730
731   REAL(kind=dp):: rcntrl(20), rstatus(20)
732   INTEGER       :: icntrl(20), istatus(20), ierr
733
734
735   icntrl(:) = 0
736   rcntrl(:) = 0.0_dp
737   istatus(:) = 0
738   rstatus(:) = 0.0_dp
739
740    !~~~> fine-tune the integrator:
741   icntrl(1) = 0      ! 0 - non- autonomous, 1 - autonomous
742   icntrl(2) = 0      ! 0 - vector tolerances, 1 - scalars
743
744   ! IF OPTIONAL PARAMETERs are given, and IF they are >0,
745   ! THEN they overwrite default settings.
746   IF (PRESENT(icntrl_u))THEN
747     WHERE(icntrl_u(:)> 0)icntrl(:) = icntrl_u(:)
748   ENDIF
749   IF (PRESENT(rcntrl_u))THEN
750     WHERE(rcntrl_u(:)> 0)rcntrl(:) = rcntrl_u(:)
751   ENDIF
752
753
754   CALL rosenbrock(nvar, var, tin, tout,  &
755         atol, rtol,               &
756         rcntrl, icntrl, rstatus, istatus, ierr)
757
758   !~~~> debug option: show no of steps
759   ! ntotal = ntotal + istatus(nstp)
760   ! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')','  O3=',VAR(ind_O3)
761
762   stepmin = rstatus(nhexit)
763   ! IF OPTIONAL PARAMETERs are given for output they
764   ! are updated with the RETURN information
765   IF (PRESENT(istatus_u))istatus_u(:) = istatus(:)
766   IF (PRESENT(rstatus_u))rstatus_u(:) = rstatus(:)
767   IF (PRESENT(ierr_u))  ierr_u       = ierr
768
769END SUBROUTINE integrate
770 
771SUBROUTINE fun(v, f, rct, vdot)
772
773! V - Concentrations of variable species (local)
774  REAL(kind=dp):: v(nvar)
775! F - Concentrations of fixed species (local)
776  REAL(kind=dp):: f(nfix)
777! RCT - Rate constants (local)
778  REAL(kind=dp):: rct(nreact)
779! Vdot - Time derivative of variable species concentrations
780  REAL(kind=dp):: vdot(nvar)
781
782
783! Following line is just to avoid compiler message about unused variables
784   IF ( f(nfix) > 0.0_dp )  CONTINUE
785
786! Computation of equation rates
787  a(1) = rct(1) * v(3)
788  a(2) = rct(2) * v(1) * v(2)
789
790! Aggregate function
791  vdot(1) = a(1) - a(2)
792  vdot(2) = a(1) - a(2)
793  vdot(3) = - a(1) + a(2)
794     
795END SUBROUTINE fun
796 
797SUBROUTINE kppsolve(jvs, x)
798
799! JVS - sparse Jacobian of variables
800  REAL(kind=dp):: jvs(lu_nonzero)
801! X - Vector for variables
802  REAL(kind=dp):: x(nvar)
803
804  x(2) = x(2) - jvs(4) * x(1)
805  x(3) = x(3) - jvs(7) * x(1) - jvs(8) * x(2)
806  x(3) = x(3) / jvs(9)
807  x(2) = (x(2) - jvs(6) * x(3)) /(jvs(5))
808  x(1) = (x(1) - jvs(2) * x(2) - jvs(3) * x(3)) /(jvs(1))
809     
810END SUBROUTINE kppsolve
811 
812SUBROUTINE jac_sp(v, f, rct, jvs)
813
814! V - Concentrations of variable species (local)
815  REAL(kind=dp):: v(nvar)
816! F - Concentrations of fixed species (local)
817  REAL(kind=dp):: f(nfix)
818! RCT - Rate constants (local)
819  REAL(kind=dp):: rct(nreact)
820! JVS - sparse Jacobian of variables
821  REAL(kind=dp):: jvs(lu_nonzero)
822
823
824! Local variables
825! B - Temporary array
826  REAL(kind=dp):: b(3)
827!
828! Following line is just to avoid compiler message about unused variables
829   IF ( f(nfix) > 0.0_dp )  CONTINUE
830
831! B(1) = dA(1)/dV(3)
832  b(1) = rct(1)
833! B(2) = dA(2)/dV(1)
834  b(2) = rct(2) * v(2)
835! B(3) = dA(2)/dV(2)
836  b(3) = rct(2) * v(1)
837
838! Construct the Jacobian terms from B's
839! JVS(1) = Jac_FULL(1,1)
840  jvs(1) = - b(2)
841! JVS(2) = Jac_FULL(1,2)
842  jvs(2) = - b(3)
843! JVS(3) = Jac_FULL(1,3)
844  jvs(3) = b(1)
845! JVS(4) = Jac_FULL(2,1)
846  jvs(4) = - b(2)
847! JVS(5) = Jac_FULL(2,2)
848  jvs(5) = - b(3)
849! JVS(6) = Jac_FULL(2,3)
850  jvs(6) = b(1)
851! JVS(7) = Jac_FULL(3,1)
852  jvs(7) = b(2)
853! JVS(8) = Jac_FULL(3,2)
854  jvs(8) = b(3)
855! JVS(9) = Jac_FULL(3,3)
856  jvs(9) = - b(1)
857     
858END SUBROUTINE jac_sp
859 
860  elemental REAL(kind=dp)FUNCTION k_arr (k_298, tdep, temp)
861    ! arrhenius FUNCTION
862
863    REAL,    INTENT(IN):: k_298 ! k at t = 298.15k
864    REAL,    INTENT(IN):: tdep  ! temperature dependence
865    REAL(kind=dp), INTENT(IN):: temp  ! temperature
866
867    intrinsic exp
868
869    k_arr = k_298 * exp(tdep* (1._dp/temp- 3.3540e-3_dp))! 1/298.15=3.3540e-3
870
871  END FUNCTION k_arr
872 
873SUBROUTINE update_rconst()
874 INTEGER         :: k
875
876  k = is
877
878! Begin INLINED RCONST
879
880
881! End INLINED RCONST
882
883  rconst(1) = (phot(j_no2))
884  rconst(2) = (arr2(1.8e-12_dp , 1370.0_dp , temp))
885     
886END SUBROUTINE update_rconst
887 
888!  END FUNCTION ARR2
889REAL(kind=dp)FUNCTION arr2( a0, b0, temp)
890    REAL(kind=dp):: temp
891    REAL(kind=dp):: a0, b0
892    arr2 = a0 * exp( - b0 / temp)
893END FUNCTION arr2
894 
895SUBROUTINE initialize_kpp_ctrl(status)
896
897
898  ! i/o
899  INTEGER,         INTENT(OUT):: status
900
901  ! local
902  REAL(dp):: tsum
903  INTEGER  :: i
904
905  ! check fixed time steps
906  tsum = 0.0_dp
907  DO i=1, nmaxfixsteps
908     IF (t_steps(i)< tiny(0.0_dp))exit
909     tsum = tsum + t_steps(i)
910  ENDDO
911
912  nfsteps = i- 1
913
914  l_fixed_step = (nfsteps > 0).and.((tsum - 1.0)< tiny(0.0_dp))
915
916  IF (l_vector)THEN
917     WRITE(*,*) ' MODE             : VECTOR (LENGTH=',VL_DIM,')'
918  ELSE
919     WRITE(*,*) ' MODE             : SCALAR'
920  ENDIF
921  !
922  WRITE(*,*) ' DE-INDEXING MODE :',I_LU_DI
923  !
924  WRITE(*,*) ' ICNTRL           : ',icntrl
925  WRITE(*,*) ' RCNTRL           : ',rcntrl
926  !
927  ! note: this is ONLY meaningful for vectorized (kp4)rosenbrock- methods
928  IF (l_vector)THEN
929     IF (l_fixed_step)THEN
930        WRITE(*,*) ' TIME STEPS       : FIXED (',t_steps(1:nfsteps),')'
931     ELSE
932        WRITE(*,*) ' TIME STEPS       : AUTOMATIC'
933     ENDIF
934  ELSE
935     WRITE(*,*) ' TIME STEPS       : AUTOMATIC '//&
936          &'(t_steps (CTRL_KPP) ignored in SCALAR MODE)'
937  ENDIF
938  ! mz_pj_20070531-
939
940  status = 0
941
942
943END SUBROUTINE initialize_kpp_ctrl
944 
945SUBROUTINE error_output(c, ierr, pe)
946
947
948  INTEGER, INTENT(IN):: ierr
949  INTEGER, INTENT(IN):: pe
950  REAL(dp), DIMENSION(:), INTENT(IN):: c
951
952  write(6,*) 'ERROR in chem_gasphase_mod ',ierr,C(1),PE
953
954
955END SUBROUTINE error_output
956 
957      SUBROUTINE wscal(n, alpha, x, incx)
958!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
959!     constant times a vector: x(1:N) <- Alpha*x(1:N)
960!     only for incX=incY=1
961!     after BLAS
962!     replace this by the function from the optimized BLAS implementation:
963!         CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
964!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
965
966      INTEGER  :: i, incx, m, mp1, n
967      REAL(kind=dp) :: x(n), alpha
968      REAL(kind=dp), PARAMETER  :: zero=0.0_dp, one=1.0_dp
969
970! Following line is just to avoid compiler message about unused variables
971   IF ( incx == 0 )  CONTINUE
972
973      IF (alpha .eq. one)RETURN
974      IF (n .le. 0)RETURN
975
976      m = mod(n, 5)
977      IF ( m .ne. 0)THEN
978        IF (alpha .eq. (- one))THEN
979          DO i = 1, m
980            x(i) = - x(i)
981          ENDDO
982        ELSEIF (alpha .eq. zero)THEN
983          DO i = 1, m
984            x(i) = zero
985          ENDDO
986        ELSE
987          DO i = 1, m
988            x(i) = alpha* x(i)
989          ENDDO
990        ENDIF
991        IF ( n .lt. 5)RETURN
992      ENDIF
993      mp1 = m + 1
994      IF (alpha .eq. (- one))THEN
995        DO i = mp1, n, 5
996          x(i)   = - x(i)
997          x(i + 1) = - x(i + 1)
998          x(i + 2) = - x(i + 2)
999          x(i + 3) = - x(i + 3)
1000          x(i + 4) = - x(i + 4)
1001        ENDDO
1002      ELSEIF (alpha .eq. zero)THEN
1003        DO i = mp1, n, 5
1004          x(i)   = zero
1005          x(i + 1) = zero
1006          x(i + 2) = zero
1007          x(i + 3) = zero
1008          x(i + 4) = zero
1009        ENDDO
1010      ELSE
1011        DO i = mp1, n, 5
1012          x(i)   = alpha* x(i)
1013          x(i + 1) = alpha* x(i + 1)
1014          x(i + 2) = alpha* x(i + 2)
1015          x(i + 3) = alpha* x(i + 3)
1016          x(i + 4) = alpha* x(i + 4)
1017        ENDDO
1018      ENDIF
1019
1020      END SUBROUTINE wscal
1021 
1022      SUBROUTINE waxpy(n, alpha, x, incx, y, incy)
1023!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1024!     constant times a vector plus a vector: y <- y + Alpha*x
1025!     only for incX=incY=1
1026!     after BLAS
1027!     replace this by the function from the optimized BLAS implementation:
1028!         CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
1029!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1030
1031      INTEGER  :: i, incx, incy, m, mp1, n
1032      REAL(kind=dp):: x(n), y(n), alpha
1033      REAL(kind=dp), PARAMETER :: zero = 0.0_dp
1034
1035
1036! Following line is just to avoid compiler message about unused variables
1037   IF ( incx == 0  .OR.  incy == 0 )  CONTINUE
1038      IF (alpha .eq. zero)RETURN
1039      IF (n .le. 0)RETURN
1040
1041      m = mod(n, 4)
1042      IF ( m .ne. 0)THEN
1043        DO i = 1, m
1044          y(i) = y(i) + alpha* x(i)
1045        ENDDO
1046        IF ( n .lt. 4)RETURN
1047      ENDIF
1048      mp1 = m + 1
1049      DO i = mp1, n, 4
1050        y(i) = y(i) + alpha* x(i)
1051        y(i + 1) = y(i + 1) + alpha* x(i + 1)
1052        y(i + 2) = y(i + 2) + alpha* x(i + 2)
1053        y(i + 3) = y(i + 3) + alpha* x(i + 3)
1054      ENDDO
1055     
1056      END SUBROUTINE waxpy
1057 
1058SUBROUTINE rosenbrock(n, y, tstart, tend, &
1059           abstol, reltol,             &
1060           rcntrl, icntrl, rstatus, istatus, ierr)
1061!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1062!
1063!    Solves the system y'=F(t,y) using a Rosenbrock method defined by:
1064!
1065!     G = 1/(H*gamma(1)) - Jac(t0,Y0)
1066!     T_i = t0 + Alpha(i)*H
1067!     Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
1068!     G *K_i = Fun( T_i,Y_i)+ \sum_{j=1}^S C(i,j)/H *K_j +
1069!         gamma(i)*dF/dT(t0,Y0)
1070!     Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
1071!
1072!    For details on Rosenbrock methods and their implementation consult:
1073!      E. Hairer and G. Wanner
1074!      "Solving ODEs II. Stiff and differential-algebraic problems".
1075!      Springer series in computational mathematics,Springer-Verlag,1996.
1076!    The codes contained in the book inspired this implementation.
1077!
1078!    (C)  Adrian Sandu,August 2004
1079!    Virginia Polytechnic Institute and State University
1080!    Contact: sandu@cs.vt.edu
1081!    Revised by Philipp Miehe and Adrian Sandu,May 2006                 
1082!    This implementation is part of KPP - the Kinetic PreProcessor
1083!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1084!
1085!~~~>   input arguments:
1086!
1087!-     y(n)  = vector of initial conditions (at t=tstart)
1088!-    [tstart, tend]  = time range of integration
1089!     (if Tstart>Tend the integration is performed backwards in time)
1090!-    reltol, abstol = user precribed accuracy
1091!- SUBROUTINE  fun( t, y, ydot) = ode FUNCTION,
1092!                       returns Ydot = Y' = F(T,Y)
1093!- SUBROUTINE  jac( t, y, jcb) = jacobian of the ode FUNCTION,
1094!                       returns Jcb = dFun/dY
1095!-    icntrl(1:20)  = INTEGER inputs PARAMETERs
1096!-    rcntrl(1:20)  = REAL inputs PARAMETERs
1097!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1098!
1099!~~~>     output arguments:
1100!
1101!-    y(n)  - > vector of final states (at t- >tend)
1102!-    istatus(1:20) - > INTEGER output PARAMETERs
1103!-    rstatus(1:20) - > REAL output PARAMETERs
1104!-    ierr            - > job status upon RETURN
1105!                        success (positive value) or
1106!                        failure (negative value)
1107!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1108!
1109!~~~>     input PARAMETERs:
1110!
1111!    Note: For input parameters equal to zero the default values of the
1112!       corresponding variables are used.
1113!
1114!    ICNTRL(1) = 1: F = F(y)   Independent of T (AUTONOMOUS)
1115!              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
1116!
1117!    ICNTRL(2) = 0: AbsTol,RelTol are N-dimensional vectors
1118!              = 1: AbsTol,RelTol are scalars
1119!
1120!    ICNTRL(3)  -> selection of a particular Rosenbrock method
1121!        = 0 :    Rodas3 (default)
1122!        = 1 :    Ros2
1123!        = 2 :    Ros3
1124!        = 3 :    Ros4
1125!        = 4 :    Rodas3
1126!        = 5 :    Rodas4
1127!
1128!    ICNTRL(4)  -> maximum number of integration steps
1129!        For ICNTRL(4) =0) the default value of 100000 is used
1130!
1131!    RCNTRL(1)  -> Hmin,lower bound for the integration step size
1132!          It is strongly recommended to keep Hmin = ZERO
1133!    RCNTRL(2)  -> Hmax,upper bound for the integration step size
1134!    RCNTRL(3)  -> Hstart,starting value for the integration step size
1135!
1136!    RCNTRL(4)  -> FacMin,lower bound on step decrease factor (default=0.2)
1137!    RCNTRL(5)  -> FacMax,upper bound on step increase factor (default=6)
1138!    RCNTRL(6)  -> FacRej,step decrease factor after multiple rejections
1139!                          (default=0.1)
1140!    RCNTRL(7)  -> FacSafe,by which the new step is slightly smaller
1141!         than the predicted value  (default=0.9)
1142!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1143!
1144!
1145!    OUTPUT ARGUMENTS:
1146!    -----------------
1147!
1148!    T           -> T value for which the solution has been computed
1149!                     (after successful return T=Tend).
1150!
1151!    Y(N)        -> Numerical solution at T
1152!
1153!    IDID        -> Reports on successfulness upon return:
1154!                    = 1 for success
1155!                    < 0 for error (value equals error code)
1156!
1157!    ISTATUS(1)  -> No. of function calls
1158!    ISTATUS(2)  -> No. of jacobian calls
1159!    ISTATUS(3)  -> No. of steps
1160!    ISTATUS(4)  -> No. of accepted steps
1161!    ISTATUS(5)  -> No. of rejected steps (except at very beginning)
1162!    ISTATUS(6)  -> No. of LU decompositions
1163!    ISTATUS(7)  -> No. of forward/backward substitutions
1164!    ISTATUS(8)  -> No. of singular matrix decompositions
1165!
1166!    RSTATUS(1)  -> Texit,the time corresponding to the
1167!                     computed Y upon return
1168!    RSTATUS(2)  -> Hexit,last accepted step before exit
1169!    RSTATUS(3)  -> Hnew,last predicted step (not yet taken)
1170!                   For multiple restarts,use Hnew as Hstart
1171!                     in the subsequent run
1172!
1173!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1174
1175
1176!~~~>  arguments
1177   INTEGER,      INTENT(IN)  :: n
1178   REAL(kind=dp), INTENT(INOUT):: y(n)
1179   REAL(kind=dp), INTENT(IN)  :: tstart, tend
1180   REAL(kind=dp), INTENT(IN)  :: abstol(n), reltol(n)
1181   INTEGER,      INTENT(IN)  :: icntrl(20)
1182   REAL(kind=dp), INTENT(IN)  :: rcntrl(20)
1183   INTEGER,      INTENT(INOUT):: istatus(20)
1184   REAL(kind=dp), INTENT(INOUT):: rstatus(20)
1185   INTEGER, INTENT(OUT) :: ierr
1186!~~~>  PARAMETERs of the rosenbrock method, up to 6 stages
1187   INTEGER ::  ros_s, rosmethod
1188   INTEGER, PARAMETER :: rs2=1, rs3=2, rs4=3, rd3=4, rd4=5, rg3=6
1189   REAL(kind=dp):: ros_a(15), ros_c(15), ros_m(6), ros_e(6), &
1190                    ros_alpha(6), ros_gamma(6), ros_elo
1191   LOGICAL :: ros_newf(6)
1192   CHARACTER(len=12):: ros_name
1193!~~~>  local variables
1194   REAL(kind=dp):: roundoff, facmin, facmax, facrej, facsafe
1195   REAL(kind=dp):: hmin, hmax, hstart
1196   REAL(kind=dp):: texit
1197   INTEGER       :: i, uplimtol, max_no_steps
1198   LOGICAL       :: autonomous, vectortol
1199!~~~>   PARAMETERs
1200   REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
1201   REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
1202
1203!~~~>  initialize statistics
1204   istatus(1:8) = 0
1205   rstatus(1:3) = zero
1206
1207!~~~>  autonomous or time dependent ode. default is time dependent.
1208   autonomous = .not.(icntrl(1) == 0)
1209
1210!~~~>  for scalar tolerances (icntrl(2).ne.0) the code uses abstol(1)and reltol(1)
1211!   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
1212   IF (icntrl(2) == 0)THEN
1213      vectortol = .TRUE.
1214      uplimtol  = n
1215   ELSE
1216      vectortol = .FALSE.
1217      uplimtol  = 1
1218   ENDIF
1219
1220!~~~>   initialize the particular rosenbrock method selected
1221   select CASE (icntrl(3))
1222     CASE (1)
1223       CALL ros2
1224     CASE (2)
1225       CALL ros3
1226     CASE (3)
1227       CALL ros4
1228     CASE (0, 4)
1229       CALL rodas3
1230     CASE (5)
1231       CALL rodas4
1232     CASE (6)
1233       CALL rang3
1234     CASE default
1235       PRINT *,'Unknown Rosenbrock method: ICNTRL(3) =',ICNTRL(3) 
1236       CALL ros_errormsg(- 2, tstart, zero, ierr)
1237       RETURN
1238   END select
1239
1240!~~~>   the maximum number of steps admitted
1241   IF (icntrl(4) == 0)THEN
1242      max_no_steps = 200000
1243   ELSEIF (icntrl(4)> 0)THEN
1244      max_no_steps=icntrl(4)
1245   ELSE
1246      PRINT *,'User-selected max no. of steps: ICNTRL(4) =',ICNTRL(4)
1247      CALL ros_errormsg(- 1, tstart, zero, ierr)
1248      RETURN
1249   ENDIF
1250
1251!~~~>  unit roundoff (1+ roundoff>1)
1252   roundoff = epsilon(one)
1253
1254!~~~>  lower bound on the step size: (positive value)
1255   IF (rcntrl(1) == zero)THEN
1256      hmin = zero
1257   ELSEIF (rcntrl(1)> zero)THEN
1258      hmin = rcntrl(1)
1259   ELSE
1260      PRINT *,'User-selected Hmin: RCNTRL(1) =',RCNTRL(1)
1261      CALL ros_errormsg(- 3, tstart, zero, ierr)
1262      RETURN
1263   ENDIF
1264!~~~>  upper bound on the step size: (positive value)
1265   IF (rcntrl(2) == zero)THEN
1266      hmax = abs(tend-tstart)
1267   ELSEIF (rcntrl(2)> zero)THEN
1268      hmax = min(abs(rcntrl(2)), abs(tend-tstart))
1269   ELSE
1270      PRINT *,'User-selected Hmax: RCNTRL(2) =',RCNTRL(2)
1271      CALL ros_errormsg(- 3, tstart, zero, ierr)
1272      RETURN
1273   ENDIF
1274!~~~>  starting step size: (positive value)
1275   IF (rcntrl(3) == zero)THEN
1276      hstart = max(hmin, deltamin)
1277   ELSEIF (rcntrl(3)> zero)THEN
1278      hstart = min(abs(rcntrl(3)), abs(tend-tstart))
1279   ELSE
1280      PRINT *,'User-selected Hstart: RCNTRL(3) =',RCNTRL(3)
1281      CALL ros_errormsg(- 3, tstart, zero, ierr)
1282      RETURN
1283   ENDIF
1284!~~~>  step size can be changed s.t.  facmin < hnew/hold < facmax
1285   IF (rcntrl(4) == zero)THEN
1286      facmin = 0.2_dp
1287   ELSEIF (rcntrl(4)> zero)THEN
1288      facmin = rcntrl(4)
1289   ELSE
1290      PRINT *,'User-selected FacMin: RCNTRL(4) =',RCNTRL(4)
1291      CALL ros_errormsg(- 4, tstart, zero, ierr)
1292      RETURN
1293   ENDIF
1294   IF (rcntrl(5) == zero)THEN
1295      facmax = 6.0_dp
1296   ELSEIF (rcntrl(5)> zero)THEN
1297      facmax = rcntrl(5)
1298   ELSE
1299      PRINT *,'User-selected FacMax: RCNTRL(5) =',RCNTRL(5)
1300      CALL ros_errormsg(- 4, tstart, zero, ierr)
1301      RETURN
1302   ENDIF
1303!~~~>   facrej: factor to decrease step after 2 succesive rejections
1304   IF (rcntrl(6) == zero)THEN
1305      facrej = 0.1_dp
1306   ELSEIF (rcntrl(6)> zero)THEN
1307      facrej = rcntrl(6)
1308   ELSE
1309      PRINT *,'User-selected FacRej: RCNTRL(6) =',RCNTRL(6)
1310      CALL ros_errormsg(- 4, tstart, zero, ierr)
1311      RETURN
1312   ENDIF
1313!~~~>   facsafe: safety factor in the computation of new step size
1314   IF (rcntrl(7) == zero)THEN
1315      facsafe = 0.9_dp
1316   ELSEIF (rcntrl(7)> zero)THEN
1317      facsafe = rcntrl(7)
1318   ELSE
1319      PRINT *,'User-selected FacSafe: RCNTRL(7) =',RCNTRL(7)
1320      CALL ros_errormsg(- 4, tstart, zero, ierr)
1321      RETURN
1322   ENDIF
1323!~~~>  check IF tolerances are reasonable
1324    DO i=1, uplimtol
1325      IF ((abstol(i)<= zero).or. (reltol(i)<= 10.0_dp* roundoff)&
1326         .or. (reltol(i)>= 1.0_dp))THEN
1327        PRINT *,' AbsTol(',i,') = ',AbsTol(i)
1328        PRINT *,' RelTol(',i,') = ',RelTol(i)
1329        CALL ros_errormsg(- 5, tstart, zero, ierr)
1330        RETURN
1331      ENDIF
1332    ENDDO
1333
1334
1335!~~~>  CALL rosenbrock method
1336   CALL ros_integrator(y, tstart, tend, texit,  &
1337        abstol, reltol,                         &
1338!  Integration parameters
1339        autonomous, vectortol, max_no_steps,    &
1340        roundoff, hmin, hmax, hstart,           &
1341        facmin, facmax, facrej, facsafe,        &
1342!  Error indicator
1343        ierr)
1344
1345!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1346CONTAINS !  SUBROUTINEs internal to rosenbrock
1347!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1348   
1349!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1350 SUBROUTINE ros_errormsg(code, t, h, ierr)
1351!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1352!    Handles all error messages
1353!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1354 
1355   REAL(kind=dp), INTENT(IN):: t, h
1356   INTEGER, INTENT(IN) :: code
1357   INTEGER, INTENT(OUT):: ierr
1358   
1359   ierr = code
1360   print * , &
1361     'Forced exit from Rosenbrock due to the following error:' 
1362     
1363   select CASE (code)
1364    CASE (- 1) 
1365      PRINT *,'--> Improper value for maximal no of steps'
1366    CASE (- 2) 
1367      PRINT *,'--> Selected Rosenbrock method not implemented'
1368    CASE (- 3) 
1369      PRINT *,'--> Hmin/Hmax/Hstart must be positive'
1370    CASE (- 4) 
1371      PRINT *,'--> FacMin/FacMax/FacRej must be positive'
1372    CASE (- 5)
1373      PRINT *,'--> Improper tolerance values'
1374    CASE (- 6)
1375      PRINT *,'--> No of steps exceeds maximum bound'
1376    CASE (- 7)
1377      PRINT *,'--> Step size too small: T + 10*H = T',&
1378            ' or H < Roundoff'
1379    CASE (- 8) 
1380      PRINT *,'--> Matrix is repeatedly singular'
1381    CASE default
1382      PRINT *,'Unknown Error code: ',Code
1383   END select
1384   
1385   print * , "t=", t, "and h=", h
1386     
1387 END SUBROUTINE ros_errormsg
1388
1389!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1390 SUBROUTINE ros_integrator (y, tstart, tend, t, &
1391        abstol, reltol,                         &
1392!~~~> integration PARAMETERs
1393        autonomous, vectortol, max_no_steps,    &
1394        roundoff, hmin, hmax, hstart,           &
1395        facmin, facmax, facrej, facsafe,        &
1396!~~~> error indicator
1397        ierr)
1398!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1399!   Template for the implementation of a generic Rosenbrock method
1400!      defined by ros_S (no of stages)
1401!      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1402!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1403
1404
1405!~~~> input: the initial condition at tstart; output: the solution at t
1406   REAL(kind=dp), INTENT(INOUT):: y(n)
1407!~~~> input: integration interval
1408   REAL(kind=dp), INTENT(IN):: tstart, tend
1409!~~~> output: time at which the solution is RETURNed (t=tendIF success)
1410   REAL(kind=dp), INTENT(OUT)::  t
1411!~~~> input: tolerances
1412   REAL(kind=dp), INTENT(IN)::  abstol(n), reltol(n)
1413!~~~> input: integration PARAMETERs
1414   LOGICAL, INTENT(IN):: autonomous, vectortol
1415   REAL(kind=dp), INTENT(IN):: hstart, hmin, hmax
1416   INTEGER, INTENT(IN):: max_no_steps
1417   REAL(kind=dp), INTENT(IN):: roundoff, facmin, facmax, facrej, facsafe
1418!~~~> output: error indicator
1419   INTEGER, INTENT(OUT):: ierr
1420! ~~~~ Local variables
1421   REAL(kind=dp):: ynew(n), fcn0(n), fcn(n)
1422   REAL(kind=dp):: k(n* ros_s), dfdt(n)
1423#ifdef full_algebra   
1424   REAL(kind=dp):: jac0(n, n), ghimj(n, n)
1425#else
1426   REAL(kind=dp):: jac0(lu_nonzero), ghimj(lu_nonzero)
1427#endif
1428   REAL(kind=dp):: h, hnew, hc, hg, fac, tau
1429   REAL(kind=dp):: err, yerr(n)
1430   INTEGER :: pivot(n), direction, ioffset, j, istage
1431   LOGICAL :: rejectlasth, rejectmoreh, singular
1432!~~~>  local PARAMETERs
1433   REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
1434   REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
1435!~~~>  locally called FUNCTIONs
1436!    REAL(kind=dp) WLAMCH
1437!    EXTERNAL WLAMCH
1438!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1439
1440
1441!~~~>  initial preparations
1442   t = tstart
1443   rstatus(nhexit) = zero
1444   h = min( max(abs(hmin), abs(hstart)), abs(hmax))
1445   IF (abs(h)<= 10.0_dp* roundoff)h = deltamin
1446
1447   IF (tend  >=  tstart)THEN
1448     direction = + 1
1449   ELSE
1450     direction = - 1
1451   ENDIF
1452   h = direction* h
1453
1454   rejectlasth=.FALSE.
1455   rejectmoreh=.FALSE.
1456
1457!~~~> time loop begins below
1458
1459timeloop: DO WHILE((direction > 0).and.((t- tend) + roundoff <= zero)&
1460       .or. (direction < 0).and.((tend-t) + roundoff <= zero))
1461
1462   IF (istatus(nstp)> max_no_steps)THEN  ! too many steps
1463      CALL ros_errormsg(- 6, t, h, ierr)
1464      RETURN
1465   ENDIF
1466   IF (((t+ 0.1_dp* h) == t).or.(h <= roundoff))THEN  ! step size too small
1467      CALL ros_errormsg(- 7, t, h, ierr)
1468      RETURN
1469   ENDIF
1470
1471!~~~>  limit h IF necessary to avoid going beyond tend
1472   h = min(h, abs(tend-t))
1473
1474!~~~>   compute the FUNCTION at current time
1475   CALL funtemplate(t, y, fcn0)
1476   istatus(nfun) = istatus(nfun) + 1
1477
1478!~~~>  compute the FUNCTION derivative with respect to t
1479   IF (.not.autonomous)THEN
1480      CALL ros_funtimederivative(t, roundoff, y, &
1481                fcn0, dfdt)
1482   ENDIF
1483
1484!~~~>   compute the jacobian at current time
1485   CALL jactemplate(t, y, jac0)
1486   istatus(njac) = istatus(njac) + 1
1487
1488!~~~>  repeat step calculation until current step accepted
1489untilaccepted: do
1490
1491   CALL ros_preparematrix(h, direction, ros_gamma(1), &
1492          jac0, ghimj, pivot, singular)
1493   IF (singular)THEN ! more than 5 consecutive failed decompositions
1494       CALL ros_errormsg(- 8, t, h, ierr)
1495       RETURN
1496   ENDIF
1497
1498!~~~>   compute the stages
1499stage: DO istage = 1, ros_s
1500
1501      ! current istage offset. current istage vector is k(ioffset+ 1:ioffset+ n)
1502       ioffset = n* (istage-1)
1503
1504      ! for the 1st istage the FUNCTION has been computed previously
1505       IF (istage == 1)THEN
1506         !slim: CALL wcopy(n, fcn0, 1, fcn, 1)
1507       fcn(1:n) = fcn0(1:n)
1508      ! istage>1 and a new FUNCTION evaluation is needed at the current istage
1509       ELSEIF(ros_newf(istage))THEN
1510         !slim: CALL wcopy(n, y, 1, ynew, 1)
1511       ynew(1:n) = y(1:n)
1512         DO j = 1, istage-1
1513           CALL waxpy(n, ros_a((istage-1) * (istage-2) /2+ j), &
1514            k(n* (j- 1) + 1), 1, ynew, 1)
1515         ENDDO
1516         tau = t + ros_alpha(istage) * direction* h
1517         CALL funtemplate(tau, ynew, fcn)
1518         istatus(nfun) = istatus(nfun) + 1
1519       ENDIF ! IF istage == 1 ELSEIF ros_newf(istage)
1520       !slim: CALL wcopy(n, fcn, 1, k(ioffset+ 1), 1)
1521       k(ioffset+ 1:ioffset+ n) = fcn(1:n)
1522       DO j = 1, istage-1
1523         hc = ros_c((istage-1) * (istage-2) /2+ j) /(direction* h)
1524         CALL waxpy(n, hc, k(n* (j- 1) + 1), 1, k(ioffset+ 1), 1)
1525       ENDDO
1526       IF ((.not. autonomous).and.(ros_gamma(istage).ne.zero))THEN
1527         hg = direction* h* ros_gamma(istage)
1528         CALL waxpy(n, hg, dfdt, 1, k(ioffset+ 1), 1)
1529       ENDIF
1530       CALL ros_solve(ghimj, pivot, k(ioffset+ 1))
1531
1532   END DO stage
1533
1534
1535!~~~>  compute the new solution
1536   !slim: CALL wcopy(n, y, 1, ynew, 1)
1537   ynew(1:n) = y(1:n)
1538   DO j=1, ros_s
1539         CALL waxpy(n, ros_m(j), k(n* (j- 1) + 1), 1, ynew, 1)
1540   ENDDO
1541
1542!~~~>  compute the error estimation
1543   !slim: CALL wscal(n, zero, yerr, 1)
1544   yerr(1:n) = zero
1545   DO j=1, ros_s
1546        CALL waxpy(n, ros_e(j), k(n* (j- 1) + 1), 1, yerr, 1)
1547   ENDDO
1548   err = ros_errornorm(y, ynew, yerr, abstol, reltol, vectortol)
1549
1550!~~~> new step size is bounded by facmin <= hnew/h <= facmax
1551   fac  = min(facmax, max(facmin, facsafe/err** (one/ros_elo)))
1552   hnew = h* fac
1553
1554!~~~>  check the error magnitude and adjust step size
1555   istatus(nstp) = istatus(nstp) + 1
1556   IF ((err <= one).or.(h <= hmin))THEN  !~~~> accept step
1557      istatus(nacc) = istatus(nacc) + 1
1558      !slim: CALL wcopy(n, ynew, 1, y, 1)
1559      y(1:n) = ynew(1:n)
1560      t = t + direction* h
1561      hnew = max(hmin, min(hnew, hmax))
1562      IF (rejectlasth)THEN  ! no step size increase after a rejected step
1563         hnew = min(hnew, h)
1564      ENDIF
1565      rstatus(nhexit) = h
1566      rstatus(nhnew) = hnew
1567      rstatus(ntexit) = t
1568      rejectlasth = .FALSE.
1569      rejectmoreh = .FALSE.
1570      h = hnew
1571      exit untilaccepted ! exit the loop: WHILE step not accepted
1572   ELSE           !~~~> reject step
1573      IF (rejectmoreh)THEN
1574         hnew = h* facrej
1575      ENDIF
1576      rejectmoreh = rejectlasth
1577      rejectlasth = .TRUE.
1578      h = hnew
1579      IF (istatus(nacc)>= 1) istatus(nrej) = istatus(nrej) + 1
1580   ENDIF ! err <= 1
1581
1582   END DO untilaccepted
1583
1584   END DO timeloop
1585
1586!~~~> succesful exit
1587   ierr = 1  !~~~> the integration was successful
1588
1589  END SUBROUTINE ros_integrator
1590
1591
1592!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1593  REAL(kind=dp)FUNCTION ros_errornorm(y, ynew, yerr, &
1594                               abstol, reltol, vectortol)
1595!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1596!~~~> computes the "scaled norm" of the error vector yerr
1597!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1598
1599! Input arguments
1600   REAL(kind=dp), INTENT(IN):: y(n), ynew(n), &
1601          yerr(n), abstol(n), reltol(n)
1602   LOGICAL, INTENT(IN)::  vectortol
1603! Local variables
1604   REAL(kind=dp):: err, scale, ymax
1605   INTEGER  :: i
1606   REAL(kind=dp), PARAMETER :: zero = 0.0_dp
1607
1608   err = zero
1609   DO i=1, n
1610     ymax = max(abs(y(i)), abs(ynew(i)))
1611     IF (vectortol)THEN
1612       scale = abstol(i) + reltol(i) * ymax
1613     ELSE
1614       scale = abstol(1) + reltol(1) * ymax
1615     ENDIF
1616     err = err+ (yerr(i) /scale) ** 2
1617   ENDDO
1618   err  = sqrt(err/n)
1619
1620   ros_errornorm = max(err, 1.0d-10)
1621
1622  END FUNCTION ros_errornorm
1623
1624
1625!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1626  SUBROUTINE ros_funtimederivative(t, roundoff, y, &
1627                fcn0, dfdt)
1628!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1629!~~~> the time partial derivative of the FUNCTION by finite differences
1630!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1631
1632!~~~> input arguments
1633   REAL(kind=dp), INTENT(IN):: t, roundoff, y(n), fcn0(n)
1634!~~~> output arguments
1635   REAL(kind=dp), INTENT(OUT):: dfdt(n)
1636!~~~> local variables
1637   REAL(kind=dp):: delta
1638   REAL(kind=dp), PARAMETER :: one = 1.0_dp, deltamin = 1.0e-6_dp
1639
1640   delta = sqrt(roundoff) * max(deltamin, abs(t))
1641   CALL funtemplate(t+ delta, y, dfdt)
1642   istatus(nfun) = istatus(nfun) + 1
1643   CALL waxpy(n, (- one), fcn0, 1, dfdt, 1)
1644   CALL wscal(n, (one/delta), dfdt, 1)
1645
1646  END SUBROUTINE ros_funtimederivative
1647
1648
1649!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1650  SUBROUTINE ros_preparematrix(h, direction, gam, &
1651             jac0, ghimj, pivot, singular)
1652! --- --- --- --- --- --- --- --- --- --- --- --- ---
1653!  Prepares the LHS matrix for stage calculations
1654!  1.  Construct Ghimj = 1/(H*ham) - Jac0
1655!      "(Gamma H) Inverse Minus Jacobian"
1656!  2.  Repeat LU decomposition of Ghimj until successful.
1657!       -half the step size if LU decomposition fails and retry
1658!       -exit after 5 consecutive fails
1659! --- --- --- --- --- --- --- --- --- --- --- --- ---
1660
1661!~~~> input arguments
1662#ifdef full_algebra   
1663   REAL(kind=dp), INTENT(IN)::  jac0(n, n)
1664#else
1665   REAL(kind=dp), INTENT(IN)::  jac0(lu_nonzero)
1666#endif   
1667   REAL(kind=dp), INTENT(IN)::  gam
1668   INTEGER, INTENT(IN)::  direction
1669!~~~> output arguments
1670#ifdef full_algebra   
1671   REAL(kind=dp), INTENT(OUT):: ghimj(n, n)
1672#else
1673   REAL(kind=dp), INTENT(OUT):: ghimj(lu_nonzero)
1674#endif   
1675   LOGICAL, INTENT(OUT)::  singular
1676   INTEGER, INTENT(OUT)::  pivot(n)
1677!~~~> inout arguments
1678   REAL(kind=dp), INTENT(INOUT):: h   ! step size is decreased when lu fails
1679!~~~> local variables
1680   INTEGER  :: i, ising, nconsecutive
1681   REAL(kind=dp):: ghinv
1682   REAL(kind=dp), PARAMETER :: one  = 1.0_dp, half = 0.5_dp
1683
1684   nconsecutive = 0
1685   singular = .TRUE.
1686
1687   DO WHILE (singular)
1688
1689!~~~>    construct ghimj = 1/(h* gam) - jac0
1690#ifdef full_algebra   
1691     !slim: CALL wcopy(n* n, jac0, 1, ghimj, 1)
1692     !slim: CALL wscal(n* n, (- one), ghimj, 1)
1693     ghimj = - jac0
1694     ghinv = one/(direction* h* gam)
1695     DO i=1, n
1696       ghimj(i, i) = ghimj(i, i) + ghinv
1697     ENDDO
1698#else
1699     !slim: CALL wcopy(lu_nonzero, jac0, 1, ghimj, 1)
1700     !slim: CALL wscal(lu_nonzero, (- one), ghimj, 1)
1701     ghimj(1:lu_nonzero) = - jac0(1:lu_nonzero)
1702     ghinv = one/(direction* h* gam)
1703     DO i=1, n
1704       ghimj(lu_diag(i)) = ghimj(lu_diag(i)) + ghinv
1705     ENDDO
1706#endif   
1707!~~~>    compute lu decomposition
1708     CALL ros_decomp( ghimj, pivot, ising)
1709     IF (ising == 0)THEN
1710!~~~>    IF successful done
1711        singular = .FALSE.
1712     ELSE ! ising .ne. 0
1713!~~~>    IF unsuccessful half the step size; IF 5 consecutive fails THEN RETURN
1714        istatus(nsng) = istatus(nsng) + 1
1715        nconsecutive = nconsecutive+1
1716        singular = .TRUE.
1717        PRINT*,'Warning: LU Decomposition returned ISING = ',ISING
1718        IF (nconsecutive <= 5)THEN ! less than 5 consecutive failed decompositions
1719           h = h* half
1720        ELSE  ! more than 5 consecutive failed decompositions
1721           RETURN
1722        ENDIF  ! nconsecutive
1723      ENDIF    ! ising
1724
1725   END DO ! WHILE singular
1726
1727  END SUBROUTINE ros_preparematrix
1728
1729
1730!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1731  SUBROUTINE ros_decomp( a, pivot, ising)
1732!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1733!  Template for the LU decomposition
1734!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1735!~~~> inout variables
1736#ifdef full_algebra   
1737   REAL(kind=dp), INTENT(INOUT):: a(n, n)
1738#else   
1739   REAL(kind=dp), INTENT(INOUT):: a(lu_nonzero)
1740#endif
1741!~~~> output variables
1742   INTEGER, INTENT(OUT):: pivot(n), ising
1743
1744#ifdef full_algebra   
1745   CALL  dgetrf( n, n, a, n, pivot, ising)
1746#else   
1747   CALL kppdecomp(a, ising)
1748   pivot(1) = 1
1749#endif
1750   istatus(ndec) = istatus(ndec) + 1
1751
1752  END SUBROUTINE ros_decomp
1753
1754
1755!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1756  SUBROUTINE ros_solve( a, pivot, b)
1757!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1758!  Template for the forward/backward substitution (using pre-computed LU decomposition)
1759!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1760!~~~> input variables
1761#ifdef full_algebra   
1762   REAL(kind=dp), INTENT(IN):: a(n, n)
1763   INTEGER :: ising
1764#else   
1765   REAL(kind=dp), INTENT(IN):: a(lu_nonzero)
1766#endif
1767   INTEGER, INTENT(IN):: pivot(n)
1768!~~~> inout variables
1769   REAL(kind=dp), INTENT(INOUT):: b(n)
1770
1771! Following line is just to avoid compiler message about unused variables
1772   IF ( pivot(1) == 0 )  CONTINUE
1773
1774#ifdef full_algebra   
1775   CALL  DGETRS( 'N',N ,1,A,N,Pivot,b,N,ISING)
1776   IF (info < 0)THEN
1777      print* , "error in dgetrs. ising=", ising
1778   ENDIF 
1779#else   
1780   CALL kppsolve( a, b)
1781#endif
1782
1783   istatus(nsol) = istatus(nsol) + 1
1784
1785  END SUBROUTINE ros_solve
1786
1787
1788
1789!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1790  SUBROUTINE ros2
1791!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1792! --- AN L-STABLE METHOD,2 stages,order 2
1793!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1794
1795   double precision g
1796
1797    g = 1.0_dp + 1.0_dp/sqrt(2.0_dp)
1798    rosmethod = rs2
1799!~~~> name of the method
1800    ros_Name = 'ROS-2'
1801!~~~> number of stages
1802    ros_s = 2
1803
1804!~~~> the coefficient matrices a and c are strictly lower triangular.
1805!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1806!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1807!   The general mapping formula is:
1808!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1809!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1810
1811    ros_a(1) = (1.0_dp) /g
1812    ros_c(1) = (- 2.0_dp) /g
1813!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1814!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1815    ros_newf(1) = .TRUE.
1816    ros_newf(2) = .TRUE.
1817!~~~> m_i = coefficients for new step solution
1818    ros_m(1) = (3.0_dp) /(2.0_dp* g)
1819    ros_m(2) = (1.0_dp) /(2.0_dp* g)
1820! E_i = Coefficients for error estimator
1821    ros_e(1) = 1.0_dp/(2.0_dp* g)
1822    ros_e(2) = 1.0_dp/(2.0_dp* g)
1823!~~~> ros_elo = estimator of local order - the minimum between the
1824!    main and the embedded scheme orders plus one
1825    ros_elo = 2.0_dp
1826!~~~> y_stage_i ~ y( t + h* alpha_i)
1827    ros_alpha(1) = 0.0_dp
1828    ros_alpha(2) = 1.0_dp
1829!~~~> gamma_i = \sum_j  gamma_{i, j}
1830    ros_gamma(1) = g
1831    ros_gamma(2) = -g
1832
1833 END SUBROUTINE ros2
1834
1835
1836!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1837  SUBROUTINE ros3
1838!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1839! --- AN L-STABLE METHOD,3 stages,order 3,2 function evaluations
1840!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1841
1842   rosmethod = rs3
1843!~~~> name of the method
1844   ros_Name = 'ROS-3'
1845!~~~> number of stages
1846   ros_s = 3
1847
1848!~~~> the coefficient matrices a and c are strictly lower triangular.
1849!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1850!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1851!   The general mapping formula is:
1852!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1853!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1854
1855   ros_a(1) = 1.0_dp
1856   ros_a(2) = 1.0_dp
1857   ros_a(3) = 0.0_dp
1858
1859   ros_c(1) = - 0.10156171083877702091975600115545e+01_dp
1860   ros_c(2) =  0.40759956452537699824805835358067e+01_dp
1861   ros_c(3) =  0.92076794298330791242156818474003e+01_dp
1862!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1863!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1864   ros_newf(1) = .TRUE.
1865   ros_newf(2) = .TRUE.
1866   ros_newf(3) = .FALSE.
1867!~~~> m_i = coefficients for new step solution
1868   ros_m(1) =  0.1e+01_dp
1869   ros_m(2) =  0.61697947043828245592553615689730e+01_dp
1870   ros_m(3) = - 0.42772256543218573326238373806514_dp
1871! E_i = Coefficients for error estimator
1872   ros_e(1) =  0.5_dp
1873   ros_e(2) = - 0.29079558716805469821718236208017e+01_dp
1874   ros_e(3) =  0.22354069897811569627360909276199_dp
1875!~~~> ros_elo = estimator of local order - the minimum between the
1876!    main and the embedded scheme orders plus 1
1877   ros_elo = 3.0_dp
1878!~~~> y_stage_i ~ y( t + h* alpha_i)
1879   ros_alpha(1) = 0.0_dp
1880   ros_alpha(2) = 0.43586652150845899941601945119356_dp
1881   ros_alpha(3) = 0.43586652150845899941601945119356_dp
1882!~~~> gamma_i = \sum_j  gamma_{i, j}
1883   ros_gamma(1) = 0.43586652150845899941601945119356_dp
1884   ros_gamma(2) = 0.24291996454816804366592249683314_dp
1885   ros_gamma(3) = 0.21851380027664058511513169485832e+01_dp
1886
1887  END SUBROUTINE ros3
1888
1889!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1890
1891
1892!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1893  SUBROUTINE ros4
1894!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1895!     L-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 4 STAGES
1896!     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
1897!
1898!      E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
1899!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1900!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1901!      SPRINGER-VERLAG (1990)
1902!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1903
1904
1905   rosmethod = rs4
1906!~~~> name of the method
1907   ros_Name = 'ROS-4'
1908!~~~> number of stages
1909   ros_s = 4
1910
1911!~~~> the coefficient matrices a and c are strictly lower triangular.
1912!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1913!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1914!   The general mapping formula is:
1915!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1916!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1917
1918   ros_a(1) = 0.2000000000000000e+01_dp
1919   ros_a(2) = 0.1867943637803922e+01_dp
1920   ros_a(3) = 0.2344449711399156_dp
1921   ros_a(4) = ros_a(2)
1922   ros_a(5) = ros_a(3)
1923   ros_a(6) = 0.0_dp
1924
1925   ros_c(1) = -0.7137615036412310e+01_dp
1926   ros_c(2) = 0.2580708087951457e+01_dp
1927   ros_c(3) = 0.6515950076447975_dp
1928   ros_c(4) = -0.2137148994382534e+01_dp
1929   ros_c(5) = -0.3214669691237626_dp
1930   ros_c(6) = -0.6949742501781779_dp
1931!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1932!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1933   ros_newf(1) = .TRUE.
1934   ros_newf(2) = .TRUE.
1935   ros_newf(3) = .TRUE.
1936   ros_newf(4) = .FALSE.
1937!~~~> m_i = coefficients for new step solution
1938   ros_m(1) = 0.2255570073418735e+01_dp
1939   ros_m(2) = 0.2870493262186792_dp
1940   ros_m(3) = 0.4353179431840180_dp
1941   ros_m(4) = 0.1093502252409163e+01_dp
1942!~~~> e_i  = coefficients for error estimator
1943   ros_e(1) = -0.2815431932141155_dp
1944   ros_e(2) = -0.7276199124938920e-01_dp
1945   ros_e(3) = -0.1082196201495311_dp
1946   ros_e(4) = -0.1093502252409163e+01_dp
1947!~~~> ros_elo  = estimator of local order - the minimum between the
1948!    main and the embedded scheme orders plus 1
1949   ros_elo  = 4.0_dp
1950!~~~> y_stage_i ~ y( t + h* alpha_i)
1951   ros_alpha(1) = 0.0_dp
1952   ros_alpha(2) = 0.1145640000000000e+01_dp
1953   ros_alpha(3) = 0.6552168638155900_dp
1954   ros_alpha(4) = ros_alpha(3)
1955!~~~> gamma_i = \sum_j  gamma_{i, j}
1956   ros_gamma(1) = 0.5728200000000000_dp
1957   ros_gamma(2) = -0.1769193891319233e+01_dp
1958   ros_gamma(3) = 0.7592633437920482_dp
1959   ros_gamma(4) = -0.1049021087100450_dp
1960
1961  END SUBROUTINE ros4
1962
1963!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1964  SUBROUTINE rodas3
1965!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1966! --- A STIFFLY-STABLE METHOD,4 stages,order 3
1967!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1968
1969
1970   rosmethod = rd3
1971!~~~> name of the method
1972   ros_Name = 'RODAS-3'
1973!~~~> number of stages
1974   ros_s = 4
1975
1976!~~~> the coefficient matrices a and c are strictly lower triangular.
1977!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1978!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1979!   The general mapping formula is:
1980!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1981!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1982
1983   ros_a(1) = 0.0_dp
1984   ros_a(2) = 2.0_dp
1985   ros_a(3) = 0.0_dp
1986   ros_a(4) = 2.0_dp
1987   ros_a(5) = 0.0_dp
1988   ros_a(6) = 1.0_dp
1989
1990   ros_c(1) = 4.0_dp
1991   ros_c(2) = 1.0_dp
1992   ros_c(3) = -1.0_dp
1993   ros_c(4) = 1.0_dp
1994   ros_c(5) = -1.0_dp
1995   ros_c(6) = -(8.0_dp/3.0_dp)
1996
1997!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1998!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1999   ros_newf(1) = .TRUE.
2000   ros_newf(2) = .FALSE.
2001   ros_newf(3) = .TRUE.
2002   ros_newf(4) = .TRUE.
2003!~~~> m_i = coefficients for new step solution
2004   ros_m(1) = 2.0_dp
2005   ros_m(2) = 0.0_dp
2006   ros_m(3) = 1.0_dp
2007   ros_m(4) = 1.0_dp
2008!~~~> e_i  = coefficients for error estimator
2009   ros_e(1) = 0.0_dp
2010   ros_e(2) = 0.0_dp
2011   ros_e(3) = 0.0_dp
2012   ros_e(4) = 1.0_dp
2013!~~~> ros_elo  = estimator of local order - the minimum between the
2014!    main and the embedded scheme orders plus 1
2015   ros_elo  = 3.0_dp
2016!~~~> y_stage_i ~ y( t + h* alpha_i)
2017   ros_alpha(1) = 0.0_dp
2018   ros_alpha(2) = 0.0_dp
2019   ros_alpha(3) = 1.0_dp
2020   ros_alpha(4) = 1.0_dp
2021!~~~> gamma_i = \sum_j  gamma_{i, j}
2022   ros_gamma(1) = 0.5_dp
2023   ros_gamma(2) = 1.5_dp
2024   ros_gamma(3) = 0.0_dp
2025   ros_gamma(4) = 0.0_dp
2026
2027  END SUBROUTINE rodas3
2028
2029!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2030  SUBROUTINE rodas4
2031!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2032!     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 6 STAGES
2033!
2034!      E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
2035!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
2036!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
2037!      SPRINGER-VERLAG (1996)
2038!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2039
2040
2041    rosmethod = rd4
2042!~~~> name of the method
2043    ros_Name = 'RODAS-4'
2044!~~~> number of stages
2045    ros_s = 6
2046
2047!~~~> y_stage_i ~ y( t + h* alpha_i)
2048    ros_alpha(1) = 0.000_dp
2049    ros_alpha(2) = 0.386_dp
2050    ros_alpha(3) = 0.210_dp
2051    ros_alpha(4) = 0.630_dp
2052    ros_alpha(5) = 1.000_dp
2053    ros_alpha(6) = 1.000_dp
2054
2055!~~~> gamma_i = \sum_j  gamma_{i, j}
2056    ros_gamma(1) = 0.2500000000000000_dp
2057    ros_gamma(2) = -0.1043000000000000_dp
2058    ros_gamma(3) = 0.1035000000000000_dp
2059    ros_gamma(4) = -0.3620000000000023e-01_dp
2060    ros_gamma(5) = 0.0_dp
2061    ros_gamma(6) = 0.0_dp
2062
2063!~~~> the coefficient matrices a and c are strictly lower triangular.
2064!   The lower triangular (subdiagonal) elements are stored in row-wise order:
2065!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
2066!   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
2067!                  C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
2068
2069    ros_a(1) = 0.1544000000000000e+01_dp
2070    ros_a(2) = 0.9466785280815826_dp
2071    ros_a(3) = 0.2557011698983284_dp
2072    ros_a(4) = 0.3314825187068521e+01_dp
2073    ros_a(5) = 0.2896124015972201e+01_dp
2074    ros_a(6) = 0.9986419139977817_dp
2075    ros_a(7) = 0.1221224509226641e+01_dp
2076    ros_a(8) = 0.6019134481288629e+01_dp
2077    ros_a(9) = 0.1253708332932087e+02_dp
2078    ros_a(10) = -0.6878860361058950_dp
2079    ros_a(11) = ros_a(7)
2080    ros_a(12) = ros_a(8)
2081    ros_a(13) = ros_a(9)
2082    ros_a(14) = ros_a(10)
2083    ros_a(15) = 1.0_dp
2084
2085    ros_c(1) = -0.5668800000000000e+01_dp
2086    ros_c(2) = -0.2430093356833875e+01_dp
2087    ros_c(3) = -0.2063599157091915_dp
2088    ros_c(4) = -0.1073529058151375_dp
2089    ros_c(5) = -0.9594562251023355e+01_dp
2090    ros_c(6) = -0.2047028614809616e+02_dp
2091    ros_c(7) = 0.7496443313967647e+01_dp
2092    ros_c(8) = -0.1024680431464352e+02_dp
2093    ros_c(9) = -0.3399990352819905e+02_dp
2094    ros_c(10) = 0.1170890893206160e+02_dp
2095    ros_c(11) = 0.8083246795921522e+01_dp
2096    ros_c(12) = -0.7981132988064893e+01_dp
2097    ros_c(13) = -0.3152159432874371e+02_dp
2098    ros_c(14) = 0.1631930543123136e+02_dp
2099    ros_c(15) = -0.6058818238834054e+01_dp
2100
2101!~~~> m_i = coefficients for new step solution
2102    ros_m(1) = ros_a(7)
2103    ros_m(2) = ros_a(8)
2104    ros_m(3) = ros_a(9)
2105    ros_m(4) = ros_a(10)
2106    ros_m(5) = 1.0_dp
2107    ros_m(6) = 1.0_dp
2108
2109!~~~> e_i  = coefficients for error estimator
2110    ros_e(1) = 0.0_dp
2111    ros_e(2) = 0.0_dp
2112    ros_e(3) = 0.0_dp
2113    ros_e(4) = 0.0_dp
2114    ros_e(5) = 0.0_dp
2115    ros_e(6) = 1.0_dp
2116
2117!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2118!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2119    ros_newf(1) = .TRUE.
2120    ros_newf(2) = .TRUE.
2121    ros_newf(3) = .TRUE.
2122    ros_newf(4) = .TRUE.
2123    ros_newf(5) = .TRUE.
2124    ros_newf(6) = .TRUE.
2125
2126!~~~> ros_elo  = estimator of local order - the minimum between the
2127!        main and the embedded scheme orders plus 1
2128    ros_elo = 4.0_dp
2129
2130  END SUBROUTINE rodas4
2131!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2132  SUBROUTINE rang3
2133!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2134! STIFFLY-STABLE W METHOD OF ORDER 3,WITH 4 STAGES
2135!
2136! J. RANG and L. ANGERMANN
2137! NEW ROSENBROCK W-METHODS OF ORDER 3
2138! FOR PARTIAL DIFFERENTIAL ALGEBRAIC
2139!        EQUATIONS OF INDEX 1                                             
2140! BIT Numerical Mathematics (2005) 45: 761-787
2141!  DOI: 10.1007/s10543-005-0035-y
2142! Table 4.1-4.2
2143!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2144
2145
2146    rosmethod = rg3
2147!~~~> name of the method
2148    ros_Name = 'RANG-3'
2149!~~~> number of stages
2150    ros_s = 4
2151
2152    ros_a(1) = 5.09052051067020d+00;
2153    ros_a(2) = 5.09052051067020d+00;
2154    ros_a(3) = 0.0d0;
2155    ros_a(4) = 4.97628111010787d+00;
2156    ros_a(5) = 2.77268164715849d-02;
2157    ros_a(6) = 2.29428036027904d-01;
2158
2159    ros_c(1) = - 1.16790812312283d+01;
2160    ros_c(2) = - 1.64057326467367d+01;
2161    ros_c(3) = - 2.77268164715850d-01;
2162    ros_c(4) = - 8.38103960500476d+00;
2163    ros_c(5) = - 8.48328409199343d-01;
2164    ros_c(6) =  2.87009860433106d-01;
2165
2166    ros_m(1) =  5.22582761233094d+00;
2167    ros_m(2) = - 5.56971148154165d-01;
2168    ros_m(3) =  3.57979469353645d-01;
2169    ros_m(4) =  1.72337398521064d+00;
2170
2171    ros_e(1) = - 5.16845212784040d+00;
2172    ros_e(2) = - 1.26351942603842d+00;
2173    ros_e(3) = - 1.11022302462516d-16;
2174    ros_e(4) =  2.22044604925031d-16;
2175
2176    ros_alpha(1) = 0.0d00;
2177    ros_alpha(2) = 2.21878746765329d+00;
2178    ros_alpha(3) = 2.21878746765329d+00;
2179    ros_alpha(4) = 1.55392337535788d+00;
2180
2181    ros_gamma(1) =  4.35866521508459d-01;
2182    ros_gamma(2) = - 1.78292094614483d+00;
2183    ros_gamma(3) = - 2.46541900496934d+00;
2184    ros_gamma(4) = - 8.05529997906370d-01;
2185
2186
2187!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2188!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2189    ros_newf(1) = .TRUE.
2190    ros_newf(2) = .TRUE.
2191    ros_newf(3) = .TRUE.
2192    ros_newf(4) = .TRUE.
2193
2194!~~~> ros_elo  = estimator of local order - the minimum between the
2195!        main and the embedded scheme orders plus 1
2196    ros_elo = 3.0_dp
2197
2198  END SUBROUTINE rang3
2199!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2200
2201!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2202!   End of the set of internal Rosenbrock subroutines
2203!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2204END SUBROUTINE rosenbrock
2205 
2206SUBROUTINE funtemplate( t, y, ydot)
2207!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2208!  Template for the ODE function call.
2209!  Updates the rate coefficients (and possibly the fixed species) at each call
2210!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2211!~~~> input variables
2212   REAL(kind=dp):: t, y(nvar)
2213!~~~> output variables
2214   REAL(kind=dp):: ydot(nvar)
2215!~~~> local variables
2216   REAL(kind=dp):: told
2217
2218   told = time
2219   time = t
2220   CALL fun( y, fix, rconst, ydot)
2221   time = told
2222
2223END SUBROUTINE funtemplate
2224 
2225SUBROUTINE jactemplate( t, y, jcb)
2226!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2227!  Template for the ODE Jacobian call.
2228!  Updates the rate coefficients (and possibly the fixed species) at each call
2229!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2230!~~~> input variables
2231    REAL(kind=dp):: t, y(nvar)
2232!~~~> output variables
2233#ifdef full_algebra   
2234    REAL(kind=dp):: jv(lu_nonzero), jcb(nvar, nvar)
2235#else
2236    REAL(kind=dp):: jcb(lu_nonzero)
2237#endif   
2238!~~~> local variables
2239    REAL(kind=dp):: told
2240#ifdef full_algebra   
2241    INTEGER :: i, j
2242#endif   
2243
2244    told = time
2245    time = t
2246#ifdef full_algebra   
2247    CALL jac_sp(y, fix, rconst, jv)
2248    DO j=1, nvar
2249      DO i=1, nvar
2250         jcb(i, j) = 0.0_dp
2251      ENDDO
2252    ENDDO
2253    DO i=1, lu_nonzero
2254       jcb(lu_irow(i), lu_icol(i)) = jv(i)
2255    ENDDO
2256#else
2257    CALL jac_sp( y, fix, rconst, jcb)
2258#endif   
2259    time = told
2260
2261END SUBROUTINE jactemplate
2262 
2263  SUBROUTINE kppdecomp( jvs, ier)                                   
2264  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2265  !        sparse lu factorization                                   
2266  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2267  ! loop expansion generated by kp4                                   
2268                                                                     
2269    INTEGER  :: ier                                                   
2270    REAL(kind=dp):: jvs(lu_nonzero), a                         
2271                                                                     
2272    a = 0.                                                           
2273    ier = 0                                                           
2274                                                                     
2275!   i = 1
2276!   i = 2
2277    jvs(4) = (jvs(4)) / jvs(1) 
2278    jvs(5) = jvs(5) - jvs(2) * jvs(4)
2279    jvs(6) = jvs(6) - jvs(3) * jvs(4)
2280!   i = 3
2281    jvs(7) = (jvs(7)) / jvs(1) 
2282    a = 0.0; a = a  - jvs(2) * jvs(7)
2283    jvs(8) = (jvs(8) + a) / jvs(5) 
2284    jvs(9) = jvs(9) - jvs(3) * jvs(7) - jvs(6) * jvs(8)
2285    RETURN                                                           
2286                                                                     
2287  END SUBROUTINE kppdecomp                                           
2288 
2289SUBROUTINE get_mechanism_name                                       
2290                                                                   
2291  IMPLICIT NONE                                                     
2292
2293! Set cs_mech for check with mechanism name from namelist
2294    cs_mech = 'phstat'
2295                                                                   
2296  RETURN                                                           
2297END SUBROUTINE get_mechanism_name                                   
2298                                                                   
2299 
2300SUBROUTINE chem_gasphase_integrate (time_step_len, conc, tempi, qvapi, fakti, photo, ierrf, xnacc, xnrej, istatus, l_debug, pe, &
2301                     icntrl_i, rcntrl_i) 
2302                                                                   
2303  IMPLICIT NONE                                                     
2304                                                                   
2305  REAL(dp), INTENT(IN)                  :: time_step_len           
2306  REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: conc                   
2307  REAL(dp), DIMENSION(:, :), INTENT(IN)   :: photo                   
2308  REAL(dp), DIMENSION(:), INTENT(IN)     :: tempi                   
2309  REAL(dp), DIMENSION(:), INTENT(IN)     :: qvapi                   
2310  REAL(dp), DIMENSION(:), INTENT(IN)     :: fakti                   
2311  INTEGER, INTENT(OUT), OPTIONAL        :: ierrf(:)               
2312  INTEGER, INTENT(OUT), OPTIONAL        :: xnacc(:)               
2313  INTEGER, INTENT(OUT), OPTIONAL        :: xnrej(:)               
2314  INTEGER, INTENT(INOUT), OPTIONAL      :: istatus(:)             
2315  INTEGER, INTENT(IN), OPTIONAL         :: pe                     
2316  LOGICAL, INTENT(IN), OPTIONAL         :: l_debug                 
2317  INTEGER, DIMENSION(nkppctrl), INTENT(IN), OPTIONAL  :: icntrl_i         
2318  REAL(dp), DIMENSION(nkppctrl), INTENT(IN), OPTIONAL  :: rcntrl_i         
2319                                                                   
2320  INTEGER                                 :: k   ! loop variable     
2321  REAL(dp)                               :: dt                     
2322  INTEGER, DIMENSION(20)                :: istatus_u               
2323  INTEGER                                :: ierr_u                 
2324  INTEGER                                :: vl_dim_lo               
2325                                                                   
2326                                                                   
2327  IF (PRESENT (istatus)) istatus = 0                             
2328  IF (PRESENT (icntrl_i)) icntrl  = icntrl_i                     
2329  IF (PRESENT (rcntrl_i)) rcntrl  = rcntrl_i                     
2330                                                                   
2331  var => c(1:nvar)                                                 
2332                                                                   
2333  vl_glo = size(tempi, 1)                                           
2334                                                                   
2335  vl_dim_lo = vl_dim                                               
2336  DO k=1, vl_glo, vl_dim_lo                                           
2337    is = k                                                         
2338    ie = min(k+ vl_dim_lo-1, vl_glo)                                 
2339    vl = ie-is+ 1                                                   
2340                                                                   
2341    c(:) = conc(is, :)                                           
2342                                                                   
2343    temp = tempi(is)                                             
2344                                                                   
2345    qvap = qvapi(is)                                             
2346                                                                   
2347    fakt = fakti(is)                                             
2348                                                                   
2349    CALL initialize                                                 
2350                                                                   
2351    phot(:) = photo(is, :)                                           
2352                                                                   
2353    CALL update_rconst                                             
2354                                                                   
2355    dt = time_step_len                                             
2356                                                                   
2357    ! integrate from t=0 to t=dt                                   
2358    CALL integrate(0._dp, dt, icntrl, rcntrl, istatus_u = istatus_u, ierr_u=ierr_u)
2359                                                                   
2360                                                                   
2361   IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN                       
2362      IF (l_debug) CALL error_output(conc(is, :), ierr_u, pe)           
2363   ENDIF                                                             
2364                                                                     
2365    conc(is, :) = c(:)                                               
2366                                                                   
2367    ! RETURN diagnostic information                                 
2368                                                                   
2369    IF (PRESENT(ierrf))   ierrf(is) = ierr_u                     
2370    IF (PRESENT(xnacc))   xnacc(is) = istatus_u(4)               
2371    IF (PRESENT(xnrej))   xnrej(is) = istatus_u(5)               
2372                                                                   
2373    IF (PRESENT (istatus)) THEN                                   
2374      istatus(1:8) = istatus(1:8) + istatus_u(1:8)               
2375    ENDIF                                                         
2376                                                                   
2377  END DO                                                           
2378 
2379                                                                   
2380! Deallocate input arrays                                           
2381                                                                   
2382                                                                   
2383  data_loaded = .FALSE.                                             
2384                                                                   
2385  RETURN                                                           
2386END SUBROUTINE chem_gasphase_integrate                             
2387
2388END MODULE chem_gasphase_mod
2389 
Note: See TracBrowser for help on using the repository browser.