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

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

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

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