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

Last change on this file since 4370 was 4370, checked in by raasch, 2 years ago

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

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