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

Last change on this file since 3187 was 2766, checked in by kanani, 7 years ago

Removal of chem directive, plus minor changes

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