source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_phstat/chem_gasphase_mod.f90 @ 3800

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

added leading blanks to dummy statements

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