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

Last change on this file since 3820 was 3820, checked in by forkel, 3 years ago

renaming of get_mechanismname, do_emiss and do_depo, sorting in chem_modules

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