Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Oct 2, 2018 12:21:11 PM (6 years ago)
Author:
kanani
Message:

Merge chemistry branch at r3297 to trunk

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/chem_gasphase_mod.f90

    r2766 r3298  
    11MODULE chem_gasphase_mod
     2 
     3!   Mechanism: passive
     4!
    25!------------------------------------------------------------------------------!
    36!
    47! ******Module chem_gasphase_mod is automatically generated by kpp4palm ******
    58!
    6 !           *********Please do NOT change this Code *********
     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.
    727!
    828!------------------------------------------------------------------------------!
    … …  
    1131! PALM is free software: you can redistribute it and/or modify it under the
    1232! terms of the GNU General Public License as published by the Free Software
    13 ! Foundation, either version 3 of the License, or (at your option) any later
     33! Foundation,either version 3 of the License,or (at your option) any later
    1434! version.
    1535!
    16 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
     36! PALM is distributed in the hope that it will be useful,but WITHOUT ANY
    1737! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    1838! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    1939!
    2040! You should have received a copy of the GNU General Public License along with
    21 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     41! PALM. If not,see <http://www.gnu.org/licenses/>.
    2242!
    2343! Copyright 1997-2018 Leibniz Universitaet Hannover
    … …  
    3050! Former revisions:
    3151! -----------------
    32 ! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
    33 ! Removed preprocessor directive __chem
    34 !
    35 ! 2460 2017-09-13 14:47:48Z forkel
    36 !
    37 !
    38 ! Variables for photolyis added
    39 !
    40 !
    41 !
     52! $Id: chem_gasphase_mod.f90 3287 2018-09-28 10:19:58Z forkel $
     53! forkel June 2018: qvap,fakt added
     54! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
     55!
     56!
     57! 3287 2018-09-28 10:19:58Z forkel
     58!
     59! forkel Sept. 2017: Variables for photolyis added
    4260!
    4361!
    … …  
    5270  USE kinds,           ONLY: dp=>wp
    5371
    54   USE pegrid,          ONLY: myid,threads_per_task
     72  USE pegrid,          ONLY: myid, threads_per_task
    5573
    5674  IMPLICIT        NONE
    5775  PRIVATE
    58   !SAVE  ! NOTE: OCCURS AGAIN IN AUTOMATICALLY GENERATED CODE ...
     76  !SAVE  ! note: occurs again in automatically generated code ...
    5977
    6078!  PUBLIC :: IERR_NAMES
    … …  
    6381!         ,REQ_MCFCT,IP_MAX,jname
    6482
    65   PUBLIC :: eqn_names, phot_names,spc_names
     83  PUBLIC :: eqn_names, phot_names, spc_names
    6684  PUBLIC :: nmaxfixsteps
    67   PUBLIC :: atol,rtol
    68   PUBLIC :: nspec,nreact
     85  PUBLIC :: atol, rtol
     86  PUBLIC :: nspec, nreact
    6987  PUBLIC :: temp
     88  PUBLIC :: qvap
     89  PUBLIC :: fakt
    7090  PUBLIC :: phot
    7191  PUBLIC :: rconst
    7292  PUBLIC :: nvar
    7393  PUBLIC :: nphot
    74  
    75   PUBLIC :: initialize,integrate,update_rconst
     94  PUBLIC :: vl_dim                     ! PUBLIC to ebable other MODULEs to distiguish between scalar and vec
     95 
     96  PUBLIC :: initialize, integrate, update_rconst
    7697  PUBLIC :: chem_gasphase_integrate
    7798  PUBLIC :: initialize_kpp_ctrl
    … …  
    82103                                                                 
    83104  LOGICAL, PARAMETER          :: l_vector = .FALSE.           
    84   INTEGER, PARAMETER          :: i_lu_di = 0
     105  INTEGER, PARAMETER          :: i_lu_di = 2
    85106  INTEGER, PARAMETER          :: vl_dim = 1
    86107  INTEGER                     :: vl                             
    87108                                                                 
    88109  INTEGER                     :: vl_glo                         
    89   INTEGER                     :: is,ie                           
     110  INTEGER                     :: is, ie                           
    90111                                                                 
    91   INTEGER, DIMENSION(vl_dim)  :: kacc,krej                     
     112                                                                 
     113  INTEGER, DIMENSION(vl_dim)  :: kacc, krej                       
    92114  INTEGER, DIMENSION(vl_dim)  :: ierrv                           
    93   LOGICAL                     :: data_loaded = .false.             
     115  LOGICAL                     :: data_loaded = .FALSE.             
    94116! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    95117!
    … …  
    107129!
    108130! File                 : chem_gasphase_mod_Parameters.f90
    109 ! Time                 : Fri Dec  8 11:54:15 2017
    110 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
     131! Time                 : Tue Sep 25 18:34:57 2018
     132! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    111133! Equation file        : chem_gasphase_mod.kpp
    112134! Output root filename : chem_gasphase_mod
    … …  
    120142
    121143! NSPEC - Number of chemical species
    122   INTEGER,PARAMETER :: nspec = 2
     144  INTEGER, PARAMETER :: nspec = 2
    123145! NVAR - Number of Variable species
    124   INTEGER,PARAMETER :: nvar = 2
     146  INTEGER, PARAMETER :: nvar = 2
    125147! NVARACT - Number of Active species
    126   INTEGER,PARAMETER :: nvaract = 2
     148  INTEGER, PARAMETER :: nvaract = 2
    127149! NFIX - Number of Fixed species
    128   INTEGER,PARAMETER :: nfix = 1
     150  INTEGER, PARAMETER :: nfix = 1
    129151! NREACT - Number of reactions
    130   INTEGER,PARAMETER :: nreact = 2
     152  INTEGER, PARAMETER :: nreact = 2
    131153! NVARST - Starting of variables in conc. vect.
    132   INTEGER,PARAMETER :: nvarst = 1
     154  INTEGER, PARAMETER :: nvarst = 1
    133155! NFIXST - Starting of fixed in conc. vect.
    134   INTEGER,PARAMETER :: nfixst = 3
     156  INTEGER, PARAMETER :: nfixst = 3
    135157! NONZERO - Number of nonzero entries in Jacobian
    136   INTEGER,PARAMETER :: nonzero = 2
     158  INTEGER, PARAMETER :: nonzero = 2
    137159! LU_NONZERO - Number of nonzero entries in LU factoriz. of Jacobian
    138   INTEGER,PARAMETER :: lu_nonzero = 2
     160  INTEGER, PARAMETER :: lu_nonzero = 2
    139161! CNVAR - (NVAR+1) Number of elements in compressed row format
    140   INTEGER,PARAMETER :: cnvar = 3
     162  INTEGER, PARAMETER :: cnvar = 3
    141163! CNEQN - (NREACT+1) Number stoicm elements in compressed col format
    142   INTEGER,PARAMETER :: cneqn = 3
     164  INTEGER, PARAMETER :: cneqn = 3
    143165! NHESS - Length of Sparse Hessian
    144   INTEGER,PARAMETER :: nhess = 1
    145 ! NLOOKAT - Number of species to look at
    146   INTEGER,PARAMETER :: nlookat = 0
    147 ! NMONITOR - Number of species to monitor
    148   INTEGER,PARAMETER :: nmonitor = 0
     166  INTEGER, PARAMETER :: nhess = 1
    149167! NMASS - Number of atoms to check mass balance
    150   INTEGER,PARAMETER :: nmass = 1
     168  INTEGER, PARAMETER :: nmass = 1
    151169
    152170! Index declaration for variable species in C and VAR
    153171!   VAR(ind_spc) = C(ind_spc)
    154172
    155   INTEGER,PARAMETER,PUBLIC :: ind_pm10 = 1
    156   INTEGER,PARAMETER,PUBLIC :: ind_pm25 = 2
     173  INTEGER, PARAMETER, PUBLIC :: ind_pm10 = 1
     174  INTEGER, PARAMETER, PUBLIC :: ind_pm25 = 2
    157175
    158176! Index declaration for fixed species in C
    … …  
    165183
    166184! NJVRP - Length of sparse Jacobian JVRP
    167   INTEGER,PARAMETER :: njvrp = 2
     185  INTEGER, PARAMETER :: njvrp = 2
    168186
    169187! NSTOICM - Length of Sparse Stoichiometric Matrix
    170   INTEGER,PARAMETER :: nstoicm = 1
     188  INTEGER, PARAMETER :: nstoicm = 1
    171189
    172190
    … …  
    186204!
    187205! File                 : chem_gasphase_mod_Global.f90
    188 ! Time                 : Fri Dec  8 11:54:15 2017
    189 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
     206! Time                 : Tue Sep 25 18:34:57 2018
     207! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    190208! Equation file        : chem_gasphase_mod.kpp
    191209! Output root filename : chem_gasphase_mod
    … …  
    207225  REAL(kind=dp):: fix(nfix)
    208226! VAR,FIX are chunks of array C
    209       equivalence( c(1),var(1))
     227      EQUIVALENCE( c(1), var(1))
    210228! RCONST - Rate constants (global)
    211229  REAL(kind=dp):: rconst(nreact)
    212230! TIME - Current integration time
    213231  REAL(kind=dp):: time
    214 ! SUN - Sunlight intensity between [0,1]
    215   REAL(kind=dp):: sun
    216232! TEMP - Temperature
    217   REAL(dp),dimension(:),allocatable             :: temp
    218 ! RTOLS - (scalar) Relative tolerance
    219   REAL(kind=dp):: rtols
     233  REAL(kind=dp):: temp
    220234! TSTART - Integration start time
    221235  REAL(kind=dp):: tstart
    222 ! TEND - Integration end time
    223   REAL(kind=dp):: tend
    224 ! DT - Integration step
    225   REAL(kind=dp):: dt
    226236! ATOL - Absolute tolerance
    227237  REAL(kind=dp):: atol(nvar)
    … …  
    230240! STEPMIN - Lower bound for integration step
    231241  REAL(kind=dp):: stepmin
    232 ! STEPMAX - Upper bound for integration step
    233   REAL(kind=dp):: stepmax
    234242! CFACTOR - Conversion factor for concentration units
    235243  REAL(kind=dp):: cfactor
    236 ! DDMTYPE - DDM sensitivity w.r.t.: 0=init.val.,1=params
    237   INTEGER :: ddmtype
    238244
    239245! INLINED global variable declarations
    240246
    241   !   declaration of global variable declarations for photolysis will come from
     247! QVAP - Water vapor
     248  REAL(kind=dp):: qvap
     249! FAKT - Conversion factor
     250  REAL(kind=dp):: fakt
     251
    242252
    243253! INLINED global variable declarations
    … …  
    260270!
    261271! File                 : chem_gasphase_mod_JacobianSP.f90
    262 ! Time                 : Fri Dec  8 11:54:15 2017
    263 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
     272! Time                 : Tue Sep 25 18:34:57 2018
     273! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    264274! Equation file        : chem_gasphase_mod.kpp
    265275! Output root filename : chem_gasphase_mod
    … …  
    275285
    276286
    277   INTEGER,PARAMETER,DIMENSION(2):: lu_irow =  (/ &
     287  INTEGER, PARAMETER, DIMENSION(2):: lu_irow =  (/ &
    278288       1, 2 /)
    279289
    280   INTEGER,PARAMETER,DIMENSION(2):: lu_icol =  (/ &
     290  INTEGER, PARAMETER, DIMENSION(2):: lu_icol =  (/ &
    281291       1, 2 /)
    282292
    283   INTEGER,PARAMETER,DIMENSION(3):: lu_crow =  (/ &
     293  INTEGER, PARAMETER, DIMENSION(3):: lu_crow =  (/ &
    284294       1, 2, 3 /)
    285295
    286   INTEGER,PARAMETER,DIMENSION(3):: lu_diag =  (/ &
     296  INTEGER, PARAMETER, DIMENSION(3):: lu_diag =  (/ &
    287297       1, 2, 3 /)
    288298
    … …  
    304314!
    305315! File                 : chem_gasphase_mod_Monitor.f90
    306 ! Time                 : Fri Dec  8 11:54:15 2017
    307 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
     316! Time                 : Tue Sep 25 18:34:57 2018
     317! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    308318! Equation file        : chem_gasphase_mod.kpp
    309319! Output root filename : chem_gasphase_mod
    … …  
    315325
    316326
    317   CHARACTER(len=15),PARAMETER,DIMENSION(2):: spc_names =  (/ &
     327  CHARACTER(len=15), PARAMETER, DIMENSION(2):: spc_names =  (/ &
    318328     'PM10           ','PM25           ' /)
    319329
    320   INTEGER,DIMENSION(1):: lookat
    321   INTEGER,DIMENSION(1):: monitor
    322   CHARACTER(len=15),DIMENSION(1):: smass
    323   CHARACTER(len=100),PARAMETER,DIMENSION(2):: eqn_names =  (/ &
     330  CHARACTER(len=100), PARAMETER, DIMENSION(2):: eqn_names =  (/ &
    324331     'PM10 --> PM10                                                                                       ',&
    325332     'PM25 --> PM25                                                                                       ' /)
    … …  
    329336  !   inline f90_data: declaration of global variables for photolysis
    330337  !   REAL(kind=dp):: phot(nphot)must eventually be moved to global later for
    331   INTEGER,PARAMETER :: nphot = 1
     338  INTEGER, PARAMETER :: nphot = 1
    332339  !   phot photolysis frequencies
    333340  REAL(kind=dp):: phot(nphot)
    334341
    335   INTEGER,PARAMETER,PUBLIC :: j_no2 = 1
    336 
    337   CHARACTER(len=15),PARAMETER,DIMENSION(nphot):: phot_names =   (/ &
     342  INTEGER, PARAMETER, PUBLIC :: j_no2 = 1
     343
     344  CHARACTER(len=15), PARAMETER, DIMENSION(nphot):: phot_names =   (/ &
    338345     'J_NO2          '/)
    339346
    … …  
    367374!
    368375! File                 : chem_gasphase_mod_Initialize.f90
    369 ! Time                 : Fri Dec  8 11:54:15 2017
    370 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
     376! Time                 : Tue Sep 25 18:34:57 2018
     377! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    371378! Equation file        : chem_gasphase_mod.kpp
    372379! Output root filename : chem_gasphase_mod
    … …  
    393400!
    394401! File                 : chem_gasphase_mod_Integrator.f90
    395 ! Time                 : Fri Dec  8 11:54:15 2017
    396 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
     402! Time                 : Tue Sep 25 18:34:57 2018
     403! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    397404! Equation file        : chem_gasphase_mod.kpp
    398405! Output root filename : chem_gasphase_mod
    … …  
    432439 
    433440!~~~>  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
     441  INTEGER, PARAMETER :: nfun=1, njac=2, nstp=3, nacc=4, &
     442                        nrej=5, ndec=6, nsol=7, nsng=8, &
     443                        ntexit=1, nhexit=2, nhnew = 3
    437444
    438445! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    451458!
    452459! File                 : chem_gasphase_mod_LinearAlgebra.f90
    453 ! Time                 : Fri Dec  8 11:54:15 2017
    454 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
     460! Time                 : Tue Sep 25 18:34:57 2018
     461! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    455462! Equation file        : chem_gasphase_mod.kpp
    456463! Output root filename : chem_gasphase_mod
    … …  
    478485!
    479486! File                 : chem_gasphase_mod_Jacobian.f90
    480 ! Time                 : Fri Dec  8 11:54:15 2017
    481 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
     487! Time                 : Tue Sep 25 18:34:57 2018
     488! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    482489! Equation file        : chem_gasphase_mod.kpp
    483490! Output root filename : chem_gasphase_mod
    … …  
    505512!
    506513! File                 : chem_gasphase_mod_Function.f90
    507 ! Time                 : Fri Dec  8 11:54:15 2017
    508 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
     514! Time                 : Tue Sep 25 18:34:57 2018
     515! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    509516! Equation file        : chem_gasphase_mod.kpp
    510517! Output root filename : chem_gasphase_mod
    … …  
    534541!
    535542! File                 : chem_gasphase_mod_Rates.f90
    536 ! Time                 : Fri Dec  8 11:54:15 2017
    537 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
     543! Time                 : Tue Sep 25 18:34:57 2018
     544! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    538545! Equation file        : chem_gasphase_mod.kpp
    539546! Output root filename : chem_gasphase_mod
    … …  
    560567!
    561568! File                 : chem_gasphase_mod_Util.f90
    562 ! Time                 : Fri Dec  8 11:54:15 2017
    563 ! Working directory    : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
     569! Time                 : Tue Sep 25 18:34:57 2018
     570! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
    564571! Equation file        : chem_gasphase_mod.kpp
    565572! Output root filename : chem_gasphase_mod
    … …  
    575582
    576583  ! 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  ! - l_vector is automatically defined by kp4
     585  ! - vl_dim is automatically defined by kp4
     586  ! - i_lu_di is automatically defined by kp4
     587  ! - wanted is automatically defined by xmecca
     588  ! - icntrl rcntrl are automatically defined by kpp
     589  ! - "USE messy_main_tools" is in MODULE_header of messy_mecca_kpp.f90
     590  ! - SAVE will be automatically added by kp4
    584591
    585592  !SAVE
    … …  
    587594  ! for fixed time step control
    588595  ! ... max. number of fixed time steps (sum must be 1)
    589   INTEGER,PARAMETER         :: nmaxfixsteps = 50
     596  INTEGER, PARAMETER         :: nmaxfixsteps = 50
    590597  ! ... switch for fixed time stepping
    591   LOGICAL,PUBLIC            :: l_fixed_step = .false.
    592   INTEGER,PUBLIC            :: nfsteps = 1
     598  LOGICAL, PUBLIC            :: l_fixed_step = .FALSE.
     599  INTEGER, PUBLIC            :: nfsteps = 1
    593600  ! ... number of kpp control PARAMETERs
    594   INTEGER,PARAMETER,PUBLIC :: nkppctrl = 20
     601  INTEGER, PARAMETER, PUBLIC :: nkppctrl = 20
    595602  !
    596   INTEGER, DIMENSION(nkppctrl),PUBLIC     :: icntrl = 0
    597   REAL(dp),DIMENSION(nkppctrl),PUBLIC     :: rcntrl = 0.0_dp
    598   REAL(dp),DIMENSION(nmaxfixsteps),PUBLIC :: t_steps = 0.0_dp
     603  INTEGER, DIMENSION(nkppctrl), PUBLIC     :: icntrl = 0
     604  REAL(dp), DIMENSION(nkppctrl), PUBLIC     :: rcntrl = 0.0_dp
     605  REAL(dp), DIMENSION(nmaxfixsteps), PUBLIC :: t_steps = 0.0_dp
    599606
    600607  ! END header MODULE initialize_kpp_ctrl_template
    … …  
    619626  END INTERFACE        kppsolve
    620627 
     628  INTERFACE            jac_sp
     629    MODULE PROCEDURE   jac_sp
     630  END INTERFACE        jac_sp
     631 
     632  INTERFACE            k_arr
     633    MODULE PROCEDURE   k_arr
     634  END INTERFACE        k_arr
     635 
     636  INTERFACE            update_rconst
     637    MODULE PROCEDURE   update_rconst
     638  END INTERFACE        update_rconst
     639 
     640  INTERFACE            arr2
     641    MODULE PROCEDURE   arr2
     642  END INTERFACE        arr2
     643 
     644  INTERFACE            initialize_kpp_ctrl
     645    MODULE PROCEDURE   initialize_kpp_ctrl
     646  END INTERFACE        initialize_kpp_ctrl
     647 
     648  INTERFACE            error_output
     649    MODULE PROCEDURE   error_output
     650  END INTERFACE        error_output
     651 
     652  INTERFACE            wscal
     653    MODULE PROCEDURE   wscal
     654  END INTERFACE        wscal
     655 
     656!INTERFACE not working  INTERFACE            waxpy
     657!INTERFACE not working    MODULE PROCEDURE   waxpy
     658!INTERFACE not working  END INTERFACE        waxpy
     659 
     660  INTERFACE            rosenbrock
     661    MODULE PROCEDURE   rosenbrock
     662  END INTERFACE        rosenbrock
     663 
     664  INTERFACE            funtemplate
     665    MODULE PROCEDURE   funtemplate
     666  END INTERFACE        funtemplate
     667 
     668  INTERFACE            jactemplate
     669    MODULE PROCEDURE   jactemplate
     670  END INTERFACE        jactemplate
     671 
    621672  INTERFACE            kppdecomp
    622673    MODULE PROCEDURE   kppdecomp
    623674  END INTERFACE        kppdecomp
    624675 
    625   INTERFACE            wlamch
    626     MODULE PROCEDURE   wlamch
    627   END INTERFACE        wlamch
    628  
    629   INTERFACE            wlamch_add
    630     MODULE PROCEDURE   wlamch_add
    631   END INTERFACE        wlamch_add
    632  
    633   INTERFACE            jac_sp
    634     MODULE PROCEDURE   jac_sp
    635   END INTERFACE        jac_sp
    636  
    637   INTERFACE            k_arr
    638     MODULE PROCEDURE   k_arr
    639   END INTERFACE        k_arr
    640  
    641   INTERFACE            update_rconst
    642     MODULE PROCEDURE   update_rconst
    643   END INTERFACE        update_rconst
    644  
    645   INTERFACE            arr2
    646     MODULE PROCEDURE   arr2
    647   END INTERFACE        arr2
    648  
    649   INTERFACE            initialize_kpp_ctrl
    650     MODULE PROCEDURE   initialize_kpp_ctrl
    651   END INTERFACE        initialize_kpp_ctrl
    652  
    653   INTERFACE            error_output
    654     MODULE PROCEDURE   error_output
    655   END INTERFACE        error_output
    656  
    657 !interface not working  INTERFACE            wcopy
    658 !interface not working    MODULE PROCEDURE   wcopy
    659 !interface not working  END INTERFACE        wcopy
    660  
    661   INTERFACE            wscal
    662     MODULE PROCEDURE   wscal
    663   END INTERFACE        wscal
    664  
    665 !interface not working  INTERFACE            waxpy
    666 !interface not working    MODULE PROCEDURE   waxpy
    667 !interface not working  END INTERFACE        waxpy
    668  
    669   INTERFACE            rosenbrock
    670     MODULE PROCEDURE   rosenbrock
    671   END INTERFACE        rosenbrock
    672  
    673   INTERFACE            funtemplate
    674     MODULE PROCEDURE   funtemplate
    675   END INTERFACE        funtemplate
    676  
    677   INTERFACE            jactemplate
    678     MODULE PROCEDURE   jactemplate
    679   END INTERFACE        jactemplate
    680  
    681676  INTERFACE            chem_gasphase_integrate
    682677    MODULE PROCEDURE   chem_gasphase_integrate
    683678  END INTERFACE        chem_gasphase_integrate
    684679 
    685   INTERFACE            fill_temp
    686     MODULE PROCEDURE   fill_temp
    687   END INTERFACE        fill_temp
    688   PUBLIC               fill_temp
    689  
    690680 
    691681 CONTAINS
    … …  
    694684
    695685
     686  INTEGER         :: j, k
    696687
    697688  INTEGER :: i
    698689  REAL(kind=dp):: x
    699 
     690  k = is
    700691  cfactor = 1.000000e+00_dp
    701692
    702   x = (0.)*cfactor
    703   DO i = 1,nvar
    704     var(i) = x
     693  x = (0.) * cfactor
     694  DO i = 1 , nvar
    705695  ENDDO
    706696
    707   x = (0.)*cfactor
    708   DO i = 1,nfix
     697  x = (0.) * cfactor
     698  DO i = 1 , nfix
    709699    fix(i) = x
    710700  ENDDO
    … …  
    720710END SUBROUTINE initialize
    721711 
    722 SUBROUTINE integrate( tin,tout,&
    723   icntrl_u,rcntrl_u,istatus_u,rstatus_u,ierr_u)
    724 
    725 
    726    REAL(kind=dp),INTENT(in):: tin  ! start time
    727    REAL(kind=dp),INTENT(in):: tout ! END time
     712SUBROUTINE integrate( tin, tout, &
     713  icntrl_u, rcntrl_u, istatus_u, rstatus_u, ierr_u)
     714
     715
     716   REAL(kind=dp), INTENT(IN):: tin  ! start time
     717   REAL(kind=dp), INTENT(IN):: tout ! END time
    728718   ! OPTIONAL input PARAMETERs and statistics
    729    INTEGER,      INTENT(in), OPTIONAL :: icntrl_u(20)
    730    REAL(kind=dp),INTENT(in), OPTIONAL :: rcntrl_u(20)
    731    INTEGER,      INTENT(out),OPTIONAL :: istatus_u(20)
    732    REAL(kind=dp),INTENT(out),OPTIONAL :: rstatus_u(20)
    733    INTEGER,      INTENT(out),OPTIONAL :: ierr_u
    734 
    735    REAL(kind=dp):: rcntrl(20),rstatus(20)
    736    INTEGER       :: icntrl(20),istatus(20),ierr
    737 
    738    INTEGER,SAVE :: ntotal = 0
     719   INTEGER,      INTENT(IN), OPTIONAL :: icntrl_u(20)
     720   REAL(kind=dp), INTENT(IN), OPTIONAL :: rcntrl_u(20)
     721   INTEGER,      INTENT(OUT), OPTIONAL :: istatus_u(20)
     722   REAL(kind=dp), INTENT(OUT), OPTIONAL :: rstatus_u(20)
     723   INTEGER,      INTENT(OUT), OPTIONAL :: ierr_u
     724
     725   REAL(kind=dp):: rcntrl(20), rstatus(20)
     726   INTEGER       :: icntrl(20), istatus(20), ierr
     727
     728   INTEGER, SAVE :: ntotal = 0
    739729
    740730   icntrl(:) = 0
    … …  
    744734
    745735    !~~~> fine-tune the integrator:
    746    icntrl(1) = 0        ! 0 -  non- autonomous,1 - autonomous
    747    icntrl(2) = 0        ! 0 -  vector tolerances,1 - scalars
    748 
    749    ! IF OPTIONAL PARAMETERs are given,and IF they are >0,
     736   icntrl(1) = 0      ! 0 - non- autonomous, 1 - autonomous
     737   icntrl(2) = 0      ! 0 - vector tolerances, 1 - scalars
     738
     739   ! IF OPTIONAL PARAMETERs are given, and IF they are >0,
    750740   ! THEN they overwrite default settings.
    751    IF (present(icntrl_u))THEN
    752      where(icntrl_u(:)> 0)icntrl(:) = icntrl_u(:)
     741   IF (PRESENT(icntrl_u))THEN
     742     WHERE(icntrl_u(:)> 0)icntrl(:) = icntrl_u(:)
    753743   ENDIF
    754    IF (present(rcntrl_u))THEN
    755      where(rcntrl_u(:)> 0)rcntrl(:) = rcntrl_u(:)
     744   IF (PRESENT(rcntrl_u))THEN
     745     WHERE(rcntrl_u(:)> 0)rcntrl(:) = rcntrl_u(:)
    756746   ENDIF
    757747
    758748
    759    CALL rosenbrock(nvar,var,tin,tout,  &
    760          atol,rtol,               &
    761          rcntrl,icntrl,rstatus,istatus,ierr)
     749   CALL rosenbrock(nvar, var, tin, tout,  &
     750         atol, rtol,               &
     751         rcntrl, icntrl, rstatus, istatus, ierr)
    762752
    763753   !~~~> debug option: show no of steps
    … …  
    768758   ! IF OPTIONAL PARAMETERs are given for output they
    769759   ! are updated with the RETURN information
    770    IF (present(istatus_u))istatus_u(:) = istatus(:)
    771    IF (present(rstatus_u))rstatus_u(:) = rstatus(:)
    772    IF (present(ierr_u))   ierr_u       = ierr
     760   IF (PRESENT(istatus_u))istatus_u(:) = istatus(:)
     761   IF (PRESENT(rstatus_u))rstatus_u(:) = rstatus(:)
     762   IF (PRESENT(ierr_u))  ierr_u       = ierr
    773763
    774764END SUBROUTINE integrate
    775765 
    776 SUBROUTINE fun(v,f,rct,vdot)
     766SUBROUTINE fun(v, f, rct, vdot)
    777767
    778768! V - Concentrations of variable species (local)
    … …  
    794784END SUBROUTINE fun
    795785 
    796 SUBROUTINE kppsolve(jvs,x)
     786SUBROUTINE kppsolve(jvs, x)
    797787
    798788! JVS - sparse Jacobian of variables
    … …  
    801791  REAL(kind=dp):: x(nvar)
    802792
    803   x(2) = x(2)/ jvs(2)
    804   x(1) = x(1)/ jvs(1)
     793  x(2) = x(2) / jvs(2)
     794  x(1) = x(1) / jvs(1)
    805795     
    806796END SUBROUTINE kppsolve
    807797 
    808 SUBROUTINE kppdecomp( jvs,ier)
    809 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    810 !        Sparse LU factorization
    811 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    812 
    813 
    814       INTEGER  :: ier
    815       REAL(kind=dp):: jvs(lu_nonzero),w(nvar),a
    816       INTEGER  :: k,kk,j,jj
    817 
    818       a = 0. ! mz_rs_20050606
    819       ier = 0
    820       DO k=1,nvar
    821         ! mz_rs_20050606: don't check if real value == 0
    822         ! IF(jvs( lu_diag(k)).eq. 0.)THEN
    823         IF(abs(jvs(lu_diag(k)))< tiny(a))THEN
    824             ier = k
    825             RETURN
    826         ENDIF
    827         DO kk = lu_crow(k),lu_crow(k+ 1)- 1
    828               w( lu_icol(kk)) = jvs(kk)
    829         ENDDO
    830         DO kk = lu_crow(k),lu_diag(k)- 1
    831             j = lu_icol(kk)
    832             a = - w(j)/ jvs( lu_diag(j))
    833             w(j) = - a
    834             DO jj = lu_diag(j)+ 1,lu_crow(j+ 1)- 1
    835                w( lu_icol(jj)) = w( lu_icol(jj))+ a*jvs(jj)
    836             ENDDO
    837          ENDDO
    838          DO kk = lu_crow(k),lu_crow(k+ 1)- 1
    839             jvs(kk) = w( lu_icol(kk))
    840          ENDDO
    841       ENDDO
    842      
    843 END SUBROUTINE kppdecomp
    844  
    845       REAL(kind=dp)FUNCTION wlamch( c)
    846 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    847 !     returns epsilon machine
    848 !     after LAPACK
    849 !     replace this by the function from the optimized LAPACK implementation:
    850 !          CALL SLAMCH('E') or CALL DLAMCH('E')
    851 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    852 !      USE chem_gasphase_mod_Precision
    853 
    854       CHARACTER ::  c
    855       INTEGER    :: i
    856       REAL(kind=dp),SAVE  ::  eps
    857       REAL(kind=dp) ::  suma
    858       REAL(kind=dp),PARAMETER  ::  one=1.0_dp, half=0.5_dp
    859       LOGICAL,SAVE   ::  first=.true.
    860      
    861       IF (first)THEN
    862         first = .false.
    863         eps = half**(16)
    864         DO i = 17,80
    865           eps = eps*half
    866           CALL wlamch_add(one,eps,suma)
    867           IF (suma.le.one)goto 10
    868         ENDDO
    869         PRINT*,'ERROR IN WLAMCH. EPS < ',Eps
    870         RETURN
    871 10      eps = eps*2
    872         i = i- 1     
    873       ENDIF
    874 
    875       wlamch = eps
    876 
    877       END FUNCTION wlamch
    878  
    879       SUBROUTINE wlamch_add( a,b,suma)
    880 !      USE chem_gasphase_mod_Precision
    881      
    882       REAL(kind=dp)a,b,suma
    883       suma = a + b
    884 
    885       END SUBROUTINE wlamch_add
    886  
    887 SUBROUTINE jac_sp(v,f,rct,jvs)
     798SUBROUTINE jac_sp(v, f, rct, jvs)
    888799
    889800! V - Concentrations of variable species (local)
    … …  
    914825END SUBROUTINE jac_sp
    915826 
    916   elemental REAL(kind=dp)FUNCTION k_arr (k_298,tdep,temp)
     827  elemental REAL(kind=dp)FUNCTION k_arr (k_298, tdep, temp)
    917828    ! arrhenius FUNCTION
    918829
    919     REAL,    INTENT(in):: k_298 ! k at t = 298.15k
    920     REAL,    INTENT(in):: tdep  ! temperature dependence
    921     REAL(kind=dp),INTENT(in):: temp  ! temperature
     830    REAL,    INTENT(IN):: k_298 ! k at t = 298.15k
     831    REAL,    INTENT(IN):: tdep  ! temperature dependence
     832    REAL(kind=dp), INTENT(IN):: temp  ! temperature
    922833
    923834    intrinsic exp
    924835
    925     k_arr = k_298 *exp(tdep*(1._dp/temp- 3.3540e-3_dp))! 1/298.15=3.3540e-3
     836    k_arr = k_298 * exp(tdep* (1._dp/temp- 3.3540e-3_dp))! 1/298.15=3.3540e-3
    926837
    927838  END FUNCTION k_arr
    928839 
    929840SUBROUTINE update_rconst()
    930  INTEGER         :: j,k
     841 INTEGER         :: k
    931842
    932843  k = is
    … …  
    942853END SUBROUTINE update_rconst
    943854 
    944 REAL(kind=dp)FUNCTION arr2( a0,b0,temp)
     855!  END FUNCTION ARR2
     856REAL(kind=dp)FUNCTION arr2( a0, b0, temp)
    945857    REAL(kind=dp):: temp
    946     REAL(kind=dp):: a0,b0
    947     arr2 = a0 *exp( - b0 / temp)
     858    REAL(kind=dp):: a0, b0
     859    arr2 = a0 * exp( - b0 / temp)
    948860END FUNCTION arr2
    949861 
    950 SUBROUTINE initialize_kpp_ctrl(status,iou,modstr)
     862SUBROUTINE initialize_kpp_ctrl(status)
    951863
    952864
    953865  ! i/o
    954   INTEGER,         INTENT(out):: status
    955   INTEGER,         INTENT(in) :: iou     ! LOGICAL i/o unit
    956   CHARACTER(len=*),INTENT(in) :: modstr  ! read <modstr>.nml
     866  INTEGER,         INTENT(OUT):: status
    957867
    958868  ! local
    … …  
    962872  ! check fixed time steps
    963873  tsum = 0.0_dp
    964   DO i=1,nmaxfixsteps
     874  DO i=1, nmaxfixsteps
    965875     IF (t_steps(i)< tiny(0.0_dp))exit
    966876     tsum = tsum + t_steps(i)
    … …  
    982892  WRITE(*,*) ' RCNTRL           : ',rcntrl
    983893  !
    984   ! note: this is only meaningful for vectorized (kp4)rosenbrock- methods
     894  ! note: this is ONLY meaningful for vectorized (kp4)rosenbrock- methods
    985895  IF (l_vector)THEN
    986896     IF (l_fixed_step)THEN
    … …  
    1000910END SUBROUTINE initialize_kpp_ctrl
    1001911 
    1002 SUBROUTINE error_output(c,ierr,pe)
    1003 
    1004 
    1005   INTEGER,INTENT(in):: ierr
    1006   INTEGER,INTENT(in):: pe
    1007   REAL(dp),DIMENSION(:),INTENT(in):: c
     912SUBROUTINE error_output(c, ierr, pe)
     913
     914
     915  INTEGER, INTENT(IN):: ierr
     916  INTEGER, INTENT(IN):: pe
     917  REAL(dp), DIMENSION(:), INTENT(IN):: c
    1008918
    1009919  write(6,*) 'ERROR in chem_gasphase_mod ',ierr,C(1)
    … …  
    1012922END SUBROUTINE error_output
    1013923 
    1014       SUBROUTINE wcopy(n,x,incx,y,incy)
    1015 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1016 !     copies a vector,x,to a vector,y:  y <- x
    1017 !     only for incX=incY=1
    1018 !     after BLAS
    1019 !     replace this by the function from the optimized BLAS implementation:
    1020 !         CALL  SCOPY(N,X,1,Y,1)   or   CALL  DCOPY(N,X,1,Y,1)
    1021 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1022 !     USE chem_gasphase_mod_Precision
    1023      
    1024       INTEGER  :: i,incx,incy,m,mp1,n
    1025       REAL(kind=dp):: x(n),y(n)
    1026 
    1027       IF (n.le.0)RETURN
    1028 
    1029       m = mod(n,8)
    1030       IF( m .ne. 0)THEN
    1031         DO i = 1,m
    1032           y(i) = x(i)
    1033         ENDDO
    1034         IF( n .lt. 8)RETURN
    1035       ENDIF   
    1036       mp1 = m+ 1
    1037       DO i = mp1,n,8
    1038         y(i) = x(i)
    1039         y(i + 1) = x(i + 1)
    1040         y(i + 2) = x(i + 2)
    1041         y(i + 3) = x(i + 3)
    1042         y(i + 4) = x(i + 4)
    1043         y(i + 5) = x(i + 5)
    1044         y(i + 6) = x(i + 6)
    1045         y(i + 7) = x(i + 7)
    1046       ENDDO
    1047 
    1048       END SUBROUTINE wcopy
    1049  
    1050       SUBROUTINE wscal(n,alpha,x,incx)
     924      SUBROUTINE wscal(n, alpha, x, incx)
    1051925!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1052926!     constant times a vector: x(1:N) <- Alpha*x(1:N)
    … …  
    1057931!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1058932
    1059       INTEGER  :: i,incx,m,mp1,n
    1060       REAL(kind=dp) :: x(n),alpha
    1061       REAL(kind=dp),PARAMETER  :: zero=0.0_dp, one=1.0_dp
     933      INTEGER  :: i, incx, m, mp1, n
     934      REAL(kind=dp) :: x(n), alpha
     935      REAL(kind=dp), PARAMETER  :: zero=0.0_dp, one=1.0_dp
    1062936
    1063937      IF (alpha .eq. one)RETURN
    1064938      IF (n .le. 0)RETURN
    1065939
    1066       m = mod(n,5)
    1067       IF( m .ne. 0)THEN
     940      m = mod(n, 5)
     941      IF ( m .ne. 0)THEN
    1068942        IF (alpha .eq. (- one))THEN
    1069           DO i = 1,m
     943          DO i = 1, m
    1070944            x(i) = - x(i)
    1071945          ENDDO
    1072946        ELSEIF (alpha .eq. zero)THEN
    1073           DO i = 1,m
     947          DO i = 1, m
    1074948            x(i) = zero
    1075949          ENDDO
    1076950        ELSE
    1077           DO i = 1,m
    1078             x(i) = alpha*x(i)
     951          DO i = 1, m
     952            x(i) = alpha* x(i)
    1079953          ENDDO
    1080954        ENDIF
    1081         IF( n .lt. 5)RETURN
     955        IF ( n .lt. 5)RETURN
    1082956      ENDIF
    1083957      mp1 = m + 1
    1084958      IF (alpha .eq. (- one))THEN
    1085         DO i = mp1,n,5
    1086           x(i)    = - x(i)
     959        DO i = mp1, n, 5
     960          x(i)   = - x(i)
    1087961          x(i + 1) = - x(i + 1)
    1088962          x(i + 2) = - x(i + 2)
    … …  
    1091965        ENDDO
    1092966      ELSEIF (alpha .eq. zero)THEN
    1093         DO i = mp1,n,5
    1094           x(i)    = zero
     967        DO i = mp1, n, 5
     968          x(i)   = zero
    1095969          x(i + 1) = zero
    1096970          x(i + 2) = zero
    … …  
    1099973        ENDDO
    1100974      ELSE
    1101         DO i = mp1,n,5
    1102           x(i)    = alpha*x(i)
    1103           x(i + 1) = alpha*x(i + 1)
    1104           x(i + 2) = alpha*x(i + 2)
    1105           x(i + 3) = alpha*x(i + 3)
    1106           x(i + 4) = alpha*x(i + 4)
     975        DO i = mp1, n, 5
     976          x(i)   = alpha* x(i)
     977          x(i + 1) = alpha* x(i + 1)
     978          x(i + 2) = alpha* x(i + 2)
     979          x(i + 3) = alpha* x(i + 3)
     980          x(i + 4) = alpha* x(i + 4)
    1107981        ENDDO
    1108982      ENDIF
    … …  
    1110984      END SUBROUTINE wscal
    1111985 
    1112       SUBROUTINE waxpy(n,alpha,x,incx,y,incy)
     986      SUBROUTINE waxpy(n, alpha, x, incx, y, incy)
    1113987!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1114988!     constant times a vector plus a vector: y <- y + Alpha*x
    … …  
    1119993!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    1120994
    1121       INTEGER  :: i,incx,incy,m,mp1,n
    1122       REAL(kind=dp):: x(n),y(n),alpha
    1123       REAL(kind=dp),PARAMETER :: zero = 0.0_dp
     995      INTEGER  :: i, incx, incy, m, mp1, n
     996      REAL(kind=dp):: x(n), y(n), alpha
     997      REAL(kind=dp), PARAMETER :: zero = 0.0_dp
    1124998
    1125999      IF (alpha .eq. zero)RETURN
    11261000      IF (n .le. 0)RETURN
    11271001
    1128       m = mod(n,4)
    1129       IF( m .ne. 0)THEN
    1130         DO i = 1,m
    1131           y(i) = y(i)+ alpha*x(i)
     1002      m = mod(n, 4)
     1003      IF ( m .ne. 0)THEN
     1004        DO i = 1, m
     1005          y(i) = y(i) + alpha* x(i)
    11321006        ENDDO
    1133         IF( n .lt. 4)RETURN
     1007        IF ( n .lt. 4)RETURN
    11341008      ENDIF
    11351009      mp1 = m + 1
    1136       DO i = mp1,n,4
    1137         y(i) = y(i)+ alpha*x(i)
    1138         y(i + 1) = y(i + 1)+ alpha*x(i + 1)
    1139         y(i + 2) = y(i + 2)+ alpha*x(i + 2)
    1140         y(i + 3) = y(i + 3)+ alpha*x(i + 3)
     1010      DO i = mp1, n, 4
     1011        y(i) = y(i) + alpha* x(i)
     1012        y(i + 1) = y(i + 1) + alpha* x(i + 1)
     1013        y(i + 2) = y(i + 2) + alpha* x(i + 2)
     1014        y(i + 3) = y(i + 3) + alpha* x(i + 3)
    11411015      ENDDO
    11421016     
    11431017      END SUBROUTINE waxpy
    11441018 
    1145 SUBROUTINE rosenbrock(n,y,tstart,tend,&
    1146            abstol,reltol,             &
    1147            rcntrl,icntrl,rstatus,istatus,ierr)
     1019SUBROUTINE rosenbrock(n, y, tstart, tend, &
     1020           abstol, reltol,             &
     1021           rcntrl, icntrl, rstatus, istatus, ierr)
    11481022!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    11491023!
    … …  
    11721046!~~~>   input arguments:
    11731047!
    1174 !-      y(n)   = vector of initial conditions (at t=tstart)
    1175 !-     [tstart,tend]  = time range of integration
     1048!-     y(n)  = vector of initial conditions (at t=tstart)
     1049!-    [tstart, tend]  = time range of integration
    11761050!     (if Tstart>Tend the integration is performed backwards in time)
    1177 !-     reltol,abstol = user precribed accuracy
    1178 !-  SUBROUTINE  fun( t,y,ydot) = ode FUNCTION,
     1051!-    reltol, abstol = user precribed accuracy
     1052!- SUBROUTINE  fun( t, y, ydot) = ode FUNCTION,
    11791053!                       returns Ydot = Y' = F(T,Y)
    1180 !-  SUBROUTINE  jac( t,y,jcb) = jacobian of the ode FUNCTION,
     1054!- SUBROUTINE  jac( t, y, jcb) = jacobian of the ode FUNCTION,
    11811055!                       returns Jcb = dFun/dY
    1182 !-     icntrl(1:20)   = INTEGER inputs PARAMETERs
    1183 !-     rcntrl(1:20)   = REAL inputs PARAMETERs
     1056!-    icntrl(1:20)  = INTEGER inputs PARAMETERs
     1057!-    rcntrl(1:20)  = REAL inputs PARAMETERs
    11841058!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    11851059!
    11861060!~~~>     output arguments:
    11871061!
    1188 !-     y(n)   - > vector of final states (at t- >tEND)
    1189 !-     istatus(1:20) - > INTEGER output PARAMETERs
    1190 !-     rstatus(1:20) - > REAL output PARAMETERs
    1191 !-     ierr            - > job status upon RETURN
     1062!-    y(n)  - > vector of final states (at t- >tend)
     1063!-    istatus(1:20) - > INTEGER output PARAMETERs
     1064!-    rstatus(1:20) - > REAL output PARAMETERs
     1065!-    ierr            - > job status upon RETURN
    11921066!                        success (positive value) or
    11931067!                        failure (negative value)
    … …  
    12621136
    12631137!~~~>  arguments
    1264    INTEGER,      INTENT(in)   :: n
    1265    REAL(kind=dp),INTENT(inout):: y(n)
    1266    REAL(kind=dp),INTENT(in)   :: tstart,tend
    1267    REAL(kind=dp),INTENT(in)   :: abstol(n),reltol(n)
    1268    INTEGER,      INTENT(in)   :: icntrl(20)
    1269    REAL(kind=dp),INTENT(in)   :: rcntrl(20)
    1270    INTEGER,      INTENT(inout):: istatus(20)
    1271    REAL(kind=dp),INTENT(inout):: rstatus(20)
    1272    INTEGER,INTENT(out) :: ierr
    1273 !~~~>  PARAMETERs of the rosenbrock method,up to 6 stages
    1274    INTEGER ::  ros_s,rosmethod
    1275    INTEGER,PARAMETER :: rs2=1,rs3=2,rs4=3,rd3=4,rd4=5,rg3=6
    1276    REAL(kind=dp):: ros_a(15),ros_c(15),ros_m(6),ros_e(6),&
    1277                     ros_alpha(6),ros_gamma(6),ros_elo
     1138   INTEGER,      INTENT(IN)  :: n
     1139   REAL(kind=dp), INTENT(INOUT):: y(n)
     1140   REAL(kind=dp), INTENT(IN)  :: tstart, tend
     1141   REAL(kind=dp), INTENT(IN)  :: abstol(n), reltol(n)
     1142   INTEGER,      INTENT(IN)  :: icntrl(20)
     1143   REAL(kind=dp), INTENT(IN)  :: rcntrl(20)
     1144   INTEGER,      INTENT(INOUT):: istatus(20)
     1145   REAL(kind=dp), INTENT(INOUT):: rstatus(20)
     1146   INTEGER, INTENT(OUT) :: ierr
     1147!~~~>  PARAMETERs of the rosenbrock method, up to 6 stages
     1148   INTEGER ::  ros_s, rosmethod
     1149   INTEGER, PARAMETER :: rs2=1, rs3=2, rs4=3, rd3=4, rd4=5, rg3=6
     1150   REAL(kind=dp):: ros_a(15), ros_c(15), ros_m(6), ros_e(6), &
     1151                    ros_alpha(6), ros_gamma(6), ros_elo
    12781152   LOGICAL :: ros_newf(6)
    12791153   CHARACTER(len=12):: ros_name
    12801154!~~~>  local variables
    1281    REAL(kind=dp):: roundoff,facmin,facmax,facrej,facsafe
    1282    REAL(kind=dp):: hmin,hmax,hstart
     1155   REAL(kind=dp):: roundoff, facmin, facmax, facrej, facsafe
     1156   REAL(kind=dp):: hmin, hmax, hstart
    12831157   REAL(kind=dp):: texit
    1284    INTEGER       :: i,uplimtol,max_no_steps
    1285    LOGICAL       :: autonomous,vectortol
     1158   INTEGER       :: i, uplimtol, max_no_steps
     1159   LOGICAL       :: autonomous, vectortol
    12861160!~~~>   PARAMETERs
    1287    REAL(kind=dp),PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
    1288    REAL(kind=dp),PARAMETER :: deltamin = 1.0e-5_dp
     1161   REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
     1162   REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
    12891163
    12901164!~~~>  initialize statistics
    … …  
    12981172!   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
    12991173   IF (icntrl(2) == 0)THEN
    1300       vectortol = .true.
     1174      vectortol = .TRUE.
    13011175      uplimtol  = n
    13021176   ELSE
    1303       vectortol = .false.
     1177      vectortol = .FALSE.
    13041178      uplimtol  = 1
    13051179   ENDIF
    … …  
    13131187     CASE (3)
    13141188       CALL ros4
    1315      CASE (0,4)
     1189     CASE (0, 4)
    13161190       CALL rodas3
    13171191     CASE (5)
    … …  
    13211195     CASE default
    13221196       PRINT *,'Unknown Rosenbrock method: ICNTRL(3) =',ICNTRL(3)
    1323        CALL ros_errormsg(- 2,tstart,zero,ierr)
     1197       CALL ros_errormsg(- 2, tstart, zero, ierr)
    13241198       RETURN
    13251199   END select
    … …  
    13321206   ELSE
    13331207      PRINT *,'User-selected max no. of steps: ICNTRL(4) =',ICNTRL(4)
    1334       CALL ros_errormsg(- 1,tstart,zero,ierr)
     1208      CALL ros_errormsg(- 1, tstart, zero, ierr)
    13351209      RETURN
    13361210   ENDIF
    13371211
    13381212!~~~>  unit roundoff (1+ roundoff>1)
    1339    Roundoff = WLAMCH('E')
     1213   roundoff = epsilon(one)
    13401214
    13411215!~~~>  lower bound on the step size: (positive value)
    … …  
    13461220   ELSE
    13471221      PRINT *,'User-selected Hmin: RCNTRL(1) =',RCNTRL(1)
    1348       CALL ros_errormsg(- 3,tstart,zero,ierr)
     1222      CALL ros_errormsg(- 3, tstart, zero, ierr)
    13491223      RETURN
    13501224   ENDIF
    … …  
    13531227      hmax = abs(tend-tstart)
    13541228   ELSEIF (rcntrl(2)> zero)THEN
    1355       hmax = min(abs(rcntrl(2)),abs(tend-tstart))
     1229      hmax = min(abs(rcntrl(2)), abs(tend-tstart))
    13561230   ELSE
    13571231      PRINT *,'User-selected Hmax: RCNTRL(2) =',RCNTRL(2)
    1358       CALL ros_errormsg(- 3,tstart,zero,ierr)
     1232      CALL ros_errormsg(- 3, tstart, zero, ierr)
    13591233      RETURN
    13601234   ENDIF
    13611235!~~~>  starting step size: (positive value)
    13621236   IF (rcntrl(3) == zero)THEN
    1363       hstart = max(hmin,deltamin)
     1237      hstart = max(hmin, deltamin)
    13641238   ELSEIF (rcntrl(3)> zero)THEN
    1365       hstart = min(abs(rcntrl(3)),abs(tend-tstart))
     1239      hstart = min(abs(rcntrl(3)), abs(tend-tstart))
    13661240   ELSE
    13671241      PRINT *,'User-selected Hstart: RCNTRL(3) =',RCNTRL(3)
    1368       CALL ros_errormsg(- 3,tstart,zero,ierr)
     1242      CALL ros_errormsg(- 3, tstart, zero, ierr)
    13691243      RETURN
    13701244   ENDIF
    … …  
    13761250   ELSE
    13771251      PRINT *,'User-selected FacMin: RCNTRL(4) =',RCNTRL(4)
    1378       CALL ros_errormsg(- 4,tstart,zero,ierr)
     1252      CALL ros_errormsg(- 4, tstart, zero, ierr)
    13791253      RETURN
    13801254   ENDIF
    … …  
    13851259   ELSE
    13861260      PRINT *,'User-selected FacMax: RCNTRL(5) =',RCNTRL(5)
    1387       CALL ros_errormsg(- 4,tstart,zero,ierr)
     1261      CALL ros_errormsg(- 4, tstart, zero, ierr)
    13881262      RETURN
    13891263   ENDIF
    … …  
    13951269   ELSE
    13961270      PRINT *,'User-selected FacRej: RCNTRL(6) =',RCNTRL(6)
    1397       CALL ros_errormsg(- 4,tstart,zero,ierr)
     1271      CALL ros_errormsg(- 4, tstart, zero, ierr)
    13981272      RETURN
    13991273   ENDIF
    … …  
    14051279   ELSE
    14061280      PRINT *,'User-selected FacSafe: RCNTRL(7) =',RCNTRL(7)
    1407       CALL ros_errormsg(- 4,tstart,zero,ierr)
     1281      CALL ros_errormsg(- 4, tstart, zero, ierr)
    14081282      RETURN
    14091283   ENDIF
    14101284!~~~>  check IF tolerances are reasonable
    1411     DO i=1,uplimtol
    1412       IF((abstol(i)<= zero).or. (reltol(i)<= 10.0_dp*roundoff)&
     1285    DO i=1, uplimtol
     1286      IF ((abstol(i)<= zero).or. (reltol(i)<= 10.0_dp* roundoff)&
    14131287         .or. (reltol(i)>= 1.0_dp))THEN
    14141288        PRINT *,' AbsTol(',i,') = ',AbsTol(i)
    14151289        PRINT *,' RelTol(',i,') = ',RelTol(i)
    1416         CALL ros_errormsg(- 5,tstart,zero,ierr)
     1290        CALL ros_errormsg(- 5, tstart, zero, ierr)
    14171291        RETURN
    14181292      ENDIF
    … …  
    14211295
    14221296!~~~>  CALL rosenbrock method
    1423    CALL ros_integrator(y,tstart,tend,texit,  &
    1424         abstol,reltol,                         &
     1297   CALL ros_integrator(y, tstart, tend, texit,  &
     1298        abstol, reltol,                         &
    14251299!  Integration parameters
    1426         autonomous,vectortol,max_no_steps,    &
    1427         roundoff,hmin,hmax,hstart,           &
    1428         facmin,facmax,facrej,facsafe,        &
     1300        autonomous, vectortol, max_no_steps,    &
     1301        roundoff, hmin, hmax, hstart,           &
     1302        facmin, facmax, facrej, facsafe,        &
    14291303!  Error indicator
    14301304        ierr)
    … …  
    14351309   
    14361310!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
    1437  SUBROUTINE ros_errormsg(code,t,h,ierr)
     1311 SUBROUTINE ros_errormsg(code, t, h, ierr)
    14381312!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
    14391313!    Handles all error messages
    14401314!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
    14411315 
    1442    REAL(kind=dp),INTENT(in):: t,h
    1443    INTEGER,INTENT(in) :: code
    1444    INTEGER,INTENT(out):: ierr
     1316   REAL(kind=dp), INTENT(IN):: t, h
     1317   INTEGER, INTENT(IN) :: code
     1318   INTEGER, INTENT(OUT):: ierr
    14451319   
    14461320   ierr = code
    1447    print *,&
     1321   print * , &
    14481322     'Forced exit from Rosenbrock due to the following error:'
    14491323     
    14501324   select CASE (code)
    1451     CASE (- 1)   
     1325    CASE (- 1) 
    14521326      PRINT *,'--> Improper value for maximal no of steps'
    1453     CASE (- 2)   
     1327    CASE (- 2) 
    14541328      PRINT *,'--> Selected Rosenbrock method not implemented'
    1455     CASE (- 3)   
     1329    CASE (- 3) 
    14561330      PRINT *,'--> Hmin/Hmax/Hstart must be positive'
    1457     CASE (- 4)   
     1331    CASE (- 4) 
    14581332      PRINT *,'--> FacMin/FacMax/FacRej must be positive'
    14591333    CASE (- 5)
    … …  
    14641338      PRINT *,'--> Step size too small: T + 10*H = T',&
    14651339            ' or H < Roundoff'
    1466     CASE (- 8)   
     1340    CASE (- 8) 
    14671341      PRINT *,'--> Matrix is repeatedly singular'
    14681342    CASE default
    … …  
    14701344   END select
    14711345   
    1472    print *,"t=",t,"and h=",h
     1346   print * , "t=", t, "and h=", h
    14731347     
    14741348 END SUBROUTINE ros_errormsg
    14751349
    14761350!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1477  SUBROUTINE ros_integrator (y,tstart,tend,t, &
    1478         abstol,reltol,                         &
     1351 SUBROUTINE ros_integrator (y, tstart, tend, t, &
     1352        abstol, reltol,                         &
    14791353!~~~> integration PARAMETERs
    1480         autonomous,vectortol,max_no_steps,    &
    1481         roundoff,hmin,hmax,hstart,           &
    1482         facmin,facmax,facrej,facsafe,        &
     1354        autonomous, vectortol, max_no_steps,    &
     1355        roundoff, hmin, hmax, hstart,           &
     1356        facmin, facmax, facrej, facsafe,        &
    14831357!~~~> error indicator
    14841358        ierr)
    … …  
    14911365
    14921366!~~~> input: the initial condition at tstart; output: the solution at t
    1493    REAL(kind=dp),INTENT(inout):: y(n)
     1367   REAL(kind=dp), INTENT(INOUT):: y(n)
    14941368!~~~> input: integration interval
    1495    REAL(kind=dp),INTENT(in):: tstart,tend
    1496 !~~~> output: time at which the solution is RETURNed (t=tENDIF success)
    1497    REAL(kind=dp),INTENT(out)::  t
     1369   REAL(kind=dp), INTENT(IN):: tstart, tend
     1370!~~~> output: time at which the solution is RETURNed (t=tendIF success)
     1371   REAL(kind=dp), INTENT(OUT)::  t
    14981372!~~~> input: tolerances
    1499    REAL(kind=dp),INTENT(in)::  abstol(n),reltol(n)
     1373   REAL(kind=dp), INTENT(IN)::  abstol(n), reltol(n)
    15001374!~~~> input: integration PARAMETERs
    1501    LOGICAL,INTENT(in):: autonomous,vectortol
    1502    REAL(kind=dp),INTENT(in):: hstart,hmin,hmax
    1503    INTEGER,INTENT(in):: max_no_steps
    1504    REAL(kind=dp),INTENT(in):: roundoff,facmin,facmax,facrej,facsafe
     1375   LOGICAL, INTENT(IN):: autonomous, vectortol
     1376   REAL(kind=dp), INTENT(IN):: hstart, hmin, hmax
     1377   INTEGER, INTENT(IN):: max_no_steps
     1378   REAL(kind=dp), INTENT(IN):: roundoff, facmin, facmax, facrej, facsafe
    15051379!~~~> output: error indicator
    1506    INTEGER,INTENT(out):: ierr
     1380   INTEGER, INTENT(OUT):: ierr
    15071381! ~~~~ Local variables
    1508    REAL(kind=dp):: ynew(n),fcn0(n),fcn(n)
    1509    REAL(kind=dp):: k(n*ros_s),dfdt(n)
     1382   REAL(kind=dp):: ynew(n), fcn0(n), fcn(n)
     1383   REAL(kind=dp):: k(n* ros_s), dfdt(n)
    15101384#ifdef full_algebra   
    1511    REAL(kind=dp):: jac0(n,n),ghimj(n,n)
     1385   REAL(kind=dp):: jac0(n, n), ghimj(n, n)
    15121386#else
    1513    REAL(kind=dp):: jac0(lu_nonzero),ghimj(lu_nonzero)
     1387   REAL(kind=dp):: jac0(lu_nonzero), ghimj(lu_nonzero)
    15141388#endif
    1515    REAL(kind=dp):: h,hnew,hc,hg,fac,tau
    1516    REAL(kind=dp):: err,yerr(n)
    1517    INTEGER :: pivot(n),direction,ioffset,j,istage
    1518    LOGICAL :: rejectlasth,rejectmoreh,singular
     1389   REAL(kind=dp):: h, hnew, hc, hg, fac, tau
     1390   REAL(kind=dp):: err, yerr(n)
     1391   INTEGER :: pivot(n), direction, ioffset, j, istage
     1392   LOGICAL :: rejectlasth, rejectmoreh, singular
    15191393!~~~>  local PARAMETERs
    1520    REAL(kind=dp),PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
    1521    REAL(kind=dp),PARAMETER :: deltamin = 1.0e-5_dp
     1394   REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
     1395   REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
    15221396!~~~>  locally called FUNCTIONs
    15231397!    REAL(kind=dp) WLAMCH
    … …  
    15291403   t = tstart
    15301404   rstatus(nhexit) = zero
    1531    h = min( max(abs(hmin),abs(hstart)),abs(hmax))
    1532    IF (abs(h)<= 10.0_dp*roundoff)h = deltamin
    1533 
    1534    IF (tEND  >=  tstart)THEN
     1405   h = min( max(abs(hmin), abs(hstart)), abs(hmax))
     1406   IF (abs(h)<= 10.0_dp* roundoff)h = deltamin
     1407
     1408   IF (tend  >=  tstart)THEN
    15351409     direction = + 1
    15361410   ELSE
    15371411     direction = - 1
    15381412   ENDIF
    1539    h = direction*h
    1540 
    1541    rejectlasth=.false.
    1542    rejectmoreh=.false.
     1413   h = direction* h
     1414
     1415   rejectlasth=.FALSE.
     1416   rejectmoreh=.FALSE.
    15431417
    15441418!~~~> time loop begins below
    15451419
    1546 timeloop: DO WHILE((direction > 0).and.((t- tEND)+ roundoff <= zero)&
    1547        .or. (direction < 0).and.((tend-t)+ roundoff <= zero))
    1548 
    1549    IF(istatus(nstp)> max_no_steps)THEN  ! too many steps
    1550       CALL ros_errormsg(- 6,t,h,ierr)
     1420timeloop: DO WHILE((direction > 0).and.((t- tend) + roundoff <= zero)&
     1421       .or. (direction < 0).and.((tend-t) + roundoff <= zero))
     1422
     1423   IF (istatus(nstp)> max_no_steps)THEN  ! too many steps
     1424      CALL ros_errormsg(- 6, t, h, ierr)
    15511425      RETURN
    15521426   ENDIF
    1553    IF(((t+ 0.1_dp*h) == t).or.(h <= roundoff))THEN  ! step size too small
    1554       CALL ros_errormsg(- 7,t,h,ierr)
     1427   IF (((t+ 0.1_dp* h) == t).or.(h <= roundoff))THEN  ! step size too small
     1428      CALL ros_errormsg(- 7, t, h, ierr)
    15551429      RETURN
    15561430   ENDIF
    15571431
    15581432!~~~>  limit h IF necessary to avoid going beyond tend
    1559    h = min(h,abs(tend-t))
     1433   h = min(h, abs(tend-t))
    15601434
    15611435!~~~>   compute the FUNCTION at current time
    1562    CALL funtemplate(t,y,fcn0)
    1563    istatus(nfun) = istatus(nfun)+ 1
     1436   CALL funtemplate(t, y, fcn0)
     1437   istatus(nfun) = istatus(nfun) + 1
    15641438
    15651439!~~~>  compute the FUNCTION derivative with respect to t
    15661440   IF (.not.autonomous)THEN
    1567       CALL ros_funtimederivative(t,roundoff,y,&
    1568                 fcn0,dfdt)
     1441      CALL ros_funtimederivative(t, roundoff, y, &
     1442                fcn0, dfdt)
    15691443   ENDIF
    15701444
    15711445!~~~>   compute the jacobian at current time
    1572    CALL jactemplate(t,y,jac0)
    1573    istatus(njac) = istatus(njac)+ 1
     1446   CALL jactemplate(t, y, jac0)
     1447   istatus(njac) = istatus(njac) + 1
    15741448
    15751449!~~~>  repeat step calculation until current step accepted
    15761450untilaccepted: do
    15771451
    1578    CALL ros_preparematrix(h,direction,ros_gamma(1),&
    1579           jac0,ghimj,pivot,singular)
     1452   CALL ros_preparematrix(h, direction, ros_gamma(1), &
     1453          jac0, ghimj, pivot, singular)
    15801454   IF (singular)THEN ! more than 5 consecutive failed decompositions
    1581        CALL ros_errormsg(- 8,t,h,ierr)
     1455       CALL ros_errormsg(- 8, t, h, ierr)
    15821456       RETURN
    15831457   ENDIF
    15841458
    15851459!~~~>   compute the stages
    1586 stage: DO istage = 1,ros_s
     1460stage: DO istage = 1, ros_s
    15871461
    15881462      ! current istage offset. current istage vector is k(ioffset+ 1:ioffset+ n)
    1589        ioffset = n*(istage-1)
     1463       ioffset = n* (istage-1)
    15901464
    15911465      ! for the 1st istage the FUNCTION has been computed previously
    1592        IF(istage == 1)THEN
    1593          !slim: CALL wcopy(n,fcn0,1,fcn,1)
    1594         fcn(1:n) = fcn0(1:n)
     1466       IF (istage == 1)THEN
     1467         !slim: CALL wcopy(n, fcn0, 1, fcn, 1)
     1468      fcn(1:n) = fcn0(1:n)
    15951469      ! istage>1 and a new FUNCTION evaluation is needed at the current istage
    15961470       ELSEIF(ros_newf(istage))THEN
    1597          !slim: CALL wcopy(n,y,1,ynew,1)
    1598         ynew(1:n) = y(1:n)
    1599          DO j = 1,istage-1
    1600            CALL waxpy(n,ros_a((istage-1)*(istage-2)/2+ j),&
    1601             k(n*(j- 1)+ 1),1,ynew,1)
     1471         !slim: CALL wcopy(n, y, 1, ynew, 1)
     1472      ynew(1:n) = y(1:n)
     1473         DO j = 1, istage-1
     1474           CALL waxpy(n, ros_a((istage-1) * (istage-2) /2+ j), &
     1475            k(n* (j- 1) + 1), 1, ynew, 1)
    16021476         ENDDO
    1603          tau = t + ros_alpha(istage)*direction*h
    1604          CALL funtemplate(tau,ynew,fcn)
    1605          istatus(nfun) = istatus(nfun)+ 1
     1477         tau = t + ros_alpha(istage) * direction* h
     1478         CALL funtemplate(tau, ynew, fcn)
     1479         istatus(nfun) = istatus(nfun) + 1
    16061480       ENDIF ! IF istage == 1 ELSEIF ros_newf(istage)
    1607        !slim: CALL wcopy(n,fcn,1,k(ioffset+ 1),1)
     1481       !slim: CALL wcopy(n, fcn, 1, k(ioffset+ 1), 1)
    16081482       k(ioffset+ 1:ioffset+ n) = fcn(1:n)
    1609        DO j = 1,istage-1
    1610          hc = ros_c((istage-1)*(istage-2)/2+ j)/(direction*h)
    1611          CALL waxpy(n,hc,k(n*(j- 1)+ 1),1,k(ioffset+ 1),1)
     1483       DO j = 1, istage-1
     1484         hc = ros_c((istage-1) * (istage-2) /2+ j) /(direction* h)
     1485         CALL waxpy(n, hc, k(n* (j- 1) + 1), 1, k(ioffset+ 1), 1)
    16121486       ENDDO
    16131487       IF ((.not. autonomous).and.(ros_gamma(istage).ne.zero))THEN
    1614          hg = direction*h*ros_gamma(istage)
    1615          CALL waxpy(n,hg,dfdt,1,k(ioffset+ 1),1)
     1488         hg = direction* h* ros_gamma(istage)
     1489         CALL waxpy(n, hg, dfdt, 1, k(ioffset+ 1), 1)
    16161490       ENDIF
    1617        CALL ros_solve(ghimj,pivot,k(ioffset+ 1))
     1491       CALL ros_solve(ghimj, pivot, k(ioffset+ 1))
    16181492
    16191493   END DO stage
    … …  
    16211495
    16221496!~~~>  compute the new solution
    1623    !slim: CALL wcopy(n,y,1,ynew,1)
     1497   !slim: CALL wcopy(n, y, 1, ynew, 1)
    16241498   ynew(1:n) = y(1:n)
    1625    DO j=1,ros_s
    1626          CALL waxpy(n,ros_m(j),k(n*(j- 1)+ 1),1,ynew,1)
     1499   DO j=1, ros_s
     1500         CALL waxpy(n, ros_m(j), k(n* (j- 1) + 1), 1, ynew, 1)
    16271501   ENDDO
    16281502
    16291503!~~~>  compute the error estimation
    1630    !slim: CALL wscal(n,zero,yerr,1)
     1504   !slim: CALL wscal(n, zero, yerr, 1)
    16311505   yerr(1:n) = zero
    1632    DO j=1,ros_s
    1633         CALL waxpy(n,ros_e(j),k(n*(j- 1)+ 1),1,yerr,1)
     1506   DO j=1, ros_s
     1507        CALL waxpy(n, ros_e(j), k(n* (j- 1) + 1), 1, yerr, 1)
    16341508   ENDDO
    1635    err = ros_errornorm(y,ynew,yerr,abstol,reltol,vectortol)
     1509   err = ros_errornorm(y, ynew, yerr, abstol, reltol, vectortol)
    16361510
    16371511!~~~> new step size is bounded by facmin <= hnew/h <= facmax
    1638    fac  = min(facmax,max(facmin,facsafe/err**(one/ros_elo)))
    1639    hnew = h*fac
     1512   fac  = min(facmax, max(facmin, facsafe/err** (one/ros_elo)))
     1513   hnew = h* fac
    16401514
    16411515!~~~>  check the error magnitude and adjust step size
    1642    istatus(nstp) = istatus(nstp)+ 1
    1643    IF((err <= one).or.(h <= hmin))THEN  !~~~> accept step
    1644       istatus(nacc) = istatus(nacc)+ 1
    1645       !slim: CALL wcopy(n,ynew,1,y,1)
     1516   istatus(nstp) = istatus(nstp) + 1
     1517   IF ((err <= one).or.(h <= hmin))THEN  !~~~> accept step
     1518      istatus(nacc) = istatus(nacc) + 1
     1519      !slim: CALL wcopy(n, ynew, 1, y, 1)
    16461520      y(1:n) = ynew(1:n)
    1647       t = t + direction*h
    1648       hnew = max(hmin,min(hnew,hmax))
     1521      t = t + direction* h
     1522      hnew = max(hmin, min(hnew, hmax))
    16491523      IF (rejectlasth)THEN  ! no step size increase after a rejected step
    1650          hnew = min(hnew,h)
     1524         hnew = min(hnew, h)
    16511525      ENDIF
    16521526      rstatus(nhexit) = h
    16531527      rstatus(nhnew) = hnew
    16541528      rstatus(ntexit) = t
    1655       rejectlasth = .false.
    1656       rejectmoreh = .false.
     1529      rejectlasth = .FALSE.
     1530      rejectmoreh = .FALSE.
    16571531      h = hnew
    16581532      exit untilaccepted ! exit the loop: WHILE step not accepted
    16591533   ELSE           !~~~> reject step
    16601534      IF (rejectmoreh)THEN
    1661          hnew = h*facrej
     1535         hnew = h* facrej
    16621536      ENDIF
    16631537      rejectmoreh = rejectlasth
    1664       rejectlasth = .true.
     1538      rejectlasth = .TRUE.
    16651539      h = hnew
    1666       IF (istatus(nacc)>= 1) istatus(nrej) = istatus(nrej)+ 1
     1540      IF (istatus(nacc)>= 1) istatus(nrej) = istatus(nrej) + 1
    16671541   ENDIF ! err <= 1
    16681542
    … …  
    16781552
    16791553!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1680   REAL(kind=dp)FUNCTION ros_errornorm(y,ynew,yerr,&
    1681                                abstol,reltol,vectortol)
     1554  REAL(kind=dp)FUNCTION ros_errornorm(y, ynew, yerr, &
     1555                               abstol, reltol, vectortol)
    16821556!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    16831557!~~~> computes the "scaled norm" of the error vector yerr
    … …  
    16851559
    16861560! Input arguments
    1687    REAL(kind=dp),INTENT(in):: y(n),ynew(n),&
    1688           yerr(n),abstol(n),reltol(n)
    1689    LOGICAL,INTENT(in)::  vectortol
     1561   REAL(kind=dp), INTENT(IN):: y(n), ynew(n), &
     1562          yerr(n), abstol(n), reltol(n)
     1563   LOGICAL, INTENT(IN)::  vectortol
    16901564! Local variables
    1691    REAL(kind=dp):: err,scale,ymax
     1565   REAL(kind=dp):: err, scale, ymax
    16921566   INTEGER  :: i
    1693    REAL(kind=dp),PARAMETER :: zero = 0.0_dp
     1567   REAL(kind=dp), PARAMETER :: zero = 0.0_dp
    16941568
    16951569   err = zero
    1696    DO i=1,n
    1697      ymax = max(abs(y(i)),abs(ynew(i)))
     1570   DO i=1, n
     1571     ymax = max(abs(y(i)), abs(ynew(i)))
    16981572     IF (vectortol)THEN
    1699        scale = abstol(i)+ reltol(i)*ymax
     1573       scale = abstol(i) + reltol(i) * ymax
    17001574     ELSE
    1701        scale = abstol(1)+ reltol(1)*ymax
     1575       scale = abstol(1) + reltol(1) * ymax
    17021576     ENDIF
    1703      err = err+(yerr(i)/scale)**2
     1577     err = err+ (yerr(i) /scale) ** 2
    17041578   ENDDO
    17051579   err  = sqrt(err/n)
    17061580
    1707    ros_errornorm = max(err,1.0d-10)
     1581   ros_errornorm = max(err, 1.0d-10)
    17081582
    17091583  END FUNCTION ros_errornorm
    … …  
    17111585
    17121586!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1713   SUBROUTINE ros_funtimederivative(t,roundoff,y,&
    1714                 fcn0,dfdt)
     1587  SUBROUTINE ros_funtimederivative(t, roundoff, y, &
     1588                fcn0, dfdt)
    17151589!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    17161590!~~~> the time partial derivative of the FUNCTION by finite differences
    … …  
    17181592
    17191593!~~~> input arguments
    1720    REAL(kind=dp),INTENT(in):: t,roundoff,y(n),fcn0(n)
     1594   REAL(kind=dp), INTENT(IN):: t, roundoff, y(n), fcn0(n)
    17211595!~~~> output arguments
    1722    REAL(kind=dp),INTENT(out):: dfdt(n)
     1596   REAL(kind=dp), INTENT(OUT):: dfdt(n)
    17231597!~~~> local variables
    17241598   REAL(kind=dp):: delta
    1725    REAL(kind=dp),PARAMETER :: one = 1.0_dp, deltamin = 1.0e-6_dp
    1726 
    1727    delta = sqrt(roundoff)*max(deltamin,abs(t))
    1728    CALL funtemplate(t+ delta,y,dfdt)
    1729    istatus(nfun) = istatus(nfun)+ 1
    1730    CALL waxpy(n,(- one),fcn0,1,dfdt,1)
    1731    CALL wscal(n,(one/delta),dfdt,1)
     1599   REAL(kind=dp), PARAMETER :: one = 1.0_dp, deltamin = 1.0e-6_dp
     1600
     1601   delta = sqrt(roundoff) * max(deltamin, abs(t))
     1602   CALL funtemplate(t+ delta, y, dfdt)
     1603   istatus(nfun) = istatus(nfun) + 1
     1604   CALL waxpy(n, (- one), fcn0, 1, dfdt, 1)
     1605   CALL wscal(n, (one/delta), dfdt, 1)
    17321606
    17331607  END SUBROUTINE ros_funtimederivative
    … …  
    17351609
    17361610!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1737   SUBROUTINE ros_preparematrix(h,direction,gam,&
    1738              jac0,ghimj,pivot,singular)
     1611  SUBROUTINE ros_preparematrix(h, direction, gam, &
     1612             jac0, ghimj, pivot, singular)
    17391613! --- --- --- --- --- --- --- --- --- --- --- --- ---
    17401614!  Prepares the LHS matrix for stage calculations
    … …  
    17481622!~~~> input arguments
    17491623#ifdef full_algebra   
    1750    REAL(kind=dp),INTENT(in)::  jac0(n,n)
     1624   REAL(kind=dp), INTENT(IN)::  jac0(n, n)
    17511625#else
    1752    REAL(kind=dp),INTENT(in)::  jac0(lu_nonzero)
     1626   REAL(kind=dp), INTENT(IN)::  jac0(lu_nonzero)
    17531627#endif   
    1754    REAL(kind=dp),INTENT(in)::  gam
    1755    INTEGER,INTENT(in)::  direction
     1628   REAL(kind=dp), INTENT(IN)::  gam
     1629   INTEGER, INTENT(IN)::  direction
    17561630!~~~> output arguments
    17571631#ifdef full_algebra   
    1758    REAL(kind=dp),INTENT(out):: ghimj(n,n)
     1632   REAL(kind=dp), INTENT(OUT):: ghimj(n, n)
    17591633#else
    1760    REAL(kind=dp),INTENT(out):: ghimj(lu_nonzero)
     1634   REAL(kind=dp), INTENT(OUT):: ghimj(lu_nonzero)
    17611635#endif   
    1762    LOGICAL,INTENT(out)::  singular
    1763    INTEGER,INTENT(out)::  pivot(n)
     1636   LOGICAL, INTENT(OUT)::  singular
     1637   INTEGER, INTENT(OUT)::  pivot(n)
    17641638!~~~> inout arguments
    1765    REAL(kind=dp),INTENT(inout):: h   ! step size is decreased when lu fails
     1639   REAL(kind=dp), INTENT(INOUT):: h   ! step size is decreased when lu fails
    17661640!~~~> local variables
    1767    INTEGER  :: i,ising,nconsecutive
     1641   INTEGER  :: i, ising, nconsecutive
    17681642   REAL(kind=dp):: ghinv
    1769    REAL(kind=dp),PARAMETER :: one  = 1.0_dp, half = 0.5_dp
     1643   REAL(kind=dp), PARAMETER :: one  = 1.0_dp, half = 0.5_dp
    17701644
    17711645   nconsecutive = 0
    1772    singular = .true.
     1646   singular = .TRUE.
    17731647
    17741648   DO WHILE (singular)
    17751649
    1776 !~~~>    construct ghimj = 1/(h*gam)- jac0
     1650!~~~>    construct ghimj = 1/(h* gam) - jac0
    17771651#ifdef full_algebra   
    1778      !slim: CALL wcopy(n*n,jac0,1,ghimj,1)
    1779      !slim: CALL wscal(n*n,(- one),ghimj,1)
     1652     !slim: CALL wcopy(n* n, jac0, 1, ghimj, 1)
     1653     !slim: CALL wscal(n* n, (- one), ghimj, 1)
    17801654     ghimj = - jac0
    1781      ghinv = one/(direction*h*gam)
    1782      DO i=1,n
    1783        ghimj(i,i) = ghimj(i,i)+ ghinv
     1655     ghinv = one/(direction* h* gam)
     1656     DO i=1, n
     1657       ghimj(i, i) = ghimj(i, i) + ghinv
    17841658     ENDDO
    17851659#else
    1786      !slim: CALL wcopy(lu_nonzero,jac0,1,ghimj,1)
    1787      !slim: CALL wscal(lu_nonzero,(- one),ghimj,1)
     1660     !slim: CALL wcopy(lu_nonzero, jac0, 1, ghimj, 1)
     1661     !slim: CALL wscal(lu_nonzero, (- one), ghimj, 1)
    17881662     ghimj(1:lu_nonzero) = - jac0(1:lu_nonzero)
    1789      ghinv = one/(direction*h*gam)
    1790      DO i=1,n
    1791        ghimj(lu_diag(i)) = ghimj(lu_diag(i))+ ghinv
     1663     ghinv = one/(direction* h* gam)
     1664     DO i=1, n
     1665       ghimj(lu_diag(i)) = ghimj(lu_diag(i)) + ghinv
    17921666     ENDDO
    17931667#endif   
    17941668!~~~>    compute lu decomposition
    1795      CALL ros_decomp( ghimj,pivot,ising)
     1669     CALL ros_decomp( ghimj, pivot, ising)
    17961670     IF (ising == 0)THEN
    17971671!~~~>    IF successful done
    1798         singular = .false.
     1672        singular = .FALSE.
    17991673     ELSE ! ising .ne. 0
    18001674!~~~>    IF unsuccessful half the step size; IF 5 consecutive fails THEN RETURN
    1801         istatus(nsng) = istatus(nsng)+ 1
     1675        istatus(nsng) = istatus(nsng) + 1
    18021676        nconsecutive = nconsecutive+1
    1803         singular = .true.
     1677        singular = .TRUE.
    18041678        PRINT*,'Warning: LU Decomposition returned ISING = ',ISING
    18051679        IF (nconsecutive <= 5)THEN ! less than 5 consecutive failed decompositions
    1806            h = h*half
     1680           h = h* half
    18071681        ELSE  ! more than 5 consecutive failed decompositions
    18081682           RETURN
    … …  
    18161690
    18171691!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1818   SUBROUTINE ros_decomp( a,pivot,ising)
     1692  SUBROUTINE ros_decomp( a, pivot, ising)
    18191693!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    18201694!  Template for the LU decomposition
    … …  
    18221696!~~~> inout variables
    18231697#ifdef full_algebra   
    1824    REAL(kind=dp),INTENT(inout):: a(n,n)
     1698   REAL(kind=dp), INTENT(INOUT):: a(n, n)
    18251699#else   
    1826    REAL(kind=dp),INTENT(inout):: a(lu_nonzero)
     1700   REAL(kind=dp), INTENT(INOUT):: a(lu_nonzero)
    18271701#endif
    18281702!~~~> output variables
    1829    INTEGER,INTENT(out):: pivot(n),ising
     1703   INTEGER, INTENT(OUT):: pivot(n), ising
    18301704
    18311705#ifdef full_algebra   
    1832    CALL  dgetrf( n,n,a,n,pivot,ising)
     1706   CALL  dgetrf( n, n, a, n, pivot, ising)
    18331707#else   
    1834    CALL kppdecomp(a,ising)
     1708   CALL kppdecomp(a, ising)
    18351709   pivot(1) = 1
    18361710#endif
    1837    istatus(ndec) = istatus(ndec)+ 1
     1711   istatus(ndec) = istatus(ndec) + 1
    18381712
    18391713  END SUBROUTINE ros_decomp
    … …  
    18411715
    18421716!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1843   SUBROUTINE ros_solve( a,pivot,b)
     1717  SUBROUTINE ros_solve( a, pivot, b)
    18441718!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    18451719!  Template for the forward/backward substitution (using pre-computed LU decomposition)
    … …  
    18471721!~~~> input variables
    18481722#ifdef full_algebra   
    1849    REAL(kind=dp),INTENT(in):: a(n,n)
     1723   REAL(kind=dp), INTENT(IN):: a(n, n)
    18501724   INTEGER :: ising
    18511725#else   
    1852    REAL(kind=dp),INTENT(in):: a(lu_nonzero)
     1726   REAL(kind=dp), INTENT(IN):: a(lu_nonzero)
    18531727#endif
    1854    INTEGER,INTENT(in):: pivot(n)
     1728   INTEGER, INTENT(IN):: pivot(n)
    18551729!~~~> inout variables
    1856    REAL(kind=dp),INTENT(inout):: b(n)
     1730   REAL(kind=dp), INTENT(INOUT):: b(n)
    18571731
    18581732#ifdef full_algebra   
    18591733   CALL  DGETRS( 'N',N ,1,A,N,Pivot,b,N,ISING)
    1860    IF(info < 0)THEN
    1861       print*,"error in dgetrs. ising=",ising
     1734   IF (info < 0)THEN
     1735      print* , "error in dgetrs. ising=", ising
    18621736   ENDIF 
    18631737#else   
    1864    CALL kppsolve( a,b)
     1738   CALL kppsolve( a, b)
    18651739#endif
    18661740
    1867    istatus(nsol) = istatus(nsol)+ 1
     1741   istatus(nsol) = istatus(nsol) + 1
    18681742
    18691743  END SUBROUTINE ros_solve
    … …  
    18931767!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
    18941768
    1895     ros_a(1) = (1.0_dp)/g
    1896     ros_c(1) = (- 2.0_dp)/g
     1769    ros_a(1) = (1.0_dp) /g
     1770    ros_c(1) = (- 2.0_dp) /g
    18971771!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
    18981772!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
    1899     ros_newf(1) = .true.
    1900     ros_newf(2) = .true.
     1773    ros_newf(1) = .TRUE.
     1774    ros_newf(2) = .TRUE.
    19011775!~~~> m_i = coefficients for new step solution
    1902     ros_m(1) = (3.0_dp)/(2.0_dp*g)
    1903     ros_m(2) = (1.0_dp)/(2.0_dp*g)
     1776    ros_m(1) = (3.0_dp) /(2.0_dp* g)
     1777    ros_m(2) = (1.0_dp) /(2.0_dp* g)
    19041778! E_i = Coefficients for error estimator
    1905     ros_e(1) = 1.0_dp/(2.0_dp*g)
    1906     ros_e(2) = 1.0_dp/(2.0_dp*g)
    1907 !~~~> ros_elo = estimator of local order -  the minimum between the
     1779    ros_e(1) = 1.0_dp/(2.0_dp* g)
     1780    ros_e(2) = 1.0_dp/(2.0_dp* g)
     1781!~~~> ros_elo = estimator of local order - the minimum between the
    19081782!    main and the embedded scheme orders plus one
    19091783    ros_elo = 2.0_dp
    1910 !~~~> y_stage_i ~ y( t + h*alpha_i)
     1784!~~~> y_stage_i ~ y( t + h* alpha_i)
    19111785    ros_alpha(1) = 0.0_dp
    19121786    ros_alpha(2) = 1.0_dp
    1913 !~~~> gamma_i = \sum_j  gamma_{i,j}
     1787!~~~> gamma_i = \sum_j  gamma_{i, j}
    19141788    ros_gamma(1) = g
    1915     ros_gamma(2) =- g
     1789    ros_gamma(2) = -g
    19161790
    19171791 END SUBROUTINE ros2
    … …  
    19461820!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
    19471821!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
    1948    ros_newf(1) = .true.
    1949    ros_newf(2) = .true.
    1950    ros_newf(3) = .false.
     1822   ros_newf(1) = .TRUE.
     1823   ros_newf(2) = .TRUE.
     1824   ros_newf(3) = .FALSE.
    19511825!~~~> m_i = coefficients for new step solution
    19521826   ros_m(1) =  0.1e+01_dp
    … …  
    19571831   ros_e(2) = - 0.29079558716805469821718236208017e+01_dp
    19581832   ros_e(3) =  0.22354069897811569627360909276199_dp
    1959 !~~~> ros_elo = estimator of local order -  the minimum between the
     1833!~~~> ros_elo = estimator of local order - the minimum between the
    19601834!    main and the embedded scheme orders plus 1
    19611835   ros_elo = 3.0_dp
    1962 !~~~> y_stage_i ~ y( t + h*alpha_i)
     1836!~~~> y_stage_i ~ y( t + h* alpha_i)
    19631837   ros_alpha(1) = 0.0_dp
    19641838   ros_alpha(2) = 0.43586652150845899941601945119356_dp
    19651839   ros_alpha(3) = 0.43586652150845899941601945119356_dp
    1966 !~~~> gamma_i = \sum_j  gamma_{i,j}
     1840!~~~> gamma_i = \sum_j  gamma_{i, j}
    19671841   ros_gamma(1) = 0.43586652150845899941601945119356_dp
    19681842   ros_gamma(2) = 0.24291996454816804366592249683314_dp
    … …  
    20071881   ros_a(6) = 0.0_dp
    20081882
    2009    ros_c(1) =- 0.7137615036412310e+01_dp
     1883   ros_c(1) = -0.7137615036412310e+01_dp
    20101884   ros_c(2) = 0.2580708087951457e+01_dp
    20111885   ros_c(3) = 0.6515950076447975_dp
    2012    ros_c(4) =- 0.2137148994382534e+01_dp
    2013    ros_c(5) =- 0.3214669691237626_dp
    2014    ros_c(6) =- 0.6949742501781779_dp
     1886   ros_c(4) = -0.2137148994382534e+01_dp
     1887   ros_c(5) = -0.3214669691237626_dp
     1888   ros_c(6) = -0.6949742501781779_dp
    20151889!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
    20161890!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
    2017    ros_newf(1) = .true.
    2018    ros_newf(2) = .true.
    2019    ros_newf(3) = .true.
    2020    ros_newf(4) = .false.
     1891   ros_newf(1) = .TRUE.
     1892   ros_newf(2) = .TRUE.
     1893   ros_newf(3) = .TRUE.
     1894   ros_newf(4) = .FALSE.
    20211895!~~~> m_i = coefficients for new step solution
    20221896   ros_m(1) = 0.2255570073418735e+01_dp
    … …  
    20251899   ros_m(4) = 0.1093502252409163e+01_dp
    20261900!~~~> e_i  = coefficients for error estimator
    2027    ros_e(1) =- 0.2815431932141155_dp
    2028    ros_e(2) =- 0.7276199124938920e-01_dp
    2029    ros_e(3) =- 0.1082196201495311_dp
    2030    ros_e(4) =- 0.1093502252409163e+01_dp
    2031 !~~~> ros_elo  = estimator of local order -  the minimum between the
     1901   ros_e(1) = -0.2815431932141155_dp
     1902   ros_e(2) = -0.7276199124938920e-01_dp
     1903   ros_e(3) = -0.1082196201495311_dp
     1904   ros_e(4) = -0.1093502252409163e+01_dp
     1905!~~~> ros_elo  = estimator of local order - the minimum between the
    20321906!    main and the embedded scheme orders plus 1
    20331907   ros_elo  = 4.0_dp
    2034 !~~~> y_stage_i ~ y( t + h*alpha_i)
     1908!~~~> y_stage_i ~ y( t + h* alpha_i)
    20351909   ros_alpha(1) = 0.0_dp
    20361910   ros_alpha(2) = 0.1145640000000000e+01_dp
    20371911   ros_alpha(3) = 0.6552168638155900_dp
    20381912   ros_alpha(4) = ros_alpha(3)
    2039 !~~~> gamma_i = \sum_j  gamma_{i,j}
     1913!~~~> gamma_i = \sum_j  gamma_{i, j}
    20401914   ros_gamma(1) = 0.5728200000000000_dp
    2041    ros_gamma(2) =- 0.1769193891319233e+01_dp
     1915   ros_gamma(2) = -0.1769193891319233e+01_dp
    20421916   ros_gamma(3) = 0.7592633437920482_dp
    2043    ros_gamma(4) =- 0.1049021087100450_dp
     1917   ros_gamma(4) = -0.1049021087100450_dp
    20441918
    20451919  END SUBROUTINE ros4
    … …  
    20741948   ros_c(1) = 4.0_dp
    20751949   ros_c(2) = 1.0_dp
    2076    ros_c(3) =- 1.0_dp
     1950   ros_c(3) = -1.0_dp
    20771951   ros_c(4) = 1.0_dp
    2078    ros_c(5) =- 1.0_dp
    2079    ros_c(6) =- (8.0_dp/3.0_dp)
     1952   ros_c(5) = -1.0_dp
     1953   ros_c(6) = -(8.0_dp/3.0_dp)
    20801954
    20811955!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
    20821956!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
    2083    ros_newf(1) = .true.
    2084    ros_newf(2) = .false.
    2085    ros_newf(3) = .true.
    2086    ros_newf(4) = .true.
     1957   ros_newf(1) = .TRUE.
     1958   ros_newf(2) = .FALSE.
     1959   ros_newf(3) = .TRUE.
     1960   ros_newf(4) = .TRUE.
    20871961!~~~> m_i = coefficients for new step solution
    20881962   ros_m(1) = 2.0_dp
    … …  
    20951969   ros_e(3) = 0.0_dp
    20961970   ros_e(4) = 1.0_dp
    2097 !~~~> ros_elo  = estimator of local order -  the minimum between the
     1971!~~~> ros_elo  = estimator of local order - the minimum between the
    20981972!    main and the embedded scheme orders plus 1
    20991973   ros_elo  = 3.0_dp
    2100 !~~~> y_stage_i ~ y( t + h*alpha_i)
     1974!~~~> y_stage_i ~ y( t + h* alpha_i)
    21011975   ros_alpha(1) = 0.0_dp
    21021976   ros_alpha(2) = 0.0_dp
    21031977   ros_alpha(3) = 1.0_dp
    21041978   ros_alpha(4) = 1.0_dp
    2105 !~~~> gamma_i = \sum_j  gamma_{i,j}
     1979!~~~> gamma_i = \sum_j  gamma_{i, j}
    21061980   ros_gamma(1) = 0.5_dp
    21071981   ros_gamma(2) = 1.5_dp
    … …  
    21292003    ros_s = 6
    21302004
    2131 !~~~> y_stage_i ~ y( t + h*alpha_i)
     2005!~~~> y_stage_i ~ y( t + h* alpha_i)
    21322006    ros_alpha(1) = 0.000_dp
    21332007    ros_alpha(2) = 0.386_dp
    … …  
    21372011    ros_alpha(6) = 1.000_dp
    21382012
    2139 !~~~> gamma_i = \sum_j  gamma_{i,j}
     2013!~~~> gamma_i = \sum_j  gamma_{i, j}
    21402014    ros_gamma(1) = 0.2500000000000000_dp
    2141     ros_gamma(2) =- 0.1043000000000000_dp
     2015    ros_gamma(2) = -0.1043000000000000_dp
    21422016    ros_gamma(3) = 0.1035000000000000_dp
    2143     ros_gamma(4) =- 0.3620000000000023e-01_dp
     2017    ros_gamma(4) = -0.3620000000000023e-01_dp
    21442018    ros_gamma(5) = 0.0_dp
    21452019    ros_gamma(6) = 0.0_dp
    … …  
    21602034    ros_a(8) = 0.6019134481288629e+01_dp
    21612035    ros_a(9) = 0.1253708332932087e+02_dp
    2162     ros_a(10) =- 0.6878860361058950_dp
     2036    ros_a(10) = -0.6878860361058950_dp
    21632037    ros_a(11) = ros_a(7)
    21642038    ros_a(12) = ros_a(8)
    … …  
    21672041    ros_a(15) = 1.0_dp
    21682042
    2169     ros_c(1) =- 0.5668800000000000e+01_dp
    2170     ros_c(2) =- 0.2430093356833875e+01_dp
    2171     ros_c(3) =- 0.2063599157091915_dp
    2172     ros_c(4) =- 0.1073529058151375_dp
    2173     ros_c(5) =- 0.9594562251023355e+01_dp
    2174     ros_c(6) =- 0.2047028614809616e+02_dp
     2043    ros_c(1) = -0.5668800000000000e+01_dp
     2044    ros_c(2) = -0.2430093356833875e+01_dp
     2045    ros_c(3) = -0.2063599157091915_dp
     2046    ros_c(4) = -0.1073529058151375_dp
     2047    ros_c(5) = -0.9594562251023355e+01_dp
     2048    ros_c(6) = -0.2047028614809616e+02_dp
    21752049    ros_c(7) = 0.7496443313967647e+01_dp
    2176     ros_c(8) =- 0.1024680431464352e+02_dp
    2177     ros_c(9) =- 0.3399990352819905e+02_dp
     2050    ros_c(8) = -0.1024680431464352e+02_dp
     2051    ros_c(9) = -0.3399990352819905e+02_dp
    21782052    ros_c(10) = 0.1170890893206160e+02_dp
    21792053    ros_c(11) = 0.8083246795921522e+01_dp
    2180     ros_c(12) =- 0.7981132988064893e+01_dp
    2181     ros_c(13) =- 0.3152159432874371e+02_dp
     2054    ros_c(12) = -0.7981132988064893e+01_dp
     2055    ros_c(13) = -0.3152159432874371e+02_dp
    21822056    ros_c(14) = 0.1631930543123136e+02_dp
    2183     ros_c(15) =- 0.6058818238834054e+01_dp
     2057    ros_c(15) = -0.6058818238834054e+01_dp
    21842058
    21852059!~~~> m_i = coefficients for new step solution
    … …  
    22012075!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
    22022076!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
    2203     ros_newf(1) = .true.
    2204     ros_newf(2) = .true.
    2205     ros_newf(3) = .true.
    2206     ros_newf(4) = .true.
    2207     ros_newf(5) = .true.
    2208     ros_newf(6) = .true.
    2209 
    2210 !~~~> ros_elo  = estimator of local order -  the minimum between the
     2077    ros_newf(1) = .TRUE.
     2078    ros_newf(2) = .TRUE.
     2079    ros_newf(3) = .TRUE.
     2080    ros_newf(4) = .TRUE.
     2081    ros_newf(5) = .TRUE.
     2082    ros_newf(6) = .TRUE.
     2083
     2084!~~~> ros_elo  = estimator of local order - the minimum between the
    22112085!        main and the embedded scheme orders plus 1
    22122086    ros_elo = 4.0_dp
    … …  
    22712145!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
    22722146!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
    2273     ros_newf(1) = .true.
    2274     ros_newf(2) = .true.
    2275     ros_newf(3) = .true.
    2276     ros_newf(4) = .true.
    2277 
    2278 !~~~> ros_elo  = estimator of local order -  the minimum between the
     2147    ros_newf(1) = .TRUE.
     2148    ros_newf(2) = .TRUE.
     2149    ros_newf(3) = .TRUE.
     2150    ros_newf(4) = .TRUE.
     2151
     2152!~~~> ros_elo  = estimator of local order - the minimum between the
    22792153!        main and the embedded scheme orders plus 1
    22802154    ros_elo = 3.0_dp
    … …  
    22882162END SUBROUTINE rosenbrock
    22892163 
    2290 SUBROUTINE funtemplate( t,y,ydot)
     2164SUBROUTINE funtemplate( t, y, ydot)
    22912165!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    22922166!  Template for the ODE function call.
    … …  
    22942168!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    22952169!~~~> input variables
    2296    REAL(kind=dp):: t,y(nvar)
     2170   REAL(kind=dp):: t, y(nvar)
    22972171!~~~> output variables
    22982172   REAL(kind=dp):: ydot(nvar)
    … …  
    23022176   told = time
    23032177   time = t
    2304    CALL fun( y,fix,rconst,ydot)
     2178   CALL fun( y, fix, rconst, ydot)
    23052179   time = told
    23062180
    23072181END SUBROUTINE funtemplate
    23082182 
    2309 SUBROUTINE jactemplate( t,y,jcb)
     2183SUBROUTINE jactemplate( t, y, jcb)
    23102184!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    23112185!  Template for the ODE Jacobian call.
    … …  
    23132187!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    23142188!~~~> input variables
    2315     REAL(kind=dp):: t,y(nvar)
     2189    REAL(kind=dp):: t, y(nvar)
    23162190!~~~> output variables
    23172191#ifdef full_algebra   
    2318     REAL(kind=dp):: jv(lu_nonzero),jcb(nvar,nvar)
     2192    REAL(kind=dp):: jv(lu_nonzero), jcb(nvar, nvar)
    23192193#else
    23202194    REAL(kind=dp):: jcb(lu_nonzero)
    … …  
    23232197    REAL(kind=dp):: told
    23242198#ifdef full_algebra   
    2325     INTEGER :: i,j
     2199    INTEGER :: i, j
    23262200#endif   
    23272201
    … …  
    23292203    time = t
    23302204#ifdef full_algebra   
    2331     CALL jac_sp(y,fix,rconst,jv)
    2332     DO j=1,nvar
    2333       DO i=1,nvar
    2334          jcb(i,j) = 0.0_dp
     2205    CALL jac_sp(y, fix, rconst, jv)
     2206    DO j=1, nvar
     2207      DO i=1, nvar
     2208         jcb(i, j) = 0.0_dp
    23352209      ENDDO
    23362210    ENDDO
    2337     DO i=1,lu_nonzero
    2338        jcb(lu_irow(i),lu_icol(i)) = jv(i)
     2211    DO i=1, lu_nonzero
     2212       jcb(lu_irow(i), lu_icol(i)) = jv(i)
    23392213    ENDDO
    23402214#else
    2341     CALL jac_sp( y,fix,rconst,jcb)
     2215    CALL jac_sp( y, fix, rconst, jcb)
    23422216#endif   
    23432217    time = told
    … …  
    23452219END SUBROUTINE jactemplate
    23462220 
    2347 SUBROUTINE chem_gasphase_integrate (time_step_len,conc,tempk,photo,ierrf,xnacc,xnrej,istatus,l_debug,pe)
     2221  SUBROUTINE kppdecomp( jvs, ier)                                   
     2222  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     2223  !        sparse lu factorization                                   
     2224  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     2225  ! loop expansion generated by kp4                                   
     2226                                                                     
     2227    INTEGER  :: ier                                                   
     2228    REAL(kind=dp):: jvs(lu_nonzero), w(nvar), a             
     2229    INTEGER  :: k, kk, j, jj                                         
     2230                                                                     
     2231    a = 0.                                                           
     2232    ier = 0                                                           
     2233                                                                     
     2234!   i = 1
     2235!   i = 2
     2236    RETURN                                                           
     2237                                                                     
     2238  END SUBROUTINE kppdecomp                                           
     2239 
     2240SUBROUTINE chem_gasphase_integrate (time_step_len, conc, tempi, qvapi, fakti, photo, ierrf, xnacc, xnrej, istatus, l_debug, pe, &
     2241                     icntrl_i, rcntrl_i)
    23482242                                                                   
    23492243  IMPLICIT NONE                                                     
    23502244                                                                   
    2351   REAL(dp), INTENT(IN)                    :: time_step_len           
    2352   REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: conc                   
    2353   REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: photo                   
    2354   REAL(dp), DIMENSION(:), INTENT(IN)      :: tempk                   
    2355   INTEGER, INTENT(OUT), OPTIONAL          :: ierrf(:)               
    2356   INTEGER, INTENT(OUT), OPTIONAL          :: xnacc(:)               
    2357   INTEGER, INTENT(OUT), OPTIONAL          :: xnrej(:)               
    2358   INTEGER, INTENT(INOUT), OPTIONAL        :: istatus(:)             
    2359   INTEGER, INTENT(IN), OPTIONAL           :: pe                     
    2360   LOGICAL, INTENT(IN), OPTIONAL           :: l_debug                 
     2245  REAL(dp), INTENT(IN)                  :: time_step_len           
     2246  REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: conc                   
     2247  REAL(dp), DIMENSION(:, :), INTENT(IN)   :: photo                   
     2248  REAL(dp), DIMENSION(:), INTENT(IN)     :: tempi                   
     2249  REAL(dp), DIMENSION(:), INTENT(IN)     :: qvapi                   
     2250  REAL(dp), DIMENSION(:), INTENT(IN)     :: fakti                   
     2251  INTEGER, INTENT(OUT), OPTIONAL        :: ierrf(:)               
     2252  INTEGER, INTENT(OUT), OPTIONAL        :: xnacc(:)               
     2253  INTEGER, INTENT(OUT), OPTIONAL        :: xnrej(:)               
     2254  INTEGER, INTENT(INOUT), OPTIONAL      :: istatus(:)             
     2255  INTEGER, INTENT(IN), OPTIONAL         :: pe                     
     2256  LOGICAL, INTENT(IN), OPTIONAL         :: l_debug                 
     2257  INTEGER, DIMENSION(nkppctrl), INTENT(IN), OPTIONAL  :: icntrl_i         
     2258  REAL(dp), DIMENSION(nkppctrl), INTENT(IN), OPTIONAL  :: rcntrl_i         
    23612259                                                                   
    23622260  INTEGER                                 :: k   ! loop variable     
    2363   REAL(dp)                                :: dt                     
    2364   INTEGER, DIMENSION(20)                  :: istatus_u               
    2365   INTEGER                                 :: ierr_u                 
    2366                                                                    
    2367                                                                    
    2368   if (present (istatus)) istatus = 0                             
    2369                                                                    
    2370   DO k=1,vl_glo,vl_dim                                             
     2261  REAL(dp)                               :: dt                     
     2262  INTEGER, DIMENSION(20)                :: istatus_u               
     2263  INTEGER                                :: ierr_u                 
     2264  INTEGER                                :: istatf                 
     2265  INTEGER                                :: vl_dim_lo               
     2266                                                                   
     2267                                                                   
     2268  IF (PRESENT (istatus)) istatus = 0                             
     2269  IF (PRESENT (icntrl_i)) icntrl  = icntrl_i                     
     2270  IF (PRESENT (rcntrl_i)) rcntrl  = rcntrl_i                     
     2271                                                                   
     2272  vl_glo = size(tempi, 1)                                           
     2273                                                                   
     2274  vl_dim_lo = vl_dim                                               
     2275  DO k=1, vl_glo, vl_dim_lo                                           
    23712276    is = k                                                         
    2372     ie = min(k+vl_dim-1,vl_glo)                                     
    2373     vl = ie-is+1                                                   
    2374                                                                    
    2375     c(:) = conc(is,:)                                             
    2376                                                                    
    2377     temp = tempk(is)                                             
    2378                                                                    
    2379     phot(:) = photo(is,:)                                             
     2277    ie = min(k+ vl_dim_lo-1, vl_glo)                                 
     2278    vl = ie-is+ 1                                                   
     2279                                                                   
     2280    c(:) = conc(is, :)                                           
     2281                                                                   
     2282    temp = tempi(is)                                             
     2283                                                                   
     2284    qvap = qvapi(is)                                             
     2285                                                                   
     2286    fakt = fakti(is)                                             
     2287                                                                   
     2288    CALL initialize                                                 
     2289                                                                   
     2290    phot(:) = photo(is, :)                                           
    23802291                                                                   
    23812292    CALL update_rconst                                             
    … …  
    23842295                                                                   
    23852296    ! integrate from t=0 to t=dt                                   
    2386     CALL integrate(0._dp, dt,icntrl,rcntrl,istatus_u = istatus_u,ierr_u=ierr_u)
    2387                                                                    
    2388                                                                    
    2389     IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN                       
    2390        IF (l_debug) CALL error_output(conc(is,:),ierr_u,pe)           
    2391     ENDIF                                                             
    2392     conc(is,:) = c(:)                                                 
    2393                                                                    
    2394     ! Return Diagnostic Information                                 
    2395                                                                    
    2396     if(PRESENT(ierrf))    ierrf(is) = ierr_u                     
    2397     if(PRESENT(xnacc))    xnacc(is) = istatus_u(4)               
    2398     if(PRESENT(xnrej))    xnrej(is) = istatus_u(5)               
    2399                                                                    
    2400     if (PRESENT (istatus)) then                                   
    2401       istatus(1:8) = istatus(1:8) + istatus_u(1:8)                 
     2297    CALL integrate(0._dp, dt, icntrl, rcntrl, istatus_u = istatus_u, ierr_u=ierr_u)
     2298                                                                   
     2299                                                                   
     2300   IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN                       
     2301      IF (l_debug) CALL error_output(conc(is, :), ierr_u, pe)         
     2302   ENDIF                                                             
     2303                                                                     
     2304    conc(is, :) = c(:)                                               
     2305                                                                   
     2306    ! RETURN diagnostic information                                 
     2307                                                                   
     2308    IF (PRESENT(ierrf))   ierrf(is) = ierr_u                     
     2309    IF (PRESENT(xnacc))   xnacc(is) = istatus_u(4)               
     2310    IF (PRESENT(xnrej))   xnrej(is) = istatus_u(5)               
     2311                                                                   
     2312    IF (PRESENT (istatus)) THEN                                   
     2313      istatus(1:8) = istatus(1:8) + istatus_u(1:8)               
    24022314    ENDIF                                                         
    24032315                                                                   
    … …  
    24072319! Deallocate input arrays                                           
    24082320                                                                   
    2409   if (ALLOCATED(temp))   DEALLOCATE(temp)   
    2410                                                                    
    2411   data_loaded = .false.                                             
     2321                                                                   
     2322  data_loaded = .FALSE.                                             
    24122323                                                                   
    24132324  RETURN                                                           
    24142325END SUBROUTINE chem_gasphase_integrate                                       
    2415  
    2416   SUBROUTINE fill_temp(status,array)
    2417  
    2418     INTEGER, INTENT(OUT)               :: status
    2419     REAL(dp), INTENT(IN), DIMENSION(:) :: array
    2420  
    2421     status = 0
    2422     IF (.not. ALLOCATED(temp)) &
    2423        ALLOCATE(temp(size(array)))
    2424  
    2425     IF (data_loaded .AND. (vl_glo /= size(array,1))) THEN
    2426        status = 1
    2427        RETURN
    2428     END IF
    2429  
    2430     vl_glo = size(array,1)
    2431     temp = array
    2432     data_loaded = .TRUE.
    2433  
    2434     RETURN
    2435  
    2436   END   SUBROUTINE fill_temp
    24372326
    24382327END MODULE chem_gasphase_mod
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |