source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsagas/chem_gasphase_mod.f90 @ 3681

Last change on this file since 3681 was 3681, checked in by hellstea, 5 years ago

Major update of pmc_interface_mod

File size: 79.4 KB
Line 
1MODULE chem_gasphase_mod
2 
3!   Mechanism: salsagas
4!
5!------------------------------------------------------------------------------!
6!
7! ******Module chem_gasphase_mod is automatically generated by kpp4palm ******
8!
9!   *********Please do NOT change this Code,it will be ovewritten *********
10!
11!------------------------------------------------------------------------------!
12! This file was created by KPP (http://people.cs.vt.edu/asandu/Software/Kpp/)
13! and kpp4palm (created by Klaus Ketelsen). kpp4palm is an adapted version
14! of KP4 (Jöckel,P.,Kerkweg,A.,Pozzer,A.,Sander,R.,Tost,H.,Riede,
15! H.,Baumgaertner,A.,Gromov,S.,and Kern,B.,2010: Development cycle 2 of
16! the Modular Earth Submodel System (MESSy2),Geosci. Model Dev.,3,717-752,
17! https://doi.org/10.5194/gmd-3-717-2010). KP4 is part of the Modular Earth
18! Submodel System (MESSy),which is is available under the  GNU General Public
19! License (GPL).
20!
21! KPP is free software; you can redistribute it and/or modify it under the terms
22! of the General Public Licence as published by the Free Software Foundation;
23! either version 2 of the License,or (at your option) any later version.
24! KPP is distributed in the hope that it will be useful,but WITHOUT ANY WARRANTY;
25! without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
26! PURPOSE. See the GNU General Public Licence for more details.
27!
28!------------------------------------------------------------------------------!
29! This file is part of the PALM model system.
30!
31! PALM is free software: you can redistribute it and/or modify it under the
32! terms of the GNU General Public License as published by the Free Software
33! Foundation,either version 3 of the License,or (at your option) any later
34! version.
35!
36! PALM is distributed in the hope that it will be useful,but WITHOUT ANY
37! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
38! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
39!
40! You should have received a copy of the GNU General Public License along with
41! PALM. If not,see <http://www.gnu.org/licenses/>.
42!
43! Copyright 1997-2019 Leibniz Universitaet Hannover
44!--------------------------------------------------------------------------------!
45!
46!
47! MODULE HEADER TEMPLATE
48!
49!  Initial version (Nov. 2016,ketelsen),for later modifications of module_header
50!  see comments in kpp4palm/src/create_kpp_module.C
51
52! Set kpp Double Precision to PALM Default Precision
53
54  USE kinds,           ONLY: dp=>wp
55
56  USE pegrid,          ONLY: myid, threads_per_task
57
58  IMPLICIT        NONE
59  PRIVATE
60  !SAVE  ! note: occurs again in automatically generated code ...
61
62!  PUBLIC :: IERR_NAMES
63 
64! PUBLIC :: SPC_NAMES,EQN_NAMES,EQN_TAGS,REQ_HET,REQ_AEROSOL,REQ_PHOTRAT &
65!         ,REQ_MCFCT,IP_MAX,jname
66
67  PUBLIC :: eqn_names, phot_names, spc_names
68  PUBLIC :: nmaxfixsteps
69  PUBLIC :: atol, rtol
70  PUBLIC :: nspec, nreact
71  PUBLIC :: temp
72  PUBLIC :: qvap
73  PUBLIC :: fakt
74  PUBLIC :: phot 
75  PUBLIC :: rconst
76  PUBLIC :: nvar
77  PUBLIC :: nphot
78  PUBLIC :: vl_dim                     ! PUBLIC to ebable other MODULEs to distiguish between scalar and vec
79 
80  PUBLIC :: initialize, integrate, update_rconst
81  PUBLIC :: chem_gasphase_integrate
82  PUBLIC :: initialize_kpp_ctrl
83
84! END OF MODULE HEADER TEMPLATE
85                                                                 
86! Variables used for vector mode                                 
87                                                                 
88  LOGICAL, PARAMETER          :: l_vector = .FALSE.           
89  INTEGER, PARAMETER          :: i_lu_di = 2
90  INTEGER, PARAMETER          :: vl_dim = 1
91  INTEGER                     :: vl                             
92                                                                 
93  INTEGER                     :: vl_glo                         
94  INTEGER                     :: is, ie                           
95                                                                 
96                                                                 
97  INTEGER, DIMENSION(vl_dim)  :: kacc, krej                       
98  INTEGER, DIMENSION(vl_dim)  :: ierrv                           
99  LOGICAL                     :: data_loaded = .FALSE.             
100! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
101!
102! Parameter Module File
103!
104! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
105!       (http://www.cs.vt.edu/~asandu/Software/KPP)
106! KPP is distributed under GPL,the general public licence
107!       (http://www.gnu.org/copyleft/gpl.html)
108! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
109! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
110!     With important contributions from:
111!        M. Damian,Villanova University,USA
112!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
113!
114! File                 : chem_gasphase_mod_Parameters.f90
115! Time                 : Fri Nov 30 13:52:19 2018
116! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
117! Equation file        : chem_gasphase_mod.kpp
118! Output root filename : chem_gasphase_mod
119!
120! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
121
122
123
124
125
126
127! NSPEC - Number of chemical species
128  INTEGER, PARAMETER :: nspec = 5 
129! NVAR - Number of Variable species
130  INTEGER, PARAMETER :: nvar = 5 
131! NVARACT - Number of Active species
132  INTEGER, PARAMETER :: nvaract = 5 
133! NFIX - Number of Fixed species
134  INTEGER, PARAMETER :: nfix = 1 
135! NREACT - Number of reactions
136  INTEGER, PARAMETER :: nreact = 5 
137! NVARST - Starting of variables in conc. vect.
138  INTEGER, PARAMETER :: nvarst = 1 
139! NFIXST - Starting of fixed in conc. vect.
140  INTEGER, PARAMETER :: nfixst = 6 
141! NONZERO - Number of nonzero entries in Jacobian
142  INTEGER, PARAMETER :: nonzero = 5 
143! LU_NONZERO - Number of nonzero entries in LU factoriz. of Jacobian
144  INTEGER, PARAMETER :: lu_nonzero = 5 
145! CNVAR - (NVAR+1) Number of elements in compressed row format
146  INTEGER, PARAMETER :: cnvar = 6 
147! CNEQN - (NREACT+1) Number stoicm elements in compressed col format
148  INTEGER, PARAMETER :: cneqn = 6 
149! NHESS - Length of Sparse Hessian
150  INTEGER, PARAMETER :: nhess = 1 
151! NMASS - Number of atoms to check mass balance
152  INTEGER, PARAMETER :: nmass = 1 
153
154! Index declaration for variable species in C and VAR
155!   VAR(ind_spc) = C(ind_spc)
156
157  INTEGER, PARAMETER, PUBLIC :: ind_hno3 = 1 
158  INTEGER, PARAMETER, PUBLIC :: ind_h2so4 = 2 
159  INTEGER, PARAMETER, PUBLIC :: ind_nh3 = 3 
160  INTEGER, PARAMETER, PUBLIC :: ind_ocnv = 4 
161  INTEGER, PARAMETER, PUBLIC :: ind_ocsv = 5 
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 = 5 
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                 : Fri Nov 30 13:52:19 2018
194! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/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):: c(nspec)
209! VAR - Concentrations of variable species (global)
210  REAL(kind=dp):: var(nvar)
211! FIX - Concentrations of fixed species (global)
212  REAL(kind=dp):: fix(nfix)
213! VAR,FIX are chunks of array C
214      EQUIVALENCE( c(1), var(1))
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! TSTART - Integration start time
222  REAL(kind=dp):: tstart
223! ATOL - Absolute tolerance
224  REAL(kind=dp):: atol(nvar)
225! RTOL - Relative tolerance
226  REAL(kind=dp):: rtol(nvar)
227! STEPMIN - Lower bound for integration step
228  REAL(kind=dp):: stepmin
229! CFACTOR - Conversion factor for concentration units
230  REAL(kind=dp):: cfactor
231
232! INLINED global variable declarations
233
234  !   declaration of global variable declarations for photolysis will come from
235
236! QVAP - Water vapor
237  REAL(kind=dp):: qvap
238! FAKT - Conversion factor
239  REAL(kind=dp):: fakt
240
241
242! INLINED global variable declarations
243
244
245
246! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
247!
248! Sparse Jacobian Data Structures File
249!
250! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
251!       (http://www.cs.vt.edu/~asandu/Software/KPP)
252! KPP is distributed under GPL,the general public licence
253!       (http://www.gnu.org/copyleft/gpl.html)
254! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
255! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
256!     With important contributions from:
257!        M. Damian,Villanova University,USA
258!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
259!
260! File                 : chem_gasphase_mod_JacobianSP.f90
261! Time                 : Fri Nov 30 13:52:19 2018
262! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
263! Equation file        : chem_gasphase_mod.kpp
264! Output root filename : chem_gasphase_mod
265!
266! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
267
268
269
270
271
272
273! Sparse Jacobian Data
274
275
276  INTEGER, PARAMETER, DIMENSION(5):: lu_irow =  (/ &
277       1, 2, 3, 4, 5 /) 
278
279  INTEGER, PARAMETER, DIMENSION(5):: lu_icol =  (/ &
280       1, 2, 3, 4, 5 /) 
281
282  INTEGER, PARAMETER, DIMENSION(6):: lu_crow =  (/ &
283       1, 2, 3, 4, 5, 6 /) 
284
285  INTEGER, PARAMETER, DIMENSION(6):: lu_diag =  (/ &
286       1, 2, 3, 4, 5, 6 /) 
287
288
289
290! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
291!
292! Utility Data Module File
293!
294! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
295!       (http://www.cs.vt.edu/~asandu/Software/KPP)
296! KPP is distributed under GPL,the general public licence
297!       (http://www.gnu.org/copyleft/gpl.html)
298! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
299! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
300!     With important contributions from:
301!        M. Damian,Villanova University,USA
302!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
303!
304! File                 : chem_gasphase_mod_Monitor.f90
305! Time                 : Fri Nov 30 13:52:19 2018
306! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
307! Equation file        : chem_gasphase_mod.kpp
308! Output root filename : chem_gasphase_mod
309!
310! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
311
312
313
314
315
316  CHARACTER(len=15), PARAMETER, DIMENSION(5):: spc_names =  (/ &
317     'HNO3           ','H2SO4          ','NH3            ',&
318     'OCNV           ','OCSV           ' /)
319
320  CHARACTER(len=100), PARAMETER, DIMENSION(5):: eqn_names =  (/ &
321     'H2SO4 --> H2SO4                                                                                     ',&
322     ' HNO3 --> HNO3                                                                                      ',&
323     '  NH3 --> NH3                                                                                       ',&
324     ' OCNV --> OCNV                                                                                      ',&
325     ' OCSV --> OCSV                                                                                      ' /)
326
327! INLINED global variables
328
329  !   inline f90_data: declaration of global variables for photolysis
330  !   REAL(kind=dp):: phot(nphot)must eventually be moved to global later for
331  INTEGER, PARAMETER :: nphot = 1
332  !   phot photolysis frequencies
333  REAL(kind=dp):: phot(nphot)
334
335  INTEGER, PARAMETER, PUBLIC :: j_no2 = 1
336
337  CHARACTER(len=15), PARAMETER, DIMENSION(nphot):: phot_names =   (/ &
338     'J_NO2          '/)
339
340! End INLINED global variables
341
342
343! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
344 
345! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
346 
347! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
348 
349! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
350 
351 
352!  variable definations from  individual module headers
353 
354! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
355!
356! Initialization File
357!
358! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
359!       (http://www.cs.vt.edu/~asandu/Software/KPP)
360! KPP is distributed under GPL,the general public licence
361!       (http://www.gnu.org/copyleft/gpl.html)
362! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
363! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
364!     With important contributions from:
365!        M. Damian,Villanova University,USA
366!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
367!
368! File                 : chem_gasphase_mod_Initialize.f90
369! Time                 : Fri Nov 30 13:52:19 2018
370! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
371! Equation file        : chem_gasphase_mod.kpp
372! Output root filename : chem_gasphase_mod
373!
374! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
375
376
377
378
379
380! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
381!
382! Numerical Integrator (Time-Stepping) File
383!
384! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
385!       (http://www.cs.vt.edu/~asandu/Software/KPP)
386! KPP is distributed under GPL,the general public licence
387!       (http://www.gnu.org/copyleft/gpl.html)
388! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
389! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
390!     With important contributions from:
391!        M. Damian,Villanova University,USA
392!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
393!
394! File                 : chem_gasphase_mod_Integrator.f90
395! Time                 : Fri Nov 30 13:52:19 2018
396! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
397! Equation file        : chem_gasphase_mod.kpp
398! Output root filename : chem_gasphase_mod
399!
400! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
401
402
403
404! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
405!
406! INTEGRATE - Integrator routine
407!   Arguments :
408!      TIN       - Start Time for Integration
409!      TOUT      - End Time for Integration
410!
411! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
412
413!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
414!  Rosenbrock - Implementation of several Rosenbrock methods:             !
415!               *Ros2                                                    !
416!               *Ros3                                                    !
417!               *Ros4                                                    !
418!               *Rodas3                                                  !
419!               *Rodas4                                                  !
420!  By default the code employs the KPP sparse linear algebra routines     !
421!  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
422!                                                                         !
423!    (C)  Adrian Sandu,August 2004                                       !
424!    Virginia Polytechnic Institute and State University                  !
425!    Contact: sandu@cs.vt.edu                                             !
426!    Revised by Philipp Miehe and Adrian Sandu,May 2006                  !                               !
427!    This implementation is part of KPP - the Kinetic PreProcessor        !
428!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
429
430
431  SAVE
432 
433!~~~>  statistics on the work performed by the rosenbrock method
434  INTEGER, PARAMETER :: nfun=1, njac=2, nstp=3, nacc=4, &
435                        nrej=5, ndec=6, nsol=7, nsng=8, &
436                        ntexit=1, nhexit=2, nhnew = 3
437
438! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
439!
440! Linear Algebra Data and Routines File
441!
442! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
443!       (http://www.cs.vt.edu/~asandu/Software/KPP)
444! KPP is distributed under GPL,the general public licence
445!       (http://www.gnu.org/copyleft/gpl.html)
446! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
447! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
448!     With important contributions from:
449!        M. Damian,Villanova University,USA
450!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
451!
452! File                 : chem_gasphase_mod_LinearAlgebra.f90
453! Time                 : Fri Nov 30 13:52:19 2018
454! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
455! Equation file        : chem_gasphase_mod.kpp
456! Output root filename : chem_gasphase_mod
457!
458! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
459
460
461
462
463
464
465! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
466!
467! The ODE Jacobian of Chemical Model File
468!
469! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
470!       (http://www.cs.vt.edu/~asandu/Software/KPP)
471! KPP is distributed under GPL,the general public licence
472!       (http://www.gnu.org/copyleft/gpl.html)
473! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
474! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
475!     With important contributions from:
476!        M. Damian,Villanova University,USA
477!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
478!
479! File                 : chem_gasphase_mod_Jacobian.f90
480! Time                 : Fri Nov 30 13:52:19 2018
481! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
482! Equation file        : chem_gasphase_mod.kpp
483! Output root filename : chem_gasphase_mod
484!
485! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
486
487
488
489
490
491
492! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
493!
494! The ODE Function of Chemical Model File
495!
496! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
497!       (http://www.cs.vt.edu/~asandu/Software/KPP)
498! KPP is distributed under GPL,the general public licence
499!       (http://www.gnu.org/copyleft/gpl.html)
500! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
501! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
502!     With important contributions from:
503!        M. Damian,Villanova University,USA
504!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
505!
506! File                 : chem_gasphase_mod_Function.f90
507! Time                 : Fri Nov 30 13:52:19 2018
508! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
509! Equation file        : chem_gasphase_mod.kpp
510! Output root filename : chem_gasphase_mod
511!
512! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
513
514
515
516
517
518! A - Rate for each equation
519  REAL(kind=dp):: a(nreact)
520
521! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
522!
523! The Reaction Rates File
524!
525! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
526!       (http://www.cs.vt.edu/~asandu/Software/KPP)
527! KPP is distributed under GPL,the general public licence
528!       (http://www.gnu.org/copyleft/gpl.html)
529! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
530! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
531!     With important contributions from:
532!        M. Damian,Villanova University,USA
533!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
534!
535! File                 : chem_gasphase_mod_Rates.f90
536! Time                 : Fri Nov 30 13:52:19 2018
537! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
538! Equation file        : chem_gasphase_mod.kpp
539! Output root filename : chem_gasphase_mod
540!
541! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
542
543
544
545
546
547! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
548!
549! Auxiliary Routines File
550!
551! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
552!       (http://www.cs.vt.edu/~asandu/Software/KPP)
553! KPP is distributed under GPL,the general public licence
554!       (http://www.gnu.org/copyleft/gpl.html)
555! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
556! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
557!     With important contributions from:
558!        M. Damian,Villanova University,USA
559!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
560!
561! File                 : chem_gasphase_mod_Util.f90
562! Time                 : Fri Nov 30 13:52:19 2018
563! Working directory    : /home/forkel-r/palmstuff/work/trunk20181130/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
564! Equation file        : chem_gasphase_mod.kpp
565! Output root filename : chem_gasphase_mod
566!
567! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
568
569
570
571
572
573
574  ! header MODULE initialize_kpp_ctrl_template
575
576  ! notes:
577  ! - l_vector is automatically defined by kp4
578  ! - vl_dim is automatically defined by kp4
579  ! - i_lu_di is automatically defined by kp4
580  ! - wanted is automatically defined by xmecca
581  ! - icntrl rcntrl are automatically defined by kpp
582  ! - "USE messy_main_tools" is in MODULE_header of messy_mecca_kpp.f90
583  ! - SAVE will be automatically added by kp4
584
585  !SAVE
586
587  ! for fixed time step control
588  ! ... max. number of fixed time steps (sum must be 1)
589  INTEGER, PARAMETER         :: nmaxfixsteps = 50
590  ! ... switch for fixed time stepping
591  LOGICAL, PUBLIC            :: l_fixed_step = .FALSE.
592  INTEGER, PUBLIC            :: nfsteps = 1
593  ! ... number of kpp control PARAMETERs
594  INTEGER, PARAMETER, PUBLIC :: nkppctrl = 20
595  !
596  INTEGER, DIMENSION(nkppctrl), PUBLIC     :: icntrl = 0
597  REAL(dp), DIMENSION(nkppctrl), PUBLIC     :: rcntrl = 0.0_dp
598  REAL(dp), DIMENSION(nmaxfixsteps), PUBLIC :: t_steps = 0.0_dp
599
600  ! END header MODULE initialize_kpp_ctrl_template
601
602 
603! Interface Block
604 
605  INTERFACE            initialize
606    MODULE PROCEDURE   initialize
607  END INTERFACE        initialize
608 
609  INTERFACE            integrate
610    MODULE PROCEDURE   integrate
611  END INTERFACE        integrate
612 
613  INTERFACE            fun
614    MODULE PROCEDURE   fun
615  END INTERFACE        fun
616 
617  INTERFACE            kppsolve
618    MODULE PROCEDURE   kppsolve
619  END INTERFACE        kppsolve
620 
621  INTERFACE            jac_sp
622    MODULE PROCEDURE   jac_sp
623  END INTERFACE        jac_sp
624 
625  INTERFACE            k_arr
626    MODULE PROCEDURE   k_arr
627  END INTERFACE        k_arr
628 
629  INTERFACE            update_rconst
630    MODULE PROCEDURE   update_rconst
631  END INTERFACE        update_rconst
632 
633  INTERFACE            arr2
634    MODULE PROCEDURE   arr2
635  END INTERFACE        arr2
636 
637  INTERFACE            initialize_kpp_ctrl
638    MODULE PROCEDURE   initialize_kpp_ctrl
639  END INTERFACE        initialize_kpp_ctrl
640 
641  INTERFACE            error_output
642    MODULE PROCEDURE   error_output
643  END INTERFACE        error_output
644 
645  INTERFACE            wscal
646    MODULE PROCEDURE   wscal
647  END INTERFACE        wscal
648 
649!INTERFACE not working  INTERFACE            waxpy
650!INTERFACE not working    MODULE PROCEDURE   waxpy
651!INTERFACE not working  END INTERFACE        waxpy
652 
653  INTERFACE            rosenbrock
654    MODULE PROCEDURE   rosenbrock
655  END INTERFACE        rosenbrock
656 
657  INTERFACE            funtemplate
658    MODULE PROCEDURE   funtemplate
659  END INTERFACE        funtemplate
660 
661  INTERFACE            jactemplate
662    MODULE PROCEDURE   jactemplate
663  END INTERFACE        jactemplate
664 
665  INTERFACE            kppdecomp
666    MODULE PROCEDURE   kppdecomp
667  END INTERFACE        kppdecomp
668 
669  INTERFACE            chem_gasphase_integrate
670    MODULE PROCEDURE   chem_gasphase_integrate
671  END INTERFACE        chem_gasphase_integrate
672 
673 
674 CONTAINS
675 
676SUBROUTINE initialize()
677
678
679  INTEGER         :: j, k
680
681  INTEGER :: i
682  REAL(kind=dp):: x
683  k = is
684  cfactor = 1.000000e+00_dp
685
686  x = (0.) * cfactor
687  DO i = 1 , nvar
688  ENDDO
689
690  x = (0.) * cfactor
691  DO i = 1 , nfix
692    fix(i) = x
693  ENDDO
694
695! constant rate coefficients
696! END constant rate coefficients
697
698! INLINED initializations
699
700! End INLINED initializations
701
702     
703END SUBROUTINE initialize
704 
705SUBROUTINE integrate( tin, tout, &
706  icntrl_u, rcntrl_u, istatus_u, rstatus_u, ierr_u)
707
708
709   REAL(kind=dp), INTENT(IN):: tin  ! start time
710   REAL(kind=dp), INTENT(IN):: tout ! END time
711   ! OPTIONAL input PARAMETERs and statistics
712   INTEGER,      INTENT(IN), OPTIONAL :: icntrl_u(20)
713   REAL(kind=dp), INTENT(IN), OPTIONAL :: rcntrl_u(20)
714   INTEGER,      INTENT(OUT), OPTIONAL :: istatus_u(20)
715   REAL(kind=dp), INTENT(OUT), OPTIONAL :: rstatus_u(20)
716   INTEGER,      INTENT(OUT), OPTIONAL :: ierr_u
717
718   REAL(kind=dp):: rcntrl(20), rstatus(20)
719   INTEGER       :: icntrl(20), istatus(20), ierr
720
721   INTEGER, SAVE :: ntotal = 0
722
723   icntrl(:) = 0
724   rcntrl(:) = 0.0_dp
725   istatus(:) = 0
726   rstatus(:) = 0.0_dp
727
728    !~~~> fine-tune the integrator:
729   icntrl(1) = 0      ! 0 - non- autonomous, 1 - autonomous
730   icntrl(2) = 0      ! 0 - vector tolerances, 1 - scalars
731
732   ! IF OPTIONAL PARAMETERs are given, and IF they are >0,
733   ! THEN they overwrite default settings.
734   IF (PRESENT(icntrl_u))THEN
735     WHERE(icntrl_u(:)> 0)icntrl(:) = icntrl_u(:)
736   ENDIF
737   IF (PRESENT(rcntrl_u))THEN
738     WHERE(rcntrl_u(:)> 0)rcntrl(:) = rcntrl_u(:)
739   ENDIF
740
741
742   CALL rosenbrock(nvar, var, tin, tout,  &
743         atol, rtol,               &
744         rcntrl, icntrl, rstatus, istatus, ierr)
745
746   !~~~> debug option: show no of steps
747   ! ntotal = ntotal + istatus(nstp)
748   ! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')','  O3=',VAR(ind_O3)
749
750   stepmin = rstatus(nhexit)
751   ! IF OPTIONAL PARAMETERs are given for output they
752   ! are updated with the RETURN information
753   IF (PRESENT(istatus_u))istatus_u(:) = istatus(:)
754   IF (PRESENT(rstatus_u))rstatus_u(:) = rstatus(:)
755   IF (PRESENT(ierr_u))  ierr_u       = ierr
756
757END SUBROUTINE integrate
758 
759SUBROUTINE fun(v, f, rct, vdot)
760
761! V - Concentrations of variable species (local)
762  REAL(kind=dp):: v(nvar)
763! F - Concentrations of fixed species (local)
764  REAL(kind=dp):: f(nfix)
765! RCT - Rate constants (local)
766  REAL(kind=dp):: rct(nreact)
767! Vdot - Time derivative of variable species concentrations
768  REAL(kind=dp):: vdot(nvar)
769
770
771! Computation of equation rates
772
773! Aggregate function
774  vdot(1) = 0
775  vdot(2) = 0
776  vdot(3) = 0
777  vdot(4) = 0
778  vdot(5) = 0
779     
780END SUBROUTINE fun
781 
782SUBROUTINE kppsolve(jvs, x)
783
784! JVS - sparse Jacobian of variables
785  REAL(kind=dp):: jvs(lu_nonzero)
786! X - Vector for variables
787  REAL(kind=dp):: x(nvar)
788
789  x(5) = x(5) / jvs(5)
790  x(4) = x(4) / jvs(4)
791  x(3) = x(3) / jvs(3)
792  x(2) = x(2) / jvs(2)
793  x(1) = x(1) / jvs(1)
794     
795END SUBROUTINE kppsolve
796 
797SUBROUTINE jac_sp(v, f, rct, jvs)
798
799! V - Concentrations of variable species (local)
800  REAL(kind=dp):: v(nvar)
801! F - Concentrations of fixed species (local)
802  REAL(kind=dp):: f(nfix)
803! RCT - Rate constants (local)
804  REAL(kind=dp):: rct(nreact)
805! JVS - sparse Jacobian of variables
806  REAL(kind=dp):: jvs(lu_nonzero)
807
808
809! Local variables
810! B - Temporary array
811  REAL(kind=dp):: b(5)
812
813! B(1) = dA(1)/dV(2)
814  b(1) = rct(1)
815! B(2) = dA(2)/dV(1)
816  b(2) = rct(2)
817! B(3) = dA(3)/dV(3)
818  b(3) = rct(3)
819! B(4) = dA(4)/dV(4)
820  b(4) = rct(4)
821! B(5) = dA(5)/dV(5)
822  b(5) = rct(5)
823
824! Construct the Jacobian terms from B's
825! JVS(1) = Jac_FULL(1,1)
826  jvs(1) = 0
827! JVS(2) = Jac_FULL(2,2)
828  jvs(2) = 0
829! JVS(3) = Jac_FULL(3,3)
830  jvs(3) = 0
831! JVS(4) = Jac_FULL(4,4)
832  jvs(4) = 0
833! JVS(5) = Jac_FULL(5,5)
834  jvs(5) = 0
835     
836END SUBROUTINE jac_sp
837 
838  elemental REAL(kind=dp)FUNCTION k_arr (k_298, tdep, temp)
839    ! arrhenius FUNCTION
840
841    REAL,    INTENT(IN):: k_298 ! k at t = 298.15k
842    REAL,    INTENT(IN):: tdep  ! temperature dependence
843    REAL(kind=dp), INTENT(IN):: temp  ! temperature
844
845    intrinsic exp
846
847    k_arr = k_298 * exp(tdep* (1._dp/temp- 3.3540e-3_dp))! 1/298.15=3.3540e-3
848
849  END FUNCTION k_arr
850 
851SUBROUTINE update_rconst()
852 INTEGER         :: k
853
854  k = is
855
856! Begin INLINED RCONST
857
858
859! End INLINED RCONST
860
861  rconst(1) = (1.0_dp)
862  rconst(2) = (1.0_dp)
863  rconst(3) = (1.0_dp)
864  rconst(4) = (1.0_dp)
865  rconst(5) = (1.0_dp)
866     
867END SUBROUTINE update_rconst
868 
869!  END FUNCTION ARR2
870REAL(kind=dp)FUNCTION arr2( a0, b0, temp)
871    REAL(kind=dp):: temp
872    REAL(kind=dp):: a0, b0
873    arr2 = a0 * exp( - b0 / temp)
874END FUNCTION arr2
875 
876SUBROUTINE initialize_kpp_ctrl(status)
877
878
879  ! i/o
880  INTEGER,         INTENT(OUT):: status
881
882  ! local
883  REAL(dp):: tsum
884  INTEGER  :: i
885
886  ! check fixed time steps
887  tsum = 0.0_dp
888  DO i=1, nmaxfixsteps
889     IF (t_steps(i)< tiny(0.0_dp))exit
890     tsum = tsum + t_steps(i)
891  ENDDO
892
893  nfsteps = i- 1
894
895  l_fixed_step = (nfsteps > 0).and.((tsum - 1.0)< tiny(0.0_dp))
896
897  IF (l_vector)THEN
898     WRITE(*,*) ' MODE             : VECTOR (LENGTH=',VL_DIM,')'
899  ELSE
900     WRITE(*,*) ' MODE             : SCALAR'
901  ENDIF
902  !
903  WRITE(*,*) ' DE-INDEXING MODE :',I_LU_DI
904  !
905  WRITE(*,*) ' ICNTRL           : ',icntrl
906  WRITE(*,*) ' RCNTRL           : ',rcntrl
907  !
908  ! note: this is ONLY meaningful for vectorized (kp4)rosenbrock- methods
909  IF (l_vector)THEN
910     IF (l_fixed_step)THEN
911        WRITE(*,*) ' TIME STEPS       : FIXED (',t_steps(1:nfsteps),')'
912     ELSE
913        WRITE(*,*) ' TIME STEPS       : AUTOMATIC'
914     ENDIF
915  ELSE
916     WRITE(*,*) ' TIME STEPS       : AUTOMATIC '//&
917          &'(t_steps (CTRL_KPP) ignored in SCALAR MODE)'
918  ENDIF
919  ! mz_pj_20070531-
920
921  status = 0
922
923
924END SUBROUTINE initialize_kpp_ctrl
925 
926SUBROUTINE error_output(c, ierr, pe)
927
928
929  INTEGER, INTENT(IN):: ierr
930  INTEGER, INTENT(IN):: pe
931  REAL(dp), DIMENSION(:), INTENT(IN):: c
932
933  write(6,*) 'ERROR in chem_gasphase_mod ',ierr,C(1)
934
935
936END SUBROUTINE error_output
937 
938      SUBROUTINE wscal(n, alpha, x, incx)
939!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
940!     constant times a vector: x(1:N) <- Alpha*x(1:N)
941!     only for incX=incY=1
942!     after BLAS
943!     replace this by the function from the optimized BLAS implementation:
944!         CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
945!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
946
947      INTEGER  :: i, incx, m, mp1, n
948      REAL(kind=dp) :: x(n), alpha
949      REAL(kind=dp), PARAMETER  :: zero=0.0_dp, one=1.0_dp
950
951      IF (alpha .eq. one)RETURN
952      IF (n .le. 0)RETURN
953
954      m = mod(n, 5)
955      IF ( m .ne. 0)THEN
956        IF (alpha .eq. (- one))THEN
957          DO i = 1, m
958            x(i) = - x(i)
959          ENDDO
960        ELSEIF (alpha .eq. zero)THEN
961          DO i = 1, m
962            x(i) = zero
963          ENDDO
964        ELSE
965          DO i = 1, m
966            x(i) = alpha* x(i)
967          ENDDO
968        ENDIF
969        IF ( n .lt. 5)RETURN
970      ENDIF
971      mp1 = m + 1
972      IF (alpha .eq. (- one))THEN
973        DO i = mp1, n, 5
974          x(i)   = - x(i)
975          x(i + 1) = - x(i + 1)
976          x(i + 2) = - x(i + 2)
977          x(i + 3) = - x(i + 3)
978          x(i + 4) = - x(i + 4)
979        ENDDO
980      ELSEIF (alpha .eq. zero)THEN
981        DO i = mp1, n, 5
982          x(i)   = zero
983          x(i + 1) = zero
984          x(i + 2) = zero
985          x(i + 3) = zero
986          x(i + 4) = zero
987        ENDDO
988      ELSE
989        DO i = mp1, n, 5
990          x(i)   = alpha* x(i)
991          x(i + 1) = alpha* x(i + 1)
992          x(i + 2) = alpha* x(i + 2)
993          x(i + 3) = alpha* x(i + 3)
994          x(i + 4) = alpha* x(i + 4)
995        ENDDO
996      ENDIF
997
998      END SUBROUTINE wscal
999 
1000      SUBROUTINE waxpy(n, alpha, x, incx, y, incy)
1001!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1002!     constant times a vector plus a vector: y <- y + Alpha*x
1003!     only for incX=incY=1
1004!     after BLAS
1005!     replace this by the function from the optimized BLAS implementation:
1006!         CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
1007!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1008
1009      INTEGER  :: i, incx, incy, m, mp1, n
1010      REAL(kind=dp):: x(n), y(n), alpha
1011      REAL(kind=dp), PARAMETER :: zero = 0.0_dp
1012
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#ifdef full_algebra   
1747   CALL  DGETRS( 'N',N ,1,A,N,Pivot,b,N,ISING)
1748   IF (info < 0)THEN
1749      print* , "error in dgetrs. ising=", ising
1750   ENDIF 
1751#else   
1752   CALL kppsolve( a, b)
1753#endif
1754
1755   istatus(nsol) = istatus(nsol) + 1
1756
1757  END SUBROUTINE ros_solve
1758
1759
1760
1761!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1762  SUBROUTINE ros2
1763!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1764! --- AN L-STABLE METHOD,2 stages,order 2
1765!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1766
1767   double precision g
1768
1769    g = 1.0_dp + 1.0_dp/sqrt(2.0_dp)
1770    rosmethod = rs2
1771!~~~> name of the method
1772    ros_Name = 'ROS-2'
1773!~~~> number of stages
1774    ros_s = 2
1775
1776!~~~> the coefficient matrices a and c are strictly lower triangular.
1777!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1778!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1779!   The general mapping formula is:
1780!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1781!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1782
1783    ros_a(1) = (1.0_dp) /g
1784    ros_c(1) = (- 2.0_dp) /g
1785!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1786!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1787    ros_newf(1) = .TRUE.
1788    ros_newf(2) = .TRUE.
1789!~~~> m_i = coefficients for new step solution
1790    ros_m(1) = (3.0_dp) /(2.0_dp* g)
1791    ros_m(2) = (1.0_dp) /(2.0_dp* g)
1792! E_i = Coefficients for error estimator
1793    ros_e(1) = 1.0_dp/(2.0_dp* g)
1794    ros_e(2) = 1.0_dp/(2.0_dp* g)
1795!~~~> ros_elo = estimator of local order - the minimum between the
1796!    main and the embedded scheme orders plus one
1797    ros_elo = 2.0_dp
1798!~~~> y_stage_i ~ y( t + h* alpha_i)
1799    ros_alpha(1) = 0.0_dp
1800    ros_alpha(2) = 1.0_dp
1801!~~~> gamma_i = \sum_j  gamma_{i, j}
1802    ros_gamma(1) = g
1803    ros_gamma(2) = -g
1804
1805 END SUBROUTINE ros2
1806
1807
1808!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1809  SUBROUTINE ros3
1810!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1811! --- AN L-STABLE METHOD,3 stages,order 3,2 function evaluations
1812!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1813
1814   rosmethod = rs3
1815!~~~> name of the method
1816   ros_Name = 'ROS-3'
1817!~~~> number of stages
1818   ros_s = 3
1819
1820!~~~> the coefficient matrices a and c are strictly lower triangular.
1821!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1822!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1823!   The general mapping formula is:
1824!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1825!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1826
1827   ros_a(1) = 1.0_dp
1828   ros_a(2) = 1.0_dp
1829   ros_a(3) = 0.0_dp
1830
1831   ros_c(1) = - 0.10156171083877702091975600115545e+01_dp
1832   ros_c(2) =  0.40759956452537699824805835358067e+01_dp
1833   ros_c(3) =  0.92076794298330791242156818474003e+01_dp
1834!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1835!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1836   ros_newf(1) = .TRUE.
1837   ros_newf(2) = .TRUE.
1838   ros_newf(3) = .FALSE.
1839!~~~> m_i = coefficients for new step solution
1840   ros_m(1) =  0.1e+01_dp
1841   ros_m(2) =  0.61697947043828245592553615689730e+01_dp
1842   ros_m(3) = - 0.42772256543218573326238373806514_dp
1843! E_i = Coefficients for error estimator
1844   ros_e(1) =  0.5_dp
1845   ros_e(2) = - 0.29079558716805469821718236208017e+01_dp
1846   ros_e(3) =  0.22354069897811569627360909276199_dp
1847!~~~> ros_elo = estimator of local order - the minimum between the
1848!    main and the embedded scheme orders plus 1
1849   ros_elo = 3.0_dp
1850!~~~> y_stage_i ~ y( t + h* alpha_i)
1851   ros_alpha(1) = 0.0_dp
1852   ros_alpha(2) = 0.43586652150845899941601945119356_dp
1853   ros_alpha(3) = 0.43586652150845899941601945119356_dp
1854!~~~> gamma_i = \sum_j  gamma_{i, j}
1855   ros_gamma(1) = 0.43586652150845899941601945119356_dp
1856   ros_gamma(2) = 0.24291996454816804366592249683314_dp
1857   ros_gamma(3) = 0.21851380027664058511513169485832e+01_dp
1858
1859  END SUBROUTINE ros3
1860
1861!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1862
1863
1864!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1865  SUBROUTINE ros4
1866!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1867!     L-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 4 STAGES
1868!     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
1869!
1870!      E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
1871!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1872!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1873!      SPRINGER-VERLAG (1990)
1874!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1875
1876
1877   rosmethod = rs4
1878!~~~> name of the method
1879   ros_Name = 'ROS-4'
1880!~~~> number of stages
1881   ros_s = 4
1882
1883!~~~> the coefficient matrices a and c are strictly lower triangular.
1884!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1885!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1886!   The general mapping formula is:
1887!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1888!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1889
1890   ros_a(1) = 0.2000000000000000e+01_dp
1891   ros_a(2) = 0.1867943637803922e+01_dp
1892   ros_a(3) = 0.2344449711399156_dp
1893   ros_a(4) = ros_a(2)
1894   ros_a(5) = ros_a(3)
1895   ros_a(6) = 0.0_dp
1896
1897   ros_c(1) = -0.7137615036412310e+01_dp
1898   ros_c(2) = 0.2580708087951457e+01_dp
1899   ros_c(3) = 0.6515950076447975_dp
1900   ros_c(4) = -0.2137148994382534e+01_dp
1901   ros_c(5) = -0.3214669691237626_dp
1902   ros_c(6) = -0.6949742501781779_dp
1903!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1904!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1905   ros_newf(1) = .TRUE.
1906   ros_newf(2) = .TRUE.
1907   ros_newf(3) = .TRUE.
1908   ros_newf(4) = .FALSE.
1909!~~~> m_i = coefficients for new step solution
1910   ros_m(1) = 0.2255570073418735e+01_dp
1911   ros_m(2) = 0.2870493262186792_dp
1912   ros_m(3) = 0.4353179431840180_dp
1913   ros_m(4) = 0.1093502252409163e+01_dp
1914!~~~> e_i  = coefficients for error estimator
1915   ros_e(1) = -0.2815431932141155_dp
1916   ros_e(2) = -0.7276199124938920e-01_dp
1917   ros_e(3) = -0.1082196201495311_dp
1918   ros_e(4) = -0.1093502252409163e+01_dp
1919!~~~> ros_elo  = estimator of local order - the minimum between the
1920!    main and the embedded scheme orders plus 1
1921   ros_elo  = 4.0_dp
1922!~~~> y_stage_i ~ y( t + h* alpha_i)
1923   ros_alpha(1) = 0.0_dp
1924   ros_alpha(2) = 0.1145640000000000e+01_dp
1925   ros_alpha(3) = 0.6552168638155900_dp
1926   ros_alpha(4) = ros_alpha(3)
1927!~~~> gamma_i = \sum_j  gamma_{i, j}
1928   ros_gamma(1) = 0.5728200000000000_dp
1929   ros_gamma(2) = -0.1769193891319233e+01_dp
1930   ros_gamma(3) = 0.7592633437920482_dp
1931   ros_gamma(4) = -0.1049021087100450_dp
1932
1933  END SUBROUTINE ros4
1934
1935!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1936  SUBROUTINE rodas3
1937!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1938! --- A STIFFLY-STABLE METHOD,4 stages,order 3
1939!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1940
1941
1942   rosmethod = rd3
1943!~~~> name of the method
1944   ros_Name = 'RODAS-3'
1945!~~~> number of stages
1946   ros_s = 4
1947
1948!~~~> the coefficient matrices a and c are strictly lower triangular.
1949!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1950!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1951!   The general mapping formula is:
1952!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1953!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1954
1955   ros_a(1) = 0.0_dp
1956   ros_a(2) = 2.0_dp
1957   ros_a(3) = 0.0_dp
1958   ros_a(4) = 2.0_dp
1959   ros_a(5) = 0.0_dp
1960   ros_a(6) = 1.0_dp
1961
1962   ros_c(1) = 4.0_dp
1963   ros_c(2) = 1.0_dp
1964   ros_c(3) = -1.0_dp
1965   ros_c(4) = 1.0_dp
1966   ros_c(5) = -1.0_dp
1967   ros_c(6) = -(8.0_dp/3.0_dp)
1968
1969!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1970!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1971   ros_newf(1) = .TRUE.
1972   ros_newf(2) = .FALSE.
1973   ros_newf(3) = .TRUE.
1974   ros_newf(4) = .TRUE.
1975!~~~> m_i = coefficients for new step solution
1976   ros_m(1) = 2.0_dp
1977   ros_m(2) = 0.0_dp
1978   ros_m(3) = 1.0_dp
1979   ros_m(4) = 1.0_dp
1980!~~~> e_i  = coefficients for error estimator
1981   ros_e(1) = 0.0_dp
1982   ros_e(2) = 0.0_dp
1983   ros_e(3) = 0.0_dp
1984   ros_e(4) = 1.0_dp
1985!~~~> ros_elo  = estimator of local order - the minimum between the
1986!    main and the embedded scheme orders plus 1
1987   ros_elo  = 3.0_dp
1988!~~~> y_stage_i ~ y( t + h* alpha_i)
1989   ros_alpha(1) = 0.0_dp
1990   ros_alpha(2) = 0.0_dp
1991   ros_alpha(3) = 1.0_dp
1992   ros_alpha(4) = 1.0_dp
1993!~~~> gamma_i = \sum_j  gamma_{i, j}
1994   ros_gamma(1) = 0.5_dp
1995   ros_gamma(2) = 1.5_dp
1996   ros_gamma(3) = 0.0_dp
1997   ros_gamma(4) = 0.0_dp
1998
1999  END SUBROUTINE rodas3
2000
2001!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2002  SUBROUTINE rodas4
2003!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2004!     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 6 STAGES
2005!
2006!      E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
2007!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
2008!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
2009!      SPRINGER-VERLAG (1996)
2010!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2011
2012
2013    rosmethod = rd4
2014!~~~> name of the method
2015    ros_Name = 'RODAS-4'
2016!~~~> number of stages
2017    ros_s = 6
2018
2019!~~~> y_stage_i ~ y( t + h* alpha_i)
2020    ros_alpha(1) = 0.000_dp
2021    ros_alpha(2) = 0.386_dp
2022    ros_alpha(3) = 0.210_dp
2023    ros_alpha(4) = 0.630_dp
2024    ros_alpha(5) = 1.000_dp
2025    ros_alpha(6) = 1.000_dp
2026
2027!~~~> gamma_i = \sum_j  gamma_{i, j}
2028    ros_gamma(1) = 0.2500000000000000_dp
2029    ros_gamma(2) = -0.1043000000000000_dp
2030    ros_gamma(3) = 0.1035000000000000_dp
2031    ros_gamma(4) = -0.3620000000000023e-01_dp
2032    ros_gamma(5) = 0.0_dp
2033    ros_gamma(6) = 0.0_dp
2034
2035!~~~> the coefficient matrices a and c are strictly lower triangular.
2036!   The lower triangular (subdiagonal) elements are stored in row-wise order:
2037!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
2038!   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
2039!                  C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
2040
2041    ros_a(1) = 0.1544000000000000e+01_dp
2042    ros_a(2) = 0.9466785280815826_dp
2043    ros_a(3) = 0.2557011698983284_dp
2044    ros_a(4) = 0.3314825187068521e+01_dp
2045    ros_a(5) = 0.2896124015972201e+01_dp
2046    ros_a(6) = 0.9986419139977817_dp
2047    ros_a(7) = 0.1221224509226641e+01_dp
2048    ros_a(8) = 0.6019134481288629e+01_dp
2049    ros_a(9) = 0.1253708332932087e+02_dp
2050    ros_a(10) = -0.6878860361058950_dp
2051    ros_a(11) = ros_a(7)
2052    ros_a(12) = ros_a(8)
2053    ros_a(13) = ros_a(9)
2054    ros_a(14) = ros_a(10)
2055    ros_a(15) = 1.0_dp
2056
2057    ros_c(1) = -0.5668800000000000e+01_dp
2058    ros_c(2) = -0.2430093356833875e+01_dp
2059    ros_c(3) = -0.2063599157091915_dp
2060    ros_c(4) = -0.1073529058151375_dp
2061    ros_c(5) = -0.9594562251023355e+01_dp
2062    ros_c(6) = -0.2047028614809616e+02_dp
2063    ros_c(7) = 0.7496443313967647e+01_dp
2064    ros_c(8) = -0.1024680431464352e+02_dp
2065    ros_c(9) = -0.3399990352819905e+02_dp
2066    ros_c(10) = 0.1170890893206160e+02_dp
2067    ros_c(11) = 0.8083246795921522e+01_dp
2068    ros_c(12) = -0.7981132988064893e+01_dp
2069    ros_c(13) = -0.3152159432874371e+02_dp
2070    ros_c(14) = 0.1631930543123136e+02_dp
2071    ros_c(15) = -0.6058818238834054e+01_dp
2072
2073!~~~> m_i = coefficients for new step solution
2074    ros_m(1) = ros_a(7)
2075    ros_m(2) = ros_a(8)
2076    ros_m(3) = ros_a(9)
2077    ros_m(4) = ros_a(10)
2078    ros_m(5) = 1.0_dp
2079    ros_m(6) = 1.0_dp
2080
2081!~~~> e_i  = coefficients for error estimator
2082    ros_e(1) = 0.0_dp
2083    ros_e(2) = 0.0_dp
2084    ros_e(3) = 0.0_dp
2085    ros_e(4) = 0.0_dp
2086    ros_e(5) = 0.0_dp
2087    ros_e(6) = 1.0_dp
2088
2089!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2090!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2091    ros_newf(1) = .TRUE.
2092    ros_newf(2) = .TRUE.
2093    ros_newf(3) = .TRUE.
2094    ros_newf(4) = .TRUE.
2095    ros_newf(5) = .TRUE.
2096    ros_newf(6) = .TRUE.
2097
2098!~~~> ros_elo  = estimator of local order - the minimum between the
2099!        main and the embedded scheme orders plus 1
2100    ros_elo = 4.0_dp
2101
2102  END SUBROUTINE rodas4
2103!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2104  SUBROUTINE rang3
2105!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2106! STIFFLY-STABLE W METHOD OF ORDER 3,WITH 4 STAGES
2107!
2108! J. RANG and L. ANGERMANN
2109! NEW ROSENBROCK W-METHODS OF ORDER 3
2110! FOR PARTIAL DIFFERENTIAL ALGEBRAIC
2111!        EQUATIONS OF INDEX 1                                             
2112! BIT Numerical Mathematics (2005) 45: 761-787
2113!  DOI: 10.1007/s10543-005-0035-y
2114! Table 4.1-4.2
2115!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2116
2117
2118    rosmethod = rg3
2119!~~~> name of the method
2120    ros_Name = 'RANG-3'
2121!~~~> number of stages
2122    ros_s = 4
2123
2124    ros_a(1) = 5.09052051067020d+00;
2125    ros_a(2) = 5.09052051067020d+00;
2126    ros_a(3) = 0.0d0;
2127    ros_a(4) = 4.97628111010787d+00;
2128    ros_a(5) = 2.77268164715849d-02;
2129    ros_a(6) = 2.29428036027904d-01;
2130
2131    ros_c(1) = - 1.16790812312283d+01;
2132    ros_c(2) = - 1.64057326467367d+01;
2133    ros_c(3) = - 2.77268164715850d-01;
2134    ros_c(4) = - 8.38103960500476d+00;
2135    ros_c(5) = - 8.48328409199343d-01;
2136    ros_c(6) =  2.87009860433106d-01;
2137
2138    ros_m(1) =  5.22582761233094d+00;
2139    ros_m(2) = - 5.56971148154165d-01;
2140    ros_m(3) =  3.57979469353645d-01;
2141    ros_m(4) =  1.72337398521064d+00;
2142
2143    ros_e(1) = - 5.16845212784040d+00;
2144    ros_e(2) = - 1.26351942603842d+00;
2145    ros_e(3) = - 1.11022302462516d-16;
2146    ros_e(4) =  2.22044604925031d-16;
2147
2148    ros_alpha(1) = 0.0d00;
2149    ros_alpha(2) = 2.21878746765329d+00;
2150    ros_alpha(3) = 2.21878746765329d+00;
2151    ros_alpha(4) = 1.55392337535788d+00;
2152
2153    ros_gamma(1) =  4.35866521508459d-01;
2154    ros_gamma(2) = - 1.78292094614483d+00;
2155    ros_gamma(3) = - 2.46541900496934d+00;
2156    ros_gamma(4) = - 8.05529997906370d-01;
2157
2158
2159!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2160!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2161    ros_newf(1) = .TRUE.
2162    ros_newf(2) = .TRUE.
2163    ros_newf(3) = .TRUE.
2164    ros_newf(4) = .TRUE.
2165
2166!~~~> ros_elo  = estimator of local order - the minimum between the
2167!        main and the embedded scheme orders plus 1
2168    ros_elo = 3.0_dp
2169
2170  END SUBROUTINE rang3
2171!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2172
2173!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2174!   End of the set of internal Rosenbrock subroutines
2175!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2176END SUBROUTINE rosenbrock
2177 
2178SUBROUTINE funtemplate( t, y, ydot)
2179!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2180!  Template for the ODE function call.
2181!  Updates the rate coefficients (and possibly the fixed species) at each call
2182!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2183!~~~> input variables
2184   REAL(kind=dp):: t, y(nvar)
2185!~~~> output variables
2186   REAL(kind=dp):: ydot(nvar)
2187!~~~> local variables
2188   REAL(kind=dp):: told
2189
2190   told = time
2191   time = t
2192   CALL fun( y, fix, rconst, ydot)
2193   time = told
2194
2195END SUBROUTINE funtemplate
2196 
2197SUBROUTINE jactemplate( t, y, jcb)
2198!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2199!  Template for the ODE Jacobian call.
2200!  Updates the rate coefficients (and possibly the fixed species) at each call
2201!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2202!~~~> input variables
2203    REAL(kind=dp):: t, y(nvar)
2204!~~~> output variables
2205#ifdef full_algebra   
2206    REAL(kind=dp):: jv(lu_nonzero), jcb(nvar, nvar)
2207#else
2208    REAL(kind=dp):: jcb(lu_nonzero)
2209#endif   
2210!~~~> local variables
2211    REAL(kind=dp):: told
2212#ifdef full_algebra   
2213    INTEGER :: i, j
2214#endif   
2215
2216    told = time
2217    time = t
2218#ifdef full_algebra   
2219    CALL jac_sp(y, fix, rconst, jv)
2220    DO j=1, nvar
2221      DO i=1, nvar
2222         jcb(i, j) = 0.0_dp
2223      ENDDO
2224    ENDDO
2225    DO i=1, lu_nonzero
2226       jcb(lu_irow(i), lu_icol(i)) = jv(i)
2227    ENDDO
2228#else
2229    CALL jac_sp( y, fix, rconst, jcb)
2230#endif   
2231    time = told
2232
2233END SUBROUTINE jactemplate
2234 
2235  SUBROUTINE kppdecomp( jvs, ier)                                   
2236  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2237  !        sparse lu factorization                                   
2238  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2239  ! loop expansion generated by kp4                                   
2240                                                                     
2241    INTEGER  :: ier                                                   
2242    REAL(kind=dp):: jvs(lu_nonzero), w(nvar), a             
2243    INTEGER  :: k, kk, j, jj                                         
2244                                                                     
2245    a = 0.                                                           
2246    ier = 0                                                           
2247                                                                     
2248!   i = 1
2249!   i = 2
2250!   i = 3
2251!   i = 4
2252!   i = 5
2253    RETURN                                                           
2254                                                                     
2255  END SUBROUTINE kppdecomp                                           
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                                :: istatf                 
2282  INTEGER                                :: vl_dim_lo               
2283                                                                   
2284                                                                   
2285  IF (PRESENT (istatus)) istatus = 0                             
2286  IF (PRESENT (icntrl_i)) icntrl  = icntrl_i                     
2287  IF (PRESENT (rcntrl_i)) rcntrl  = rcntrl_i                     
2288                                                                   
2289  vl_glo = size(tempi, 1)                                           
2290                                                                   
2291  vl_dim_lo = vl_dim                                               
2292  DO k=1, vl_glo, vl_dim_lo                                           
2293    is = k                                                         
2294    ie = min(k+ vl_dim_lo-1, vl_glo)                                 
2295    vl = ie-is+ 1                                                   
2296                                                                   
2297    c(:) = conc(is, :)                                           
2298                                                                   
2299    temp = tempi(is)                                             
2300                                                                   
2301    qvap = qvapi(is)                                             
2302                                                                   
2303    fakt = fakti(is)                                             
2304                                                                   
2305    CALL initialize                                                 
2306                                                                   
2307    phot(:) = photo(is, :)                                           
2308                                                                   
2309    CALL update_rconst                                             
2310                                                                   
2311    dt = time_step_len                                             
2312                                                                   
2313    ! integrate from t=0 to t=dt                                   
2314    CALL integrate(0._dp, dt, icntrl, rcntrl, istatus_u = istatus_u, ierr_u=ierr_u)
2315                                                                   
2316                                                                   
2317   IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN                       
2318      IF (l_debug) CALL error_output(conc(is, :), ierr_u, pe)         
2319   ENDIF                                                             
2320                                                                     
2321    conc(is, :) = c(:)                                               
2322                                                                   
2323    ! RETURN diagnostic information                                 
2324                                                                   
2325    IF (PRESENT(ierrf))   ierrf(is) = ierr_u                     
2326    IF (PRESENT(xnacc))   xnacc(is) = istatus_u(4)               
2327    IF (PRESENT(xnrej))   xnrej(is) = istatus_u(5)               
2328                                                                   
2329    IF (PRESENT (istatus)) THEN                                   
2330      istatus(1:8) = istatus(1:8) + istatus_u(1:8)               
2331    ENDIF                                                         
2332                                                                   
2333  END DO                                                           
2334 
2335                                                                   
2336! Deallocate input arrays                                           
2337                                                                   
2338                                                                   
2339  data_loaded = .FALSE.                                             
2340                                                                   
2341  RETURN                                                           
2342END SUBROUTINE chem_gasphase_integrate                                       
2343
2344END MODULE chem_gasphase_mod
2345 
Note: See TracBrowser for help on using the repository browser.