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