source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_phstatp/chem_gasphase_mod.f90

Last change on this file was 4841, checked in by forkel, 4 years ago

updated copyright statements for chem_gasphase_mod.f90. This must be done in UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header and NOT in SOURCE/chem_gasphase_mod.f90.

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