source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsa+phstat/chem_gasphase_mod.f90 @ 3833

Last change on this file since 3833 was 3833, checked in by forkel, 2 years ago

removed USE chem_gasphase_mod from chem_modules, apply USE chem_gasphase for nvar, nspec, cs_mech and spc_names instead

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