source: palm/trunk/SOURCE/chem_gasphase_mod.f90 @ 3788

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