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

Last change on this file since 4370 was 4370, checked in by raasch, 23 months ago

bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added

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