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

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

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

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