source: palm/trunk/SOURCE/chem_gasphase_mod.f90 @ 3797

Last change on this file since 3797 was 3797, checked in by forkel, 2 years ago

Modifications for OpenMP version by Klaus Ketelsen

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