source: palm/trunk/SOURCE/chem_gasphase_mod.f90 @ 2718

Last change on this file since 2718 was 2718, checked in by maronga, 3 years ago

deleting of deprecated files; headers updated where needed

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