source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_passive1/chem_gasphase_mod.f90 @ 3944

Last change on this file since 3944 was 3944, checked in by maronga, 5 years ago

fix for last commit

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