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

Last change on this file since 3780 was 3780, checked in by forkel, 3 years ago

removed read from unit 10 in chemistry_model_mod.f90, added get_mechanismname

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