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

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

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

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