source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_passive/chem_gasphase_mod.f90 @ 3833

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