source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_passive1/chem_gasphase_mod.f90 @ 3458

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

Reintegrated fixes/changes from branch chemistry

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