source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_phstatp/chem_gasphase_mod.f90 @ 3799

Last change on this file since 3799 was 3799, checked in by forkel, 5 years ago

editing in kpp4palm: add statements for avoiding unused variables, remove $Id

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