Home


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

Changeset 3298 for palm/trunk/SOURCE


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:
19 edited
5 copied

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/Makefile

    r3294 r3298  
    2525# -----------------
    2626# $Id$
     27# Added missing dependencies and replaced blanks with tabs (forkel)
     28# Added chemistry emission module (Russo)
     29# Added chemistry_model_mod to flow_statistics (basit)
     30#
     31# 3294 2018-10-01 02:37:10Z raasch
    2732# changes related to modularization of the ocean mode,
    2833# bugfix: dependency to advec_ws was missed in chemistry_model_mod
    … …  
    473478        check_open.f90 \
    474479        check_parameters.f90 \
     480        chem_emissions_mod.f90 \
    475481        chem_gasphase_mod.f90 \
    476482        chemistry_model_mod.f90 \
    … …  
    576582        random_gauss.f90 \
    577583        random_generator_parallel_mod.f90 \
    578         read_restart_data_mod.f90 \
     584        read_restart_data_mod.f90 \
    579585        run_control.f90 \
    580586        set_slicer_attributes_dvrp.f90 \
    … …  
    636642        wind_turbine_model_mod.f90 \
    637643        wrd_write_string.f90 \
    638         write_restart_data_mod.f90
     644        write_restart_data_mod.f90
    639645
    640646
    … …  
    757763        bulk_cloud_model_mod.o \
    758764        chemistry_model_mod.o \
     765        chem_emissions_mod.o \
    759766        gust_mod.o \
    760767        init_vertical_profiles.o \
    … …  
    778785        vertical_nesting_mod.o \
    779786        wind_turbine_model_mod.o
     787chem_emissions_mod.o: \
     788        chem_gasphase_mod.o \
     789        chem_modules.o \
     790        date_and_time_mod.o \
     791        mod_kinds.o \
     792        modules.o \
     793        netcdf_data_input_mod.o \
     794        surface_mod.o
    780795chemistry_model_mod.o: \
    781796        advec_ws.o \
    782797        chem_gasphase_mod.o \
    783798        chem_modules.o \
     799        date_and_time_mod.o \
    784800        diffusion_s.o \
    785801        mod_kinds.o \
    786802        modules.o \
    787803        netcdf_data_input_mod.o \
    788         surface_mod.o
     804        surface_mod.o \
     805        user_module.o
    789806chem_gasphase_mod.o: \
    790807        mod_kinds.o \
    … …  
    831848        basic_constants_and_equations_mod.o \
    832849        bulk_cloud_model_mod.o \
     850        chemistry_model_mod.o \
    833851        cpulog_mod.o \
    834852        mod_kinds.o \
    … …  
    868886        basic_constants_and_equations_mod.o \
    869887        bulk_cloud_model_mod.o \
     888        chemistry_model_mod.o\
    870889        cpulog_mod.o \
    871890        gust_mod.o \
    … …  
    942961        basic_constants_and_equations_mod.o \
    943962        bulk_cloud_model_mod.o \
     963        chemistry_model_mod.o \
    944964        cpulog_mod.o \
    945965        gust_mod.o \
    … …  
    960980        basic_constants_and_equations_mod.o \
    961981        bulk_cloud_model_mod.o \
     982        chemistry_model_mod.o \
    962983        cpulog_mod.o \
    963984        date_and_time_mod.o \
    … …  
    9861007        basic_constants_and_equations_mod.o \
    9871008        bulk_cloud_model_mod.o \
    988         chemistry_model_mod.o \
     1009        chem_emissions_mod.o \
    9891010        cpulog_mod.o \
    9901011        disturb_heatflux.o \
    … …  
    12171238        basic_constants_and_equations_mod.o \
    12181239        mod_kinds.o \
    1219         modules.o
     1240        modules.o \
     1241        time_to_string.o
    12201242modules.o: \
    12211243        mod_kinds.o
    … …  
    12331255        mod_kinds.o
    12341256netcdf_data_input_mod.o: \
     1257        chem_modules.o \
    12351258        cpulog_mod.o \
    12361259        mod_kinds.o \
    … …  
    12741297        chem_photolysis_mod.o \
    12751298        cpulog_mod.o \
     1299        date_and_time_mod.o \
    12761300        land_surface_model_mod.o \
    12771301        mod_kinds.o \
    … …  
    12821306        pmc_particle_interface.o \
    12831307        surface_layer_fluxes_mod.o \
     1308        time_to_string.o \
    12841309        write_restart_data_mod.o
    12851310parin.o: \
    … …  
    15451570        buoyancy.o \
    15461571        calc_mean_profile.o \
     1572        chem_emissions_mod.o \
    15471573        chemistry_model_mod.o \
    15481574        chem_modules.o \
    … …  
    15581584        modules.o \
    15591585        multi_agent_system_mod.o \
    1560  ocean_mod.o \
     1586        ocean_mod.o \
    15611587        pmc_interface_mod.o \
    15621588        prognostic_equations.o \
    … …  
    15671593        surface_mod.o \
    15681594        synthetic_turbulence_generator_mod.o \
     1595        time_to_string.o \
    15691596        turbulence_closure_mod.o\
    15701597        urban_surface_mod.o \
    … …  
    15841611        radiation_model_mod.o \
    15851612        surface_layer_fluxes_mod.o \
     1613        time_to_string.o \
    15861614        urban_surface_mod.o
    15871615time_to_string.o: \

  • palm/trunk/SOURCE/advec_ws.f90

    r3294 r3298  
    2525! -----------------
    2626! $Id$
     27! Add formatted note from ketelsen about statistics for chemistry
     28!
     29! 3294 2018-10-01 02:37:10Z raasch
    2730! ocean renamed ocean_mode
    2831!
    … …  
    16971700                    ) * weight_substep(intermediate_timestep_count)
    16981701             ENDDO
     1702
     1703!          CASE ( 'kc' )
     1704          !kk Has to be implemented for kpp chemistry
     1705
    16991706
    17001707         END SELECT

  • palm/trunk/SOURCE/check_parameters.f90

    r3294 r3298  
    2525! -----------------
    2626! $Id: check_parameters.f90 2520 2017-10-05 13:50:26Z gronemeier &
     27! - Minor formatting and remarks
     28! - Call of chem_emissions_check_parameters commented
     29!   (preliminary, as there is still a bug in this subroutine) (forkel)
     30! - Corrected call for chemistry profile output added (basit)
     31! - Call for init_vertical_profiles for chemistry removed (basit)
     32! - Call for chem_init_profiles removed (basit)
     33! - Setting of initial profiles for chemstry (constant and call of
     34!   init_vertical_profiles) commented (forkel)
     35!
     36! 2520 2017-10-05 13:50:26Z gronemeier
    2737! changes concerning modularization of ocean option,
    2838! init_vertical_profiles moved to separate file to avoid circular dependency
    … …  
    679689               bcm_check_data_output_pr
    680690
     691    USE chem_emissions_mod,                                                    &
     692        ONLY:  chem_emissions_check_parameters
     693
     694    USE chem_modules
     695
    681696    USE chemistry_model_mod,                                                   &
    682697        ONLY:  chem_boundary_conds, chem_check_data_output,                    &
    683                chem_check_data_output_pr, chem_species
    684 
    685     USE chem_modules
     698               chem_check_data_output_pr, chem_check_parameters, chem_species
    686699
    687700    USE control_parameters
    … …  
    14431456    ENDIF
    14441457
    1445 !
     1458!   Check for chem_emission_mod parameters setting
     1459!   IF ( air_chemistry ) CALL chem_emissions_check_parameters   ! forkel preliminary
     1460
     1461
     1462!   Check for chemitry_model_mod parameters setting
     1463    IF ( air_chemistry )  CALL chem_check_parameters
     1464
     1465
    14461466!-- Check the module settings
    14471467    IF ( ocean_mode )  CALL ocean_check_parameters
    … …  
    14681488       IF ( humidity       )  q_init  = q_surface
    14691489       IF ( passive_scalar )  s_init  = s_surface
    1470        IF ( air_chemistry )  THEN
    1471          DO lsp = 1, nvar
    1472             chem_species(lsp)%conc_pr_init = cs_surface(lsp)     
    1473          ENDDO
    1474        ENDIF
     1490!
     1491!--    TODO
     1492!--    Russo: Is done in chem_init and would overwrite what is done there
     1493!--    --> kanani: Revise
     1494!      IF ( air_chemistry )  THEN
     1495!        DO lsp = 1, nvar
     1496!           chem_species(lsp)%conc_pr_init = cs_surface(lsp)     
     1497!        ENDDO
     1498!      ENDIF
    14751499!
    14761500!--
    … …  
    16831707       ENDIF
    16841708!
     1709!--    TODO
    16851710!--    Compute initial chemistry profile using the given chemical species gradients
    1686        IF ( air_chemistry )  THEN                                                         
    1687           DO  lsp = 1, nvar
    1688             CALL init_vertical_profiles( cs_vertical_gradient_level_ind(lsp,:),&
    1689                                          cs_vertical_gradient_level(lsp,:),    &
    1690                                          cs_vertical_gradient(lsp,:),          &
    1691                                          chem_species(lsp)%conc_pr_init,       &
    1692                                          cs_surface(lsp), bc_cs_t_val(lsp) )
    1693           ENDDO
    1694        ENDIF
     1711!--    Russo: Is done in chem_init --> kanani: Revise
    16951712
    16961713    ENDIF

  • 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

  • palm/trunk/SOURCE/data_output_2d.f90

    r3294 r3298  
    2525! -----------------
    2626! $Id$
     27! minor formatting (kanani)
     28! chem_data_output_2d subroutine added (basit)
     29!
     30! 3294 2018-10-01 02:37:10Z raasch
    2731! changes concerning modularization of ocean option
    2832!
    … …  
    252256        ONLY:  bulk_cloud_model, bcm_data_output_2d
    253257
     258    USE chemistry_model_mod,                                                   &
     259        ONLY:  chem_data_output_2d
     260
    254261    USE control_parameters,                                                    &
    255         ONLY:  data_output_2d_on_each_pe, data_output_xy,                      &
     262        ONLY:  air_chemistry, data_output_2d_on_each_pe, data_output_xy,       &
    256263               data_output_xz, data_output_yz, do2d,                           &
    257264               do2d_xy_last_time, do2d_xy_time_count,                          &
    … …  
    265272    USE cpulog,                                                                &
    266273        ONLY:  cpu_log, log_point
    267 
     274       
    268275    USE gust_mod,                                                              &
    269276        ONLY:  gust_data_output_2d, gust_module_enabled
    … …  
    13261333                ENDIF
    13271334
     1335                IF ( .NOT. found  .AND.  air_chemistry )  THEN
     1336                   CALL chem_data_output_2d( av, do2d(av,if), found, grid, mode, &
     1337                                             local_pf, two_d, nzb_do, nzt_do, fill_value )
     1338                ENDIF
    13281339!
    13291340!--             User defined quantities

  • palm/trunk/SOURCE/data_output_mask.f90

    r3274 r3298  
    2525! -----------------
    2626! $Id$
     27! Minor formatting (kanani)
     28! Call for chem_data_output_mask added (basit)
     29!
     30! 3274 2018-09-24 15:42:55Z knoop
    2731! Modularization of all bulk cloud physics code components
    2832!
    … …  
    133137        ONLY:  e, nc, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa,  &
    134138               tend, u, v, vpt, w, d_exner
    135    
     139
    136140    USE averaging,                                                             &
    137141        ONLY:  e_av, lpt_av, nc_av, nr_av, p_av, pc_av, pr_av, pt_av, q_av,    &
    … …  
    141145    USE basic_constants_and_equations_mod,                                     &
    142146        ONLY:  lv_d_cp
    143    
     147
     148    USE chemistry_model_mod,                                                   &
     149        ONLY:  chem_data_output_mask
     150
    144151    USE control_parameters,                                                    &
    145         ONLY:  domask, domask_no, domask_time_count, mask_i,                   &
     152        ONLY:  air_chemistry, domask, domask_no, domask_time_count, mask_i,    &
    146153               mask_j, mask_k, mask_size, mask_size_l, mask_start_l,           &
    147154               max_masks, message_string, mid, nz_do3d, simulated_time
    … …  
    151158    USE indices,                                                               &
    152159        ONLY:  nbgp, nxl, nxr, nyn, nys, nzb
    153        
     160
    154161    USE kinds
    155162
    … …  
    524531
    525532          CASE DEFAULT
    526 
    527533!
    528534!--          Radiation quantity
    … …  
    532538             ENDIF
    533539
     540             IF ( air_chemistry )  THEN
     541                CALL chem_data_output_mask(av, domask(mid,av,ivar), found,     &
     542                                           local_pf )
     543             ENDIF
    534544!
    535545!--          User defined quantity

  • palm/trunk/SOURCE/date_and_time_mod.f90

    r2718 r3298  
    2525! -----------------
    2626! $Id$
     27! - Minor formatting (kanani)
     28! - Added Routines for DEFAULT mode of chemistry emissions (Russo)
     29! - Added routine for reading-in real dates: format has to be in DDMMYYYY and
     30!   passed in the namelist parameter date_init (Russo)
     31! - Added calculation of several time indices useful in several routines
     32!   of the model (Russo)
     33!
     34! 2718 2018-01-02 08:49:38Z maronga
    2735! Corrected "Former revisions" section
    2836!
    … …  
    4654!> This routine calculates all needed information on date and time used by
    4755!> other modules
     56!> @todo Further testing and revision of routines for updating indices of
     57!>       emissions in the default mode
     58!> @todo Add routine for recognizing leap years
     59!> @todo Add recognition of exact days of week (Monday, Tuesday, etc.)
     60!> @todo Reconsider whether to remove day_of_year_init from the namelist: we
     61!>       already implemented changes for calculating it from date_init in
     62!>       calc_date_and_time
    4863!------------------------------------------------------------------------------!
    4964 MODULE date_and_time_mod
    50  
    51     USE control_parameters,                                                    &
    52         ONLY: time_since_reference_point
    53  
     65
     66    USE control_parameters,                                                   &
     67        ONLY:  coupling_start_time, days_since_reference_point,               &
     68               message_string, simulated_time, time_since_reference_point
     69
    5470    USE kinds
    5571
     72
    5673    IMPLICIT NONE
    5774
    5875    PRIVATE
    59    
    60     PUBLIC   calc_date_and_time, d_hours_day, d_seconds_hour, d_seconds_year,  &
    61              day_of_year, day_of_year_init, time_utc, time_utc_init
    62 
    63 
    64     INTEGER(iwp) ::  day_of_year              !< day of the year (DOY)
    65     INTEGER(iwp) ::  day_of_year_init = 172   !< DOY at model start (default: 21 June)
    66 
    67     REAL(wp) ::  time_utc                     !< current model time in UTC
    68     REAL(wp) ::  time_utc_init = 43200.0_wp   !< UTC time at model start
    69 
     76
     77!-- Variables Declaration
     78
     79    INTEGER(iwp)        ::  day_of_year      = 0   !< day of the year (DOY)
     80    INTEGER(iwp)        ::  day_of_year_init = 0   !< DOY at model start (default: 0)
     81
     82    ! --- Most of these indices are updated by the routine calc_date_and_time according to the current date and time of the simulation
     83    INTEGER(iwp)        ::  hour_of_year = 1                        !< hour of the current year (1:8760(8784))
     84    INTEGER(iwp)        ::  hour_of_day=1                           !< hour of the current day (1:24)
     85    INTEGER(iwp)        ::  day_of_month=0                          !< day of the current month (1:31)
     86    INTEGER(iwp)        ::  month_of_year=0                         !< month of the current year (1:12)
     87    INTEGER(iwp)        ::  current_year=0                          !< current year
     88    INTEGER(iwp)        ::  hour_call_emis=0                        !< index used to call the emissions just once every hour
     89
     90    INTEGER(iwp)        ::  ihour           
     91    INTEGER(iwp)        ::  index_mm                                !< index months of the default emission mode
     92    INTEGER(iwp)        ::  index_dd                                !< index days of the default emission mode
     93    INTEGER(iwp)        ::  index_hh                                !< index hours of the default emission mode
     94
     95    REAL(wp)            ::  time_utc                     !< current model time in UTC
     96    REAL(wp)            ::  time_utc_emis                !< current emission module time in UTC
     97    REAL(wp)            ::  time_utc_init = 43200.0_wp   !< UTC time at model start
     98    REAL(wp)            ::  time_update                  !< used to calculate actual second of the simulation
     99 
    70100    REAL(wp), PARAMETER ::  d_hours_day    = 1.0_wp / 24.0_wp       !< inverse of hours per day (1/24)
    71101    REAL(wp), PARAMETER ::  d_seconds_hour = 1.0_wp / 3600.0_wp     !< inverse of seconds per hour (1/3600)
    72102    REAL(wp), PARAMETER ::  d_seconds_year = 1.0_wp / 31536000.0_wp !< inverse of the seconds per year (1/(365*86400))
    73103   
     104    CHARACTER(len=8)    ::  date_init = "31122018"                  !< Starting date of simulation: We selected this because it was a monday
     105 
     106    !> --- Parameters
     107    INTEGER, PARAMETER, DIMENSION(12) :: days = (/31,28,31,30,31,30,31,31,30,31,30,31/) ! total number of days for each month (no leap year)
     108
    74109    SAVE
    75110
     111!-- INTERFACES PART
     112    !-- Read initial day and time of simulation
     113    INTERFACE init_date_and_time
     114       MODULE PROCEDURE init_date_and_time
     115    END INTERFACE init_date_and_time
     116
     117    !-- Get hour index in the DEAFULT case of chemistry emissions :
     118    INTERFACE time_default_indices
     119       MODULE PROCEDURE time_mdh_indices
     120       MODULE PROCEDURE time_hour_indices
     121    END INTERFACE time_default_indices
     122
     123    !-- Calculate current date and time
     124    INTERFACE calc_date_and_time
     125       MODULE PROCEDURE calc_date_and_time
     126    END INTERFACE
     127
     128
     129    !-- Public Interfaces
     130    PUBLIC calc_date_and_time, time_default_indices, init_date_and_time
     131
     132    !-- Public Variables
     133    PUBLIC date_init, d_hours_day, d_seconds_hour, d_seconds_year,               &
     134           day_of_year, day_of_year_init, time_utc, time_utc_init, day_of_month, &
     135           month_of_year, index_mm, index_dd, index_hh, hour_of_day, hour_of_year, &
     136           hour_call_emis
     137
    76138 CONTAINS
     139
     140
     141 !------------------------------------------------------------------------------!
     142 !> Reads starting date from namelist
     143 !------------------------------------------------------------------------------!
     144 
     145    SUBROUTINE init_date_and_time
     146
     147       IMPLICIT NONE
     148
     149       IF  (day_of_year_init == 0) THEN
     150          ! Day of the month at starting time
     151          READ(UNIT=date_init(1:2),fmt=*)day_of_month
     152
     153          ! Month of the year at starting time
     154          READ(UNIT=date_init(3:4),fmt=*)month_of_year
     155
     156          ! Year at starting time
     157          READ(UNIT=date_init(5:8),fmt=*)current_year
     158       
     159       ENDIF
     160
     161    END SUBROUTINE init_date_and_time
    77162
    78163!------------------------------------------------------------------------------!
    79164! Description:
    80165! ------------
    81 !> Calculate current day of the year and time in UTC
     166!> Calculate current date and time of the simulation
    82167!------------------------------------------------------------------------------!
    83168 
    … …  
    86171       IMPLICIT NONE
    87172
    88 !
    89 !--    Calculate current day of the year
    90        day_of_year = day_of_year_init + INT(FLOOR( (time_utc_init + time_since_reference_point)&
    91                                / 86400.0_wp ), KIND=iwp)
    92                        
     173!--    Variables Definition
     174       INTEGER                          :: i_mon       !< Index for going through the different months
     175
     176       !> Update simulation time in seconds
     177       time_update = simulated_time-coupling_start_time
     178
     179!--    Calculate current day of the simulated time
     180       days_since_reference_point=INT(FLOOR( (time_utc_init + time_update) &
     181                               / 86400.0_wp ) )
     182
     183!--    Calculate actual UTC time                       
    93184       time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
    94185       
     186!sB    PRILIMINARY workaround for time_utc bug concerning time_since_reference_point:
     187       time_utc_emis = MOD((time_utc_init + time_update), 86400.0_wp)     
     188
     189!--    Calculate initial day of the year: it is calculated only once. In fact, day_of_year_init is initialized to 0 and then a positive value is passed. This condition is also called only when day_of_year_init is not given in the namelist.
     190
     191       IF ( day_of_year_init == 0 ) THEN
     192
     193          !> Condition for printing an error when date_init is not provided when day_of_year_init is not given in the namelist or when the format of the date is not the one required by PALM.
     194          IF ( day_of_month .GT. 0 .AND. day_of_month .LE. 31 .AND. month_of_year .GT. 0 .AND. month_of_year .LE. 12) THEN
     195       
     196             IF ( month_of_year == 1 ) THEN  !!month of year is read in input
     197
     198                day_of_year_init = day_of_month
     199
     200             ELSE
     201
     202                day_of_year_init= SUM(days( 1:(month_of_year-1) )) + day_of_month  !day_of_month is read in input in this case
     203
     204             ENDIF
     205!kanani: Revise, we cannot force users to provide date_init, maybe set a default value?
     206!           ELSE
     207!
     208!              message_string = 'date_init not provided in the namelist or'             //          &
     209!                               ' given in the wrong format: MUST BE DDMMYYYY'                       
     210!              CALL message( 'calc_date_and_time', 'DT0100', 2, 2, 0, 6, 0 )
     211     
     212          ENDIF
     213
     214       ENDIF
     215
     216      !-- Calculate actual hour of the day: the first hour of the day is from 00:00:00 to 00:59:59.
     217
     218      hour_of_day = INT( FLOOR( time_utc_emis/3600.0_wp ) ) + 1
     219
     220!--    Calculate current day of the year !TBD: considetr leap years
     221       IF ( (day_of_year_init + days_since_reference_point)  .GT. 365 ) THEN
     222
     223          day_of_year=INT(MOD((day_of_year_init + days_since_reference_point), 365.0_wp))
     224
     225       ELSE
     226         
     227          day_of_year = day_of_year_init + days_since_reference_point
     228
     229       ENDIF
     230
     231!
     232!--    Calculate current hour of the year
     233       hour_of_year = ( (day_of_year-1) * 24 ) + hour_of_day  !> actual hour of the year
     234       
     235
     236!
     237!--    UPDATE actual day of the month and month of the year
     238       !> --------------------------------------------------------------------------------
     239       !> The first case is when date_init is not provided: we only know day_of_year_init     
     240       IF ( month_of_year == 0 .AND. day_of_month == 0) THEN
     241
     242          !> The first case is when date_init is not provided: we only know day_of_year_init
     243          !DO i_mon=1,12
     244             !IF (day_of_year .LE. SUM(days(1:i_mon))) THEN
     245          IF ( day_of_year .LE. 31 ) THEN
     246
     247             month_of_year=1
     248             day_of_month=day_of_year
     249
     250          ELSE
     251
     252             DO i_mon=2,12   !january is considered in the first case
     253                IF ( day_of_year .LE. SUM(days(1:i_mon)) .AND. day_of_year .GT. SUM(days(1:(i_mon-1))) ) THEN
     254           
     255                   month_of_year=i_mon
     256
     257                   day_of_month=INT(MOD(day_of_year, SUM(days(1:(i_mon-1)))))
     258
     259                   GOTO 38
     260
     261                ENDIF
     262
     263             38 ENDDO
     264          ENDIF
     265       !> --------------------------------------------------------------------------------
     266       !> in the second condition both day of month and month_of_year are either given in input (passed to date_init) or we are in some day successive to the initial one, so that day_of_month has already be computed in previous step
     267       !>TBD: something to calculate the current year is missing
     268       ELSEIF ( day_of_month .GT. 0 .AND. day_of_month .LE. 31 .AND. month_of_year .GT. 0 .AND. month_of_year .LE. 12) THEN
     269 
     270          !> calculate month_of_year. TBD: test the condition when day_of_year==31
     271 
     272          IF (day_of_year==1) THEN  !> this allows to turn from december to January when passing from a year to another
     273 
     274             month_of_year = 1
     275       
     276          ELSE IF (day_of_year .GT. 1 .AND. day_of_year .GT. SUM(days(1:month_of_year))) THEN
     277
     278             month_of_year = month_of_year + 1
     279
     280          ENDIF
     281
     282          !> calculate day_of_month
     283          IF ( month_of_year == 1 ) THEN
     284           
     285            day_of_month=day_of_year
     286
     287          ELSE
     288
     289            day_of_month=INT(MOD(day_of_year, SUM(days(1:(month_of_year-1)))))
     290
     291          ENDIF
     292
     293
     294
     295
     296      ELSE
     297
     298          !> Condition when date_init is provided but it is given in the wrong format
     299          message_string = 'date_init not provided in the namelist or'            //          &
     300                              ' given in the wrong format: MUST BE DDMMYYYY'                 
     301          CALL message( 'calc_date_and_time', 'DT0101', 2, 2, 0, 6, 0 )
     302
     303      ENDIF       
     304     
     305    END SUBROUTINE calc_date_and_time
     306
     307!TBD: check the routines used for update of emission indices in the DEFAULT MODE
     308!------------------------------------------------------------------------------!
     309! Description:
     310! ------------
     311!> This routine determines the time factor index in the mdh case of the DEFAULT chemistry emissions mode.
     312!------------------------------------------------------------------------------!
     313
     314 SUBROUTINE time_mdh_indices(daytype_mdh,mo, dd, hh, index_mm, index_dd, index_hh)
     315
     316    USE indices
     317
     318    IMPLICIT NONE
     319
     320    !> IN/OUTPUT
     321    INTEGER, INTENT(INOUT) :: mo          !> Month of year
     322    INTEGER, INTENT(INOUT) :: dd          !> Day of month
     323    INTEGER, INTENT(INOUT) :: hh          !> Hour of day
     324    INTEGER, INTENT(INOUT) :: index_mm    !> Index Month
     325    INTEGER, INTENT(INOUT) :: index_dd    !> Index Day
     326    INTEGER, INTENT(INOUT) :: index_hh    !> Index Hour
     327
     328    CHARACTER(len=80), INTENT(INOUT) :: daytype_mdh !> type of the day in mdh mode: one of 1-WORKDAY
     329                                                    !                                      2-WEEKEND
     330                                                    !                                      3-HOLIDAY
     331
     332    REAL(wp)                         :: frac_day=0 
     333
     334    !> ------------------------------------------------------------------------
     335
     336    INTEGER                          :: weekday
     337
     338    !> CONSTANTS
     339    INTEGER, PARAMETER               :: nmonth = 12
     340    INTEGER, PARAMETER               :: nday = 7
     341    INTEGER, PARAMETER               :: nhour = 24
     342
     343    frac_day= (dd-1)/nday    !> indicates the week of the month, supposing the month starts on monday
     344
     345   ! 1:7      1:31  7                   (0:30)/7
     346    weekday = dd-( nday * (INT( CEILING( frac_day ) ) ) ) ! for now we let the year start on Monday.
     347
     348    !TBD: set weekday correct based on date
     349    index_mm = mo
     350    index_dd = nmonth + weekday  !> Index of the days in the mdh mode (13:20)
     351
     352    SELECT CASE(TRIM(daytype_mdh))
     353
     354                   CASE ("workday")
     355               
     356                   index_hh = nmonth+ nday + hh
     357
     358                   CASE ("weekend")
     359               
     360                   index_hh = nmonth+ nday + nhour + hh
     361
     362                   CASE ("holiday")
     363               
     364                   index_hh = nmonth+ nday + 2*nhour + hh
     365    END SELECT
     366
     367
     368 END SUBROUTINE time_mdh_indices
     369
     370!------------------------------------------------------------------------------!
     371! Description:
     372! ------------
     373!> This routine determines the time factor index in the HOURLY case of the DEFAULT emissions mode.
     374!------------------------------------------------------------------------------!
     375
     376 SUBROUTINE time_hour_indices(mo,dd,hh,index_hh)
     377
     378    USE indices
     379
     380    IMPLICIT NONE
     381
     382    !> IN/OUTPUT
     383    INTEGER, INTENT(INOUT)              :: mo          !> Month
     384    INTEGER, INTENT(INOUT)              :: hh          !> Hour
     385    INTEGER, INTENT(INOUT)              :: dd          !> Day
     386    INTEGER, INTENT(INOUT)              :: index_hh    !> Index Hour
     387 
     388    !> Additional Variables for calculateing indices
     389    INTEGER                             :: index_mm    !> Index Month
     390    INTEGER                             :: index_dd    !> Index Day
     391    INTEGER                             :: i_mon       !> Index for going through the different months
     392    INTEGER                             :: sum_dd      !> Sum days
     393
     394    !> CONSTANTS
     395    INTEGER, PARAMETER                  :: nmonth=12
     396    INTEGER, PARAMETER                  :: nday = 7
     397    INTEGER, PARAMETER                  :: nhour = 24
     398    INTEGER, PARAMETER, DIMENSION(12)   :: days = (/31,28,31,30,31,30,31,31,30,31,30,31/) ! no leap year
     399   
     400
     401    index_mm = mo-1
     402    index_dd = dd-1
     403    sum_dd=0
     404
     405    IF (mo == 1) THEN
     406
     407       index_hh=(index_dd*nhour)+hh
     408
     409    ELSE
     410
     411        DO i_mon=1,index_mm
    95412       
    96 
    97     END SUBROUTINE calc_date_and_time
    98 
    99 
    100 
    101  END MODULE date_and_time_mod
     413            sum_dd=sum_dd+days(i_mon)
     414
     415        ENDDO
     416     
     417    index_hh=(sum_dd*nhour)+(index_dd*nhour)+(hh)
     418
     419    ENDIF
     420 
     421
     422 END SUBROUTINE time_hour_indices
     423
     424
     425END MODULE date_and_time_mod

  • palm/trunk/SOURCE/flow_statistics.f90

    r3294 r3298  
    1919!
    2020! Current revisions:
    21 ! -----------------
     21! ------------------
    2222!
    2323!
    … …  
    2525! -----------------
    2626! $Id$
     27! - Minor formatting (kanani)
     28! - Added  .AND. max_pr_cs > 0 before MPI_ALLREDUCE call (forkel)
     29! - Data arrays, sums, sums_l, hom, hom_sum updated for chem species (basit)
     30! - Directives of parallelism added for chemistry (basit)
     31! - Call for chem_statistics added (basit)
     32!
     33! 3294 2018-10-01 02:37:10Z raasch
    2734! ocean renamed ocean_mode
    2835!
    … …  
    283290!------------------------------------------------------------------------------!
    284291 SUBROUTINE flow_statistics
    285  
     292
    286293
    287294    USE arrays_3d,                                                             &
    … …  
    291298               sa, u, ug, v, vg, vpt, w, w_subs, waterflux_output_conversion,  &
    292299               zw, d_exner
    293        
     300
    294301    USE basic_constants_and_equations_mod,                                     &
    295         ONLY:   g, lv_d_cp
    296        
     302        ONLY:  g, lv_d_cp
     303
     304
     305    USE chem_modules,                                                          &
     306        ONLY:  max_pr_cs
     307
     308    USE chemistry_model_mod,                                                   &
     309        ONLY:  chem_species, chem_statistics
     310
    297311    USE control_parameters,                                                    &
    298         ONLY:   average_count_pr, cloud_droplets, do_sum,                      &
     312        ONLY:   air_chemistry, average_count_pr, cloud_droplets, do_sum,       &
    299313                dt_3d, humidity, initializing_actions, land_surface,           &
    300314                large_scale_forcing, large_scale_subsidence, max_pr_user,      &
    … …  
    303317                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
    304318                ws_scheme_mom, ws_scheme_sca
    305        
     319
    306320    USE cpulog,                                                                &
    307321        ONLY:   cpu_log, log_point
    308        
     322
    309323    USE grid_variables,                                                        &
    310324        ONLY:   ddx, ddy
    … …  
    17801794       ENDIF
    17811795!
     1796!--    Calculate the chemistry module profiles
     1797       IF ( air_chemistry ) THEN
     1798          CALL chem_statistics( 'profiles', sr, tn )
     1799       ENDIF
     1800!
    17821801!--    Calculate the user-defined profiles
    17831802       CALL user_statistics( 'profiles', sr, tn )
    … …  
    17971816                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
    17981817             ENDIF
     1818
     1819             IF ( air_chemistry )  THEN
     1820                IF ( max_pr_cs > 0 )  THEN                                 
     1821                     sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+ max_pr_cs,0) =          &
     1822                               sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs,0) + &
     1823                               sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs,i)
     1824
     1825                ENDIF
     1826             ENDIF
    17991827          ENDDO
    18001828       ENDIF
    … …  
    18111839                              MPI_REAL, MPI_SUM, comm2d, ierr )
    18121840       ENDIF
     1841
     1842       IF ( air_chemistry .AND. max_pr_cs > 0 )  THEN
     1843          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     1844             CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+1,0), sums(nzb,pr_palm+1), &
     1845                                 nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
     1846       ENDIF
     1847
    18131848#else
    18141849       sums = sums_l(:,:,0)
    … …  
    18751910                                    ngp_2dh_s_inner(k,sr)
    18761911          ENDDO
     1912       ENDIF
     1913
     1914       IF ( air_chemistry ) THEN
     1915          IF ( max_pr_cs > 0 )  THEN                 
     1916             DO k = nzb, nzt+1
     1917                sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) = &
     1918                                 sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) / &
     1919                                 ngp_2dh_s_inner(k,sr)
     1920             ENDDO
     1921          ENDIF 
    18771922       ENDIF
    18781923
    … …  
    20182063       ENDIF
    20192064
     2065       IF ( air_chemistry )  THEN
     2066          IF ( max_pr_cs > 0 )  THEN    ! chem_spcs profiles     
     2067             hom(:, 1, pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs, sr) = &
     2068                               sums(:, pr_palm+max_pr_user+1:pr_palm+max_pr_user+max_pr_cs)
     2069          ENDIF
     2070       ENDIF
    20202071!
    20212072!--    Determine the boundary layer height using two different schemes.

  • palm/trunk/SOURCE/header.f90

    r3294 r3298  
    2525! -----------------
    2626! $Id$
     27! Minor formatting (kanani)
     28! Add chemistry header (basit)
     29!
     30! 3294 2018-10-01 02:37:10Z raasch
    2731! changes concerning modularization of ocean option
    2832!
    … …  
    391395        ONLY:  bulk_cloud_model, bcm_header
    392396
     397    USE chemistry_model_mod,                                                   &
     398        ONLY:  chem_header
     399
    393400    USE control_parameters
    394401
    395402    USE cpulog,                                                                &
    396403        ONLY:  log_point_s
    397            
     404
    398405    USE date_and_time_mod,                                                     &
    399406        ONLY:  day_of_year_init, time_utc_init
    … …  
    11231130    IF ( gust_module_enabled )  CALL gust_header ( io )
    11241131
     1132    IF ( air_chemistry )  CALL chem_header ( io )
    11251133!
    11261134!-- Boundary conditions

  • palm/trunk/SOURCE/init_3d_model.f90

    r3294 r3298  
    2525! -----------------
    2626! $Id$
     27! Minor formatting (kanani)
     28! Added CALL to the chem_emissions module (Russo)
     29!
     30! 3294 2018-10-01 02:37:10Z raasch
    2731! changes concerning modularization of ocean option
    2832!
    … …  
    501505!------------------------------------------------------------------------------!
    502506 SUBROUTINE init_3d_model
    503  
     507
    504508
    505509    USE advec_ws
    … …  
    514518        ONLY:  bulk_cloud_model, bcm_init, bcm_init_arrays
    515519
    516     USE chemistry_model_mod,                                                   &
    517         ONLY:  chem_emissions
     520    USE chem_emissions_mod,                                                    &
     521        ONLY:  chem_emissions_init
     522
     523    USE chem_modules,                                                          &
     524        ONLY:  do_emis, max_pr_cs, nspec_out
    518525
    519526    USE control_parameters
    520    
     527
    521528    USE flight_mod,                                                            &
    522529        ONLY:  flight_init
    523    
     530
    524531    USE grid_variables,                                                        &
    525532        ONLY:  dx, dy, ddx2_mg, ddy2_mg
    … …  
    551558
    552559    USE netcdf_data_input_mod,                                                 &
    553         ONLY:  init_3d, netcdf_data_input_init_3d
     560        ONLY:  chem_emis, chem_emis_att, init_3d,                              &
     561               netcdf_data_input_init_3d, netcdf_data_input_interpolate
    554562
    555563    USE ocean_mod,                                                             &
    556564        ONLY:  ocean_init, ocean_init_arrays
    557    
     565
    558566    USE particle_attributes,                                                   &
    559567        ONLY:  particle_advection
    560    
     568
    561569    USE pegrid
    562    
     570
    563571    USE plant_canopy_model_mod,                                                &
    564572        ONLY:  pcm_init
    … …  
    22612269!
    22622270!-- If required, set chemical emissions
    2263 !-- (todo(FK): This should later on be CALLed time-dependently in init_3d_model)
    2264     IF ( air_chemistry )  THEN
    2265        CALL chem_emissions
    2266     ENDIF
    2267 
     2271!-- Initialize values of cssws according to chemistry emission values
     2272    IF ( air_chemistry  .AND.  do_emis )  THEN
     2273       CALL chem_emissions_init( chem_emis_att, chem_emis, nspec_out )
     2274    ENDIF
    22682275!
    22692276!-- Initialize radiation processes

  • palm/trunk/SOURCE/modules.f90

    r3294 r3298  
    2525! -----------------
    2626! $Id$
     27! Minor formatting/clean-up (kanani)
     28! Added some variables for time treatment (Russo)
     29!
     30! 3294 2018-10-01 02:37:10Z raasch
    2731! ocean renamed ocean_mode
    2832!
    … …  
    14011405    REAL(wp) ::  cos_alpha_surface                             !< cosine of alpha_surface
    14021406    REAL(wp) ::  coupling_start_time = 0.0_wp                  !< namelist parameter
     1407    REAL(wp) ::  days_since_reference_point = 0.0_wp           !< days after atmosphere-ocean coupling has been activated,
     1408                                                               !< or after spinup phase of LSM has been finished
    14031409    REAL(wp) ::  disturbance_amplitude = 0.25_wp               !< namelist parameter
    14041410    REAL(wp) ::  disturbance_energy_limit = 0.01_wp            !< namelist parameter

  • palm/trunk/SOURCE/netcdf_data_input_mod.f90

    r3257 r3298  
    2525! -----------------
    2626! $Id$
     27! Modified EXPERT mode to PRE-PROCESSED mode (Russo)
     28! Introduced Chemistry static netcdf file (Russo)
     29! Added the routine for reading-in netcdf data for chemistry (Russo)
     30! Added routines to get_variable interface specific for chemistry files (Russo)
     31!
     32! 3257 2018-09-17 17:11:46Z suehring
    2733! Adjust checks for building_type and building_id, which is necessary after
    2834! topography filtering (building_type and id can be modified by the filtering).
    … …  
    189195!> Modulue contains routines to input data according to Palm input data
    190196!> standart using dynamic and static input files.
    191 !>
     197!> @todo - Review Reading of netcdf files for chemistry
    192198!> @todo - Order input alphabetically
    193199!> @todo - Revise error messages and error numbers
    … …  
    341347
    342348    END TYPE init_type
     349
     350!-- Data type for the general information of chemistry emissions, do not dependent on the particular chemical species
     351    TYPE chem_emis_att_type
     352
     353       !-DIMENSIONS
     354       INTEGER(iwp)                                 :: nspec                     !< number of chem species for which emission values are provided
     355       INTEGER(iwp)                                 :: ncat                      !< number of emission categories
     356       INTEGER(iwp)                                 :: nvoc                      !< number of VOCs components
     357       INTEGER(iwp)                                 :: npm                       !< number of PMs components
     358       INTEGER(iwp)                                 :: nnox=2                    !< number of NOx components: NO and NO2
     359       INTEGER(iwp)                                 :: nsox=2                    !< number of SOx components: SO and SO4
     360       INTEGER(iwp)                                 :: nhoursyear                !< number of hours of a specific year in the HOURLY mode
     361                                                                                 !  of the default mode
     362       INTEGER(iwp)                                 :: nmonthdayhour             !< number of month days and hours in the MDH mode
     363                                                                                 !  of the default mode
     364       INTEGER(iwp)                                 :: dt_emission               !< Number of emissions timesteps for one year
     365                                                                                 !  in the pre-processed emissions case
     366       !-- 1d emission input variables
     367       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: pm_name                   !< Names of PMs components
     368       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: cat_name                  !< Emission categories names
     369       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: species_name              !< Names of emission chemical species
     370       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: voc_name                  !< Names of VOCs components
     371       CHARACTER (LEN=25)                           :: units                     !< Units
     372
     373       INTEGER(iwp)                                 :: i_hour                    !< indices for assigning the emission values at different timesteps
     374       INTEGER(iwp),ALLOCATABLE, DIMENSION(:)       :: cat_index                 !< Index of emission categories
     375       INTEGER(iwp),ALLOCATABLE, DIMENSION(:)       :: species_index             !< Index of emission chem species
     376
     377       REAL(wp), DIMENSION(24)                      :: par_emis_time_factor      !< time factors
     378                                                                                 !  for the parameterized mode: these are fixed for each hour
     379                                                                                 !  of a single day.
     380       REAL(wp),ALLOCATABLE, DIMENSION(:)           :: xm                        !< Molecular masses of emission chem species
     381
     382       !-- 2d emission input variables
     383       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: hourly_emis_time_factor   !< Time factors for HOURLY emissions (DEFAULT mode)
     384       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: mdh_emis_time_factor      !< Time factors for MDH emissions (DEFAULT mode)
     385       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: nox_comp                  !< Composition of NO and NO2
     386       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: sox_comp                  !< Composition of SO2 and SO4
     387       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: voc_comp                  !< Composition of different VOC components (number not fixed)
     388
     389       !-- 3d emission input variables
     390       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)       :: pm_comp                   !< Composition of different PMs components (number not fixed)
     391 
     392    END TYPE chem_emis_att_type
     393
     394
     395!-- Data type for the values of chemistry emissions ERUSSO
     396    TYPE chem_emis_val_type
     397
     398       !REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: stack_height              !< stack height
     399
     400       !-- 3d emission input variables
     401       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)     :: default_emission_data     !< Input Values emissions DEFAULT mode
     402
     403       !-- 4d emission input variables
     404       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:,:)   :: preproc_emission_data      !< Input Values emissions PRE-PROCESSED mode
     405
     406    END TYPE chem_emis_val_type
    343407
    344408!
    … …  
    546610    TYPE(pars)  ::  water_pars_f               !< input variable for water parameters
    547611
     612    TYPE(chem_emis_att_type)                             ::  chem_emis_att    !< Input Information of Chemistry Emission Data from netcdf 
     613    TYPE(chem_emis_val_type), ALLOCATABLE, DIMENSION(:)  ::  chem_emis        !< Input Chemistry Emission Data from netcdf 
    548614
    549615    CHARACTER(LEN=3)  ::  char_lod  = 'lod'         !< name of level-of-detail attribute in NetCDF file
    … …  
    553619    CHARACTER(LEN=100) ::  input_file_static  = 'PIDS_STATIC'  !< Name of file which comprises static input data
    554620    CHARACTER(LEN=100) ::  input_file_dynamic = 'PIDS_DYNAMIC' !< Name of file which comprises dynamic input data
     621    CHARACTER(LEN=100) ::  input_file_chem    = 'PIDS_CHEM'    !< Name of file which comprises chemistry input data
     622
     623    CHARACTER (LEN=25), ALLOCATABLE, DIMENSION(:)    :: string_values  !< output of string variables read from netcdf input files
     624 
     625    INTEGER(iwp)                                     :: id_emis        !< NetCDF id of input file for chemistry emissions: TBD: It has to be removed
    555626
    556627    INTEGER(iwp) ::  nc_stat         !< return value of nf90 function call
    … …  
    558629    LOGICAL ::  input_pids_static  = .FALSE.   !< Flag indicating whether Palm-input-data-standard file containing static information exists
    559630    LOGICAL ::  input_pids_dynamic = .FALSE.   !< Flag indicating whether Palm-input-data-standard file containing dynamic information exists
     631    LOGICAL ::  input_pids_chem = .FALSE.      !< Flag indicating whether Palm-input-data-standard file containing chemistry information exists
    560632
    561633    LOGICAL ::  collective_read = .FALSE.      !< Enable NetCDF collective read
    … …  
    581653       MODULE PROCEDURE netcdf_data_input_check_static
    582654    END INTERFACE netcdf_data_input_check_static
     655
     656    INTERFACE netcdf_data_input_chemistry_data                       
     657       MODULE PROCEDURE netcdf_data_input_chemistry_data
     658    END INTERFACE netcdf_data_input_chemistry_data
    583659
    584660    INTERFACE netcdf_data_input_inquire_file
    … …  
    611687
    612688    INTERFACE get_variable
     689       MODULE PROCEDURE get_variable_string
    613690       MODULE PROCEDURE get_variable_1d_int
    614691       MODULE PROCEDURE get_variable_1d_real
    … …  
    619696       MODULE PROCEDURE get_variable_3d_real
    620697       MODULE PROCEDURE get_variable_3d_real_dynamic
     698       MODULE PROCEDURE get_variable_4d_to_3d_real
    621699       MODULE PROCEDURE get_variable_4d_real
     700       MODULE PROCEDURE get_variable_5d_to_4d_real
    622701    END INTERFACE get_variable
    623702
    … …  
    636715!-- Public variables
    637716    PUBLIC albedo_pars_f, albedo_type_f, basal_area_density_f, buildings_f,    &
    638            building_id_f, building_pars_f, building_type_f, init_3d,           &
    639            init_model, input_file_static, input_pids_static,                   &
     717           building_id_f, building_pars_f, building_type_f,                    &
     718           chem_emis, chem_emis_att, chem_emis_att_type, chem_emis_val_type,   &
     719           init_3d, init_model, input_file_static, input_pids_static,          &
    640720           input_pids_dynamic, leaf_area_density_f, nest_offl,                 &
    641721           pavement_pars_f, pavement_subsurface_pars_f, pavement_type_f,       &
    … …  
    648728!-- Public subroutines
    649729    PUBLIC netcdf_data_input_check_dynamic, netcdf_data_input_check_static,    &
    650            netcdf_data_input_inquire_file,                                     &
     730           netcdf_data_input_chemistry_data, netcdf_data_input_inquire_file,                                     &
    651731           netcdf_data_input_init, netcdf_data_input_init_lsm,                 &
    652732           netcdf_data_input_init_3d,                                          &
    … …  
    674754       INQUIRE( FILE = TRIM( input_file_dynamic ) // TRIM( coupling_char ),    &
    675755                EXIST = input_pids_dynamic )
     756       INQUIRE( FILE = TRIM( input_file_chem ) // TRIM( coupling_char ),    &
     757                EXIST = input_pids_chem )
     758
    676759#endif
    677760
    … …  
    756839
    757840    END SUBROUTINE netcdf_data_input_init
     841
     842!------------------------------------------------------------------------------!
     843! Description:
     844! ------------
     845!> Reads Chemistry NETCDF Input data, such as emission values, emission species, etc. .
     846!------------------------------------------------------------------------------!
     847    SUBROUTINE netcdf_data_input_chemistry_data(emt_att,emt)
     848
     849       USE chem_modules,                                       &
     850           ONLY:  do_emis, mode_emis, time_fac_type,           &
     851                  surface_csflux_name
     852
     853       USE control_parameters,                                 &
     854           ONLY:  message_string
     855
     856       USE indices,                                            &
     857           ONLY:  nz, nx, ny, nxl, nxr, nys, nyn, nzb, nzt
     858
     859       IMPLICIT NONE
     860
     861       TYPE(chem_emis_att_type), INTENT(INOUT)                                        :: emt_att
     862       TYPE(chem_emis_val_type), ALLOCATABLE, DIMENSION(:), INTENT(INOUT)             :: emt
     863   
     864       CHARACTER (LEN=80)                               :: units=''              !< units of chemistry inputs
     865 
     866       INTEGER(iwp)                                     :: ispec                 !< index for number of emission species in input
     867
     868       INTEGER(iwp), ALLOCATABLE, DIMENSION(:)          :: dum_var               !< value of variable read from netcdf input
     869       INTEGER(iwp)                                     :: errno                 !< error number NF90_???? function
     870       INTEGER(iwp)                                     :: id_var                !< variable id
     871!       INTEGER(iwp)                                     :: id_emis               !< NetCDF id of input file
     872       INTEGER(iwp)                                     :: num_vars              !< number of variables in netcdf input file
     873       INTEGER(iwp)                                     :: len_dims,len_dims_2   !< Length of dimensions
     874
     875       INTEGER(iwp)                                     :: max_string_length=25  !< Variable for the maximum length of a string
     876 
     877       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE    :: var_names             !< Name of Variables
     878
     879       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)          :: dum_var_3d            !< variable for storing temporary data of 3-dimensional
     880                                                                                 !  variables read from netcdf for chemistry emissions
     881
     882       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)        :: dum_var_4d            !< variable for storing temporary data of 5-dimensional
     883                                                                                 !< variables read from netcdf for chemistry emissions
     884!--
     885       !> Start the processing of the data
     886       CALL location_message( 'starting allocation of chemistry emissions arrays', .FALSE. )
     887
     888       !> Parameterized mode of the emissions
     889       IF (TRIM(mode_emis)=="PARAMETERIZED" .OR. TRIM(mode_emis)=="parameterized") THEN
     890
     891           ispec=1
     892           emt_att%nspec=0
     893
     894          !number of species
     895           DO  WHILE (TRIM( surface_csflux_name( ispec ) ) /= 'novalue' )
     896
     897             emt_att%nspec=emt_att%nspec+1
     898             ispec=ispec+1
     899
     900           ENDDO
     901
     902          !-- allocate emission values data type arrays
     903          ALLOCATE(emt(emt_att%nspec))
     904
     905          !-- Read EMISSION SPECIES NAMES
     906
     907          !Assign values
     908          ALLOCATE(emt_att%species_name(emt_att%nspec))
     909 
     910         DO ispec=1,emt_att%nspec
     911            emt_att%species_name(ispec) = TRIM(surface_csflux_name(ispec))
     912         ENDDO
     913
     914          !Assign Constant Values of time factors:
     915          emt_att%par_emis_time_factor( : ) = (/ 0.01, 0.01, 0.01, 0.02, 0.03, 0.07, 0.09, 0.09, 0.05, 0.03, 0.03, 0.03, &
     916                                                 0.03, 0.03, 0.03, 0.04, 0.05, 0.09, 0.09, 0.09, 0.04, 0.02, 0.01, 0.01 /)
     917
     918       !> DEFAULT AND PRE-PROCESSED MODE
     919       ELSE
     920
     921#if defined ( __netcdf )       
     922          IF ( .NOT. input_pids_chem )  RETURN
     923
     924          !-- Open file in read-only mode
     925          CALL open_read_file( TRIM( input_file_chem ) //                       &
     926                               TRIM( coupling_char ), id_emis )
     927          !-- inquire number of variables
     928          CALL inquire_num_variables( id_emis, num_vars )
     929
     930          !-- Get General Dimension Lengths: only number of species and number of categories.
     931          !                                  the other dimensions depend on the mode of the emissions or on the presence of specific components
     932          !nspecies
     933          CALL get_dimension_length( id_emis, emt_att%nspec, 'nspecies' )
     934
     935 
     936          !-- Allocate emission values data type arrays
     937          ALLOCATE(emt(1:emt_att%nspec))
     938
     939
     940          !-- Read EMISSION SPECIES NAMES
     941          !Allocate Arrays
     942          ALLOCATE(emt_att%species_name(emt_att%nspec))
     943
     944          !Call get Variable
     945          CALL get_variable( id_emis, 'emission_name', string_values, emt_att%nspec )
     946          emt_att%species_name=string_values
     947          ! If allocated, Deallocate var_string, an array only used for reading-in strings
     948          IF (ALLOCATED(string_values)) DEALLOCATE(string_values)
     949
     950          !-- Read EMISSION SPECIES INDEX
     951          !Allocate Arrays
     952          ALLOCATE(emt_att%species_index(emt_att%nspec))
     953          !Call get Variable
     954          CALL get_variable( id_emis, 'emission_index', emt_att%species_index )
     955
     956
     957          !-- Now the routine has to distinguish between DEFAULT and PRE-PROCESSED chemistry emission modes
     958
     959          IF (TRIM(mode_emis)=="DEFAULT" .OR. TRIM(mode_emis)=="default") THEN
     960 
     961             !number of categories
     962             CALL get_dimension_length( id_emis, emt_att%ncat, 'ncat' )
     963
     964             !-- Read EMISSION CATEGORIES INDEX
     965             !Allocate Arrays
     966             ALLOCATE(emt_att%cat_index(emt_att%ncat))
     967             !Call get Variable
     968             CALL get_variable( id_emis, 'emission_cat_index', emt_att%cat_index )
     969
     970 
     971             DO ispec=1,emt_att%nspec
     972                !-- EMISSION_VOC_NAME (1-DIMENSIONAL)
     973                IF (TRIM(emt_att%species_name(ispec))=="VOC" .OR. TRIM(emt_att%species_name(ispec))=="voc") THEN
     974                   !Allocate Array
     975                   CALL get_dimension_length( id_emis, emt_att%nvoc, 'nvoc' )
     976                   ALLOCATE(emt_att%voc_name(1:emt_att%nvoc))
     977                   !Read-in Variable
     978                   CALL get_variable( id_emis,"emission_voc_name",string_values, emt_att%nvoc)
     979                   emt_att%voc_name=string_values
     980                   IF (ALLOCATED(string_values)) DEALLOCATE(string_values)
     981 
     982                !-- COMPOSITION VOC (2-DIMENSIONAL)
     983                   !Allocate Array
     984                   ALLOCATE(emt_att%voc_comp(1:emt_att%ncat,1:emt_att%nvoc))
     985                   !Read-in Variable
     986!               CALL get_variable(id_emis,"composition_voc",emt%voc_comp,1,1,emt%ncat,emt%nvoc)
     987                   CALL get_variable(id_emis,"composition_voc",emt_att%voc_comp,1,emt_att%ncat,1,emt_att%nvoc)
     988                ENDIF
     989
     990                !-- EMISSION_PM_NAME (1-DIMENSIONAL)
     991                IF (TRIM(emt_att%species_name(ispec))=="PM" .OR. TRIM(emt_att%species_name(ispec))=="pm") THEN
     992                   CALL get_dimension_length( id_emis, emt_att%npm, 'npm' )
     993                   ALLOCATE(emt_att%pm_name(1:emt_att%npm))
     994                   !Read-in Variable
     995                   CALL get_variable( id_emis,"pm_name",string_values, emt_att%npm)
     996                   emt_att%pm_name=string_values
     997                   IF (ALLOCATED(string_values)) DEALLOCATE(string_values)     
     998
     999                !-- COMPOSITION PM (3-DIMENSIONAL)
     1000                   !Allocate
     1001                   len_dims=3  !> number of PMs: PM1, PM2.5 and PM10
     1002                   ALLOCATE(emt_att%pm_comp(1:emt_att%ncat,1:emt_att%npm,1:len_dims))
     1003                   !Read-in Variable
     1004                   CALL get_variable(id_emis,"composition_pm",emt_att%pm_comp,1,emt_att%ncat,1,emt_att%npm,1,len_dims)                   
     1005                ENDIF
     1006
     1007                !-- COMPOSITION_NOX (2-DIMENSIONAL)
     1008                IF (TRIM(emt_att%species_name(ispec))=="NOx" .OR. TRIM(emt_att%species_name(ispec))=="nox") THEN
     1009                   !Allocate array
     1010                   ALLOCATE(emt_att%nox_comp(1:emt_att%ncat,1:emt_att%nnox))
     1011                   !Read-in Variable
     1012                   CALL get_variable(id_emis,"composition_nox",emt_att%nox_comp,1,emt_att%ncat,1,emt_att%nnox)
     1013                ENDIF
     1014
     1015                !-- COMPOSITION-SOX (2-DIMENSIONAL)
     1016                IF (TRIM(emt_att%species_name(ispec))=="SOx" .OR. TRIM(emt_att%species_name(ispec))=="sox") THEN
     1017                   ALLOCATE(emt_att%sox_comp(1:emt_att%ncat,1:emt_att%nsox))
     1018                   !Read-in Variable
     1019                   CALL get_variable(id_emis,"composition_sox",emt_att%sox_comp,1,emt_att%ncat,1,emt_att%nsox)
     1020                ENDIF
     1021             ENDDO !>ispec
     1022
     1023!-- For reading the emission time factors, the distinction between HOUR and MDH data is necessary
     1024     
     1025             !-- EMISSION_TIME_SCALING_FACTORS
     1026                !-- HOUR   
     1027             IF (TRIM(time_fac_type)=="HOUR" .OR. TRIM(time_fac_type)=="hour") THEN
     1028                !-- Allocate Array
     1029                CALL get_dimension_length( id_emis, emt_att%nhoursyear, 'nhoursyear' )
     1030                ALLOCATE(emt_att%hourly_emis_time_factor(1:emt_att%ncat,1:emt_att%nhoursyear))
     1031                !Read-in Variable
     1032                CALL get_variable(id_emis,"emission_time_factors",emt_att%hourly_emis_time_factor,1,emt_att%ncat,1,emt_att%nhoursyear)
     1033
     1034                !-- MDH
     1035             ELSE IF (TRIM(time_fac_type) == "MDH" .OR. TRIM(time_fac_type) == "mdh") THEN
     1036                !-- Allocate Array
     1037                CALL get_dimension_length( id_emis, emt_att%nmonthdayhour, 'nmonthdayhour' )
     1038                ALLOCATE(emt_att%mdh_emis_time_factor(1:emt_att%ncat,1:emt_att%nmonthdayhour))
     1039                !-- Read-in Variable
     1040                CALL get_variable(id_emis,"emission_time_factors",emt_att%mdh_emis_time_factor,1,emt_att%ncat,1,emt_att%nmonthdayhour)
     1041
     1042             ELSE
     1043
     1044             message_string = 'We are in the DEFAULT chemistry emissions mode: '            //          &
     1045                              '     !no time-factor type specified!'                        //          &
     1046                              'Please, specify the value of time_fac_type:'                 //          &
     1047                              '         either "MDH" or "HOUR"'                 
     1048             CALL message( 'netcdf_data_input_chemistry_data', 'CM0200', 2, 2, 0, 6, 0 )
     1049 
     1050
     1051             ENDIF
     1052
     1053             !-- Finally read-in the emission values and their units (DEFAULT mode)
     1054
     1055             DO ispec=1,emt_att%nspec
     1056
     1057                IF ( .NOT. ALLOCATED( emt(ispec)%default_emission_data ) ) ALLOCATE(emt(ispec)%default_emission_data(1:emt_att%ncat,1:ny+1,1:nx+1))
     1058
     1059                ALLOCATE(dum_var_3d(1:emt_att%ncat,nys+1:nyn+1,nxl+1:nxr+1))
     1060
     1061                CALL get_variable(id_emis,"emission_values",dum_var_3d,ispec,1,emt_att%ncat,nys,nyn,nxl,nxr)         
     1062
     1063                emt(ispec)%default_emission_data(:,nys+1:nyn+1,nxl+1:nxr+1)=dum_var_3d(1:emt_att%ncat,nys+1:nyn+1,nxl+1:nxr+1)
     1064
     1065                DEALLOCATE (dum_var_3d)
     1066
     1067             ENDDO
     1068
     1069             !-- UNITS
     1070             CALL get_attribute(id_emis,"units",emt_att%units,.FALSE.,"emission_values")
     1071
     1072
     1073          !-- PRE-PROCESSED MODE --
     1074
     1075          ELSE IF (TRIM(mode_emis)=="PRE-PROCESSED" .OR. TRIM(mode_emis)=="pre-processed") THEN
     1076          !-- In the PRE-PROCESSED mode, only the VOC names, the VOC_composition, the emission values and their units remain to be read at this point
     1077
     1078             DO ispec=1,emt_att%nspec
     1079
     1080             !-- EMISSION_VOC_NAME (1-DIMENSIONAL)
     1081                IF (TRIM(emt_att%species_name(ispec))=="VOC" .OR. TRIM(emt_att%species_name(ispec))=="voc") THEN
     1082                   !Allocate Array
     1083                   CALL get_dimension_length( id_emis, emt_att%nvoc, 'nvoc' )
     1084                   ALLOCATE(emt_att%voc_name(1:emt_att%nvoc))
     1085                   !Read-in Variable
     1086                   CALL get_variable( id_emis,"emission_voc_name",string_values, emt_att%nvoc)
     1087                   emt_att%voc_name=string_values
     1088                   IF (ALLOCATED(string_values)) DEALLOCATE(string_values)
     1089 
     1090             !-- COMPOSITION VOC (2-DIMENSIONAL)
     1091                   !Allocate Array
     1092                   ALLOCATE(emt_att%voc_comp(1:emt_att%ncat,1:emt_att%nvoc))
     1093                   !Read-in Variable
     1094                   CALL get_variable(id_emis,"composition_voc",emt_att%voc_comp,1,emt_att%ncat,1,emt_att%nvoc)
     1095                ENDIF
     1096 
     1097             ENDDO !> ispec
     1098
     1099             !-- EMISSION_VALUES (4-DIMENSIONAL)
     1100             !Calculate temporal dimension length
     1101             CALL get_dimension_length( id_emis, emt_att%dt_emission, 'time' )   
     1102         
     1103
     1104             DO ispec=1,emt_att%nspec
     1105
     1106                !Allocation for the entire domain has to be done only for the first processor between all the subdomains     
     1107                IF ( .NOT. ALLOCATED( emt(ispec)%preproc_emission_data ) ) ALLOCATE(emt(ispec)%preproc_emission_data(emt_att%dt_emission,1,1:ny+1,1:nx+1))
     1108
     1109                !> allocate variable where to pass emission values read from netcdf
     1110                ALLOCATE(dum_var_4d(1:emt_att%dt_emission,1,nys+1:nyn+1,nxl+1:nxr+1))
     1111
     1112                !Read-in Variable
     1113                CALL get_variable(id_emis,"emission_values",dum_var_4d,ispec,1,emt_att%dt_emission,1,1,nys,nyn,nxl,nxr)         
     1114
     1115     
     1116                emt(ispec)%preproc_emission_data(:,1,nys+1:nyn+1,nxl+1:nxr+1)=dum_var_4d(1:emt_att%dt_emission,1,nys+1:nyn+1,nxl+1:nxr+1)
     1117
     1118                DEALLOCATE ( dum_var_4d )
     1119
     1120             ENDDO
     1121
     1122             !-- UNITS
     1123             CALL get_attribute(id_emis,"units",emt_att%units,.FALSE.,"emission_values")
     1124       
     1125          ENDIF
     1126
     1127       CALL close_input_file( id_emis )
     1128
     1129#endif
     1130       ENDIF
     1131
     1132    END SUBROUTINE netcdf_data_input_chemistry_data
    7581133
    7591134!------------------------------------------------------------------------------!
    … …  
    40944469! Description:
    40954470! ------------
     4471!> Routine for reading-in a character string from the chem emissions netcdf input file. 
     4472!------------------------------------------------------------------------------!
     4473 
     4474 SUBROUTINE get_variable_string( id, variable_name, var_string, names_number)
     4475#if defined( __netcdf )
     4476
     4477    USE indices
     4478    USE pegrid
     4479
     4480    IMPLICIT NONE
     4481
     4482    CHARACTER (LEN=25), ALLOCATABLE, DIMENSION(:), INTENT(INOUT)  :: var_string
     4483
     4484    CHARACTER(LEN=*)                                              :: variable_name          !> variable name
     4485
     4486    CHARACTER (LEN=1), ALLOCATABLE, DIMENSION(:,:)                :: tmp_var_string         !> variable to be read
     4487
     4488
     4489    INTEGER(iwp), INTENT(IN)                                      :: id                     !> file id
     4490
     4491    INTEGER(iwp), INTENT(IN)                                      :: names_number           !> number of names
     4492
     4493    INTEGER(iwp)                                                  :: id_var                 !> variable id
     4494
     4495    INTEGER(iwp)                                                  :: i,j                    !> index to go through the length of the dimensions
     4496
     4497    INTEGER(iwp)                                                  :: max_string_length=25   !> this is both the maximum length of a name, but also 
     4498                                                                                            ! the number of the components of the first dimensions
     4499                                                                                            ! (rows)
     4500
     4501
     4502    ALLOCATE(tmp_var_string(max_string_length,names_number))
     4503
     4504    ALLOCATE(var_string(names_number))
     4505
     4506    !-- Inquire variable id
     4507    nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
     4508
     4509
     4510    !-- Get variable
     4511    !-- Start cycle over the emission species
     4512    DO i = 1, names_number
     4513       !-- read the first letter of each component
     4514       nc_stat = NF90_GET_VAR( id, id_var, var_string(i), start = (/ 1,i /), &
     4515                                 count = (/ 1,1 /) )
     4516       CALL handle_error( 'get_variable_string', 701 )
     4517
     4518       !-- Start cycle over charachters
     4519       DO j = 1, max_string_length
     4520                       
     4521          !-- read the rest of the components of the name
     4522          nc_stat = NF90_GET_VAR( id, id_var, tmp_var_string(j,i), start = (/ j,i /),&
     4523                                     count = (/ 1,1 /) )
     4524          CALL handle_error( 'get_variable_string', 702 )
     4525
     4526          IF ( iachar(tmp_var_string(j,i) ) == 0 ) THEN
     4527               tmp_var_string(j,i)=''
     4528          ENDIF
     4529
     4530          IF ( j>1 ) THEN
     4531          !-- Concatenate first letter of the name and the others
     4532             var_string(i)=TRIM(var_string(i)) // TRIM(tmp_var_string(j,i))
     4533
     4534          ENDIF
     4535       ENDDO
     4536    ENDDO
     4537
     4538#endif
     4539 END SUBROUTINE get_variable_string
     4540
     4541
     4542!------------------------------------------------------------------------------!
     4543! Description:
     4544! ------------
    40964545!> Reads a 1D integer variable from file.
    40974546!------------------------------------------------------------------------------!
    … …  
    42334682#endif
    42344683       ENDIF
     4684
     4685
     4686       !Temporary solution for reading emission chemistry files: TBD: we should discuss whether remove it or not
     4687       IF ( id==id_emis ) THEN
     4688
     4689          !--    Allocate temporary variable according to memory access on file.
     4690          ALLOCATE( tmp(is:ie,js:je) )
     4691
     4692          !--    Get variable
     4693          nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
     4694                                  start = (/ is,      js /),                  &
     4695                                  count = (/ ie-is+1 , je-js+1 /) )
     4696
     4697          var=tmp
     4698
     4699          CALL handle_error( 'get_variable_2d_real', 530, variable_name ) !TBD: the error number shuld be changed, but since the solution is
     4700                                                                          ! provisory, we give the same as below
     4701 
     4702          DEALLOCATE( tmp )
     4703       
     4704       !>  Original Subroutine part
     4705       ELSE
    42354706!
    42364707!--    Allocate temporary variable according to memory access on file.
    … …  
    42384709!
    42394710!--    Get variable
    4240        nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
     4711          nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
    42414712                               start = (/ is+1,      js+1 /),                  &
    42424713                               count = (/ ie-is + 1, je-js+1 /) )   
    42434714                               
    4244        CALL handle_error( 'get_variable_2d_real', 530, variable_name )
     4715          CALL handle_error( 'get_variable_2d_real', 530, variable_name )
    42454716!
    42464717!--    Resort data. Please note, dimension subscripts of var all start at 1.
    4247        DO  i = is, ie
    4248           DO  j = js, je
    4249              var(j-js+1,i-is+1) = tmp(i,j)
     4718          DO  i = is, ie
     4719             DO  j = js, je
     4720                var(j-js+1,i-is+1) = tmp(i,j)
     4721             ENDDO
    42504722          ENDDO
    4251        ENDDO
    42524723       
    4253        DEALLOCATE( tmp )
    4254 
     4724          DEALLOCATE( tmp )
     4725
     4726       ENDIF
    42554727#endif
    42564728    END SUBROUTINE get_variable_2d_real
    … …  
    46475119#endif
    46485120       ENDIF
     5121
     5122      !Temporary solution for reading emission chemistry files:
     5123       IF ( id==id_emis ) THEN
     5124
     5125          !--    Allocate temporary variable according to memory access on file.
     5126          ALLOCATE( tmp(is:ie,js:je,k1s:k1e,k2s:k2e) )
     5127
     5128          !--    Get variable
     5129          nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
     5130                                  start = (/ is,   js,   k1s+1,   k2s+1 /),                  &
     5131                                  count = (/ ie-is+1 , je-js+1, k1e-k1s+1, k2e-k2s+1 /) )
     5132
     5133          var=tmp
     5134
     5135          CALL handle_error( 'get_variable_4d_real', 535, variable_name )
     5136 
     5137          DEALLOCATE( tmp )
     5138
     5139       !> Original subroutine part
     5140       ELSE
    46495141!
    46505142!--    Allocate temporary variable according to memory access on file.
    … …  
    46525144!
    46535145!--    Get variable
    4654        nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
     5146          nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
    46555147                               start = (/ is+1,    js+1,   k1s+1, k2s+1 /),    &
    46565148                               count = (/ ie-is+1, je-js+1,                    &
    46575149                                          k1e-k1s+1, k2e-k2s+1 /) )
    46585150
    4659        CALL handle_error( 'get_variable_4d_real', 535, variable_name )
     5151          CALL handle_error( 'get_variable_4d_real', 535, variable_name )
    46605152!
    46615153!--    Resort data. Please note, dimension subscripts of var all start at 1.
    4662        DO  i = is, ie
    4663           DO  j = js, je
    4664              DO  k1 = k1s, k1e
    4665                 DO  k2 = k2s, k2e
    4666                    var(k2-k2s+1,k1-k1s+1,j-js+1,i-is+1) = tmp(i,j,k1,k2)
     5154          DO  i = is, ie
     5155             DO  j = js, je
     5156                DO  k1 = k1s, k1e
     5157                   DO  k2 = k2s, k2e
     5158                      var(k2-k2s+1,k1-k1s+1,j-js+1,i-is+1) = tmp(i,j,k1,k2)
     5159                   ENDDO
    46675160                ENDDO
    46685161             ENDDO
    46695162          ENDDO
    4670        ENDDO
    46715163       
    4672        DEALLOCATE( tmp )
     5164          DEALLOCATE( tmp )
     5165       ENDIF
    46735166#endif
    46745167    END SUBROUTINE get_variable_4d_real
    46755168
    4676 
     5169!------------------------------------------------------------------------------!
     5170! Description:
     5171! ------------
     5172!> Reads a 4D float variable from file and store it to a 3-d variable.
     5173!------------------------------------------------------------------------------!
     5174    SUBROUTINE get_variable_4d_to_3d_real( id, variable_name, var, ns, is, ie, js, je,   &
     5175                                           ks, ke )
     5176
     5177       USE indices
     5178       USE pegrid
     5179
     5180       IMPLICIT NONE
     5181
     5182       CHARACTER(LEN=*)              ::  variable_name   !< variable name
     5183
     5184       INTEGER(iwp)                  ::  ns              !< start index for subdomain input along n dimension: ns coincides here with ne, since, we select only one value along the 1st dimension n
     5185
     5186       INTEGER(iwp)                  ::  i               !< index along x direction
     5187       INTEGER(iwp)                  ::  ie              !< end index for subdomain input along x direction
     5188       INTEGER(iwp)                  ::  is              !< start index for subdomain input along x direction
     5189       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
     5190       INTEGER(iwp)                  ::  id_var          !< variable id
     5191       INTEGER(iwp)                  ::  j               !< index along y direction
     5192       INTEGER(iwp)                  ::  je              !< end index for subdomain input along y direction
     5193       INTEGER(iwp)                  ::  js              !< start index for subdomain input along y direction
     5194       INTEGER(iwp)                  ::  k               !< index along any 4th dimension
     5195       INTEGER(iwp)                  ::  ke              !< end index of 4th dimension
     5196       INTEGER(iwp)                  ::  ks              !< start index of 4th dimension
     5197       
     5198       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp !< temporary variable to read data from file according
     5199                                                         !< to its reverse memory access
     5200
     5201       REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) ::  var  !< variable where the read data have to be stored: one dimension is reduced in the process
     5202#if defined( __netcdf )
     5203
     5204!
     5205!--    Inquire variable id
     5206       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var ) 
     5207!
     5208!--    Check for collective read-operation and set respective NetCDF flags if
     5209!--    required.
     5210       IF ( collective_read )  THEN
     5211          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
     5212       ENDIF
     5213
     5214      !Temporary solution for reading emission chemistry files:
     5215       IF ( id==id_emis ) THEN
     5216
     5217          !--    Allocate temporary variable according to memory access on file.
     5218          ALLOCATE( tmp(is:ie,js:je,ks:ke) )
     5219
     5220          !--    Get variable
     5221          nc_stat = NF90_GET_VAR( id, id_var, tmp(is:ie,js:je,ks:ke),                                &
     5222                                  start = (/ ns, is,   js+1,   ks+1 /),                  &
     5223                                  count = (/ 1, ie-is+1 , je-js+1, ke-ks+1 /) )
     5224
     5225          var=tmp(:,:,:)
     5226
     5227          CALL handle_error( 'get_variable_4d_to_3d_real', 536, variable_name )
     5228 
     5229          DEALLOCATE( tmp )
     5230
     5231       ELSE
     5232!
     5233!--    Allocate temporary variable according to memory access on file.
     5234          ALLOCATE( tmp(is:ie,js:je,ks:ke) )
     5235!
     5236!--    Get variable
     5237          nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
     5238                                  start = (/ ns+1, is+1,    js+1,    ks+1 /),           &
     5239                                  count = (/ 1, ie-is+1, je-js+1, ke-ks+1 /) )   
     5240                               
     5241          CALL handle_error( 'get_variable_4d_to_3d_real', 536, variable_name )
     5242!
     5243!--    Resort data. Please note, dimension subscripts of var all start at 1.
     5244          DO  i = is, ie
     5245             DO  j = js, je
     5246                DO  k = ks, ke
     5247                   var(k-ks+1,j-js+1,i-is+1) = tmp(i,j,k)
     5248                ENDDO
     5249             ENDDO
     5250          ENDDO
     5251       
     5252         DEALLOCATE( tmp )
     5253
     5254       ENDIF
     5255#endif
     5256    END SUBROUTINE get_variable_4d_to_3d_real
    46775257
    46785258!------------------------------------------------------------------------------!
    … …  
    47555335                               count = (/ count_1, count_2, count_3 /) )
    47565336
    4757        CALL handle_error( 'get_variable_3d_real_dynamic', 536, variable_name )
     5337       CALL handle_error( 'get_variable_3d_real_dynamic', 537, variable_name )
    47585338!
    47595339!--    Resort data. Please note, dimension subscripts of var all start at 1.
    … …  
    47705350    END SUBROUTINE get_variable_3d_real_dynamic
    47715351
     5352!------------------------------------------------------------------------------!
     5353! Description:
     5354! ------------
     5355!> Reads a 5D float variable from file and store it to a 4-d variable.
     5356!------------------------------------------------------------------------------!
     5357    SUBROUTINE get_variable_5d_to_4d_real( id, variable_name, var,              &
     5358                                           ns, ts, te, is, ie, js, je, ks, ke )
     5359
     5360       USE indices
     5361       USE pegrid
     5362
     5363       IMPLICIT NONE
     5364
     5365       CHARACTER(LEN=*)              ::  variable_name   !< variable name
     5366
     5367       INTEGER(iwp)                  ::  ns              !< start index for subdomain input along n dimension: ns coincides here with ne, since, we select only one value along the 1st dimension n
     5368
     5369       INTEGER(iwp)                  ::  t               !< index along t direction
     5370       INTEGER(iwp)                  ::  te              !< end index for subdomain input along t direction
     5371       INTEGER(iwp)                  ::  ts              !< start index for subdomain input along t direction
     5372
     5373       INTEGER(iwp)                  ::  i               !< index along x direction
     5374       INTEGER(iwp)                  ::  ie              !< end index for subdomain input along x direction
     5375       INTEGER(iwp)                  ::  is              !< start index for subdomain input along x direction
     5376       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
     5377       INTEGER(iwp)                  ::  id_var          !< variable id
     5378       INTEGER(iwp)                  ::  j               !< index along y direction
     5379       INTEGER(iwp)                  ::  je              !< end index for subdomain input along y direction
     5380       INTEGER(iwp)                  ::  js              !< start index for subdomain input along y direction
     5381       INTEGER(iwp)                  ::  k               !< index along any 5th dimension
     5382       INTEGER(iwp)                  ::  ke              !< end index of 5th dimension
     5383       INTEGER(iwp)                  ::  ks              !< start index of 5th dimension
     5384       
     5385       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE   ::  tmp !< temporary variable to read data from file according
     5386                                                           ! to its reverse memory access
     5387       REAL(wp), DIMENSION(:,:,:,:), INTENT(INOUT) ::  var !< variable to be read
     5388#if defined( __netcdf )
     5389!
     5390!--    Inquire variable id
     5391       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var ) 
     5392!
     5393!--    Check for collective read-operation and set respective NetCDF flags if
     5394!--    required.
     5395       IF ( collective_read )  THEN
     5396          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
     5397       ENDIF
     5398
     5399      !Temporary solution for reading emission chemistry files:
     5400       IF ( id==id_emis ) THEN
     5401
     5402          !--    Allocate temporary variable according to memory access on file.
     5403          ALLOCATE( tmp(ts:te,1,js+1:je+1,ks+1:ke+1) )
     5404
     5405          !--    Get variable
     5406          nc_stat = NF90_GET_VAR( id, id_var, tmp(ts:te,1,js+1:je+1,ks+1:ke+1),               &
     5407                                  start = (/ ns, ts,  1,   js+1,   ks+1 /),                  &
     5408                                  count = (/ 1, te-ts+1, 1, je-js+1, ke-ks+1 /) )
     5409
     5410          var=tmp
     5411
     5412          CALL handle_error( 'get_variable_5d_to_4d_real', 538, variable_name )
     5413 
     5414          DEALLOCATE( tmp )
     5415
     5416       !>  Original Subroutine part
     5417       ELSE
     5418!
     5419!--    Allocate temporary variable according to memory access on file.
     5420          ALLOCATE( tmp(ks:ke,js:je,is:is,ts:te) )
     5421!
     5422!--    Get variable
     5423          nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
     5424                                  start = (/ ks+1, js+1, is+1, ts+1, ns /),           &
     5425                                  count = (/ ke-ks+1, je-js+1, ie-is+1, te-ts+1, 1 /) )   
     5426                               
     5427          CALL handle_error( 'get_variable_5d_to_4d_real', 538, variable_name )
     5428!
     5429!--    Resort data. Please note, dimension subscripts of var all start at 1.
     5430
     5431          DO  t = ts, te
     5432             DO  i = is, ie
     5433                DO  j = js, je
     5434                   DO  k = ks, ke
     5435                      var(t-ts+1,i-is+1,j-js+1,k-ks+1) = tmp(k,j,i,t)
     5436                   ENDDO
     5437                ENDDO
     5438             ENDDO
     5439          ENDDO
     5440
     5441         DEALLOCATE( tmp )
     5442
     5443       ENDIF
     5444#endif
     5445    END SUBROUTINE get_variable_5d_to_4d_real
    47725446
    47735447
    … …  
    47895463
    47905464       nc_stat = NF90_INQUIRE( id, NVARIABLES = num_vars )
    4791        CALL handle_error( 'inquire_num_variables', 537 )
     5465       CALL handle_error( 'inquire_num_variables', 539 )
    47925466
    47935467#endif
    … …  
    48165490       ALLOCATE( varids(1:SIZE(var_names)) )
    48175491       nc_stat = NF90_INQ_VARIDS( id, NVARS = num_vars, VARIDS = varids )
    4818        CALL handle_error( 'inquire_variable_names', 538 )
     5492       CALL handle_error( 'inquire_variable_names', 540 )
    48195493
    48205494       DO  i = 1, SIZE(var_names)
    48215495          nc_stat = NF90_INQUIRE_VARIABLE( id, varids(i), NAME = var_names(i) )
    4822           CALL handle_error( 'inquire_variable_names', 538 )
     5496          CALL handle_error( 'inquire_variable_names', 540 )
    48235497       ENDDO
    48245498

  • palm/trunk/SOURCE/palm.f90

    r3274 r3298  
    2525! -----------------
    2626! $Id$
     27! - Minor formatting (kanani)
     28! - Added Call of date_and_time_init (Russo)
     29! - Added Call of calc_date_and_time before call of init_3d where emissions
     30!   are initialized:
     31!   we have to know the time indices to initialize emission values (Russo)
     32! - Added Call of netcdf_data_input_chemistry_data (Russo)
     33!
     34! 3274 2018-09-24 15:42:55Z knoop
    2735! Modularization of all bulk cloud physics code components
    2836!
    … …  
    241249    USE arrays_3d
    242250
     251    USE bulk_cloud_model_mod,                                                  &
     252        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
     253
     254    USE chem_modules,                                                          &
     255        ONLY:  do_emis
     256
    243257    USE chemistry_model_mod,                                                   &
    244         ONLY:  chem_init
    245 
    246     USE chem_photolysis_mod,                                                   &
    247         ONLY:  photolysis_init
     258        ONLY:  chem_check_data_output_pr, chem_init
     259
     260!    USE chem_photolysis_mod,                                                   &
     261!        ONLY:  photolysis_init
    248262
    249263    USE control_parameters,                                                    &
    … …  
    259273        ONLY:  cpu_log, log_point, log_point_s, cpu_statistics
    260274
     275    USE date_and_time_mod,                                                     &
     276        ONLY:  calc_date_and_time, init_date_and_time
     277
    261278    USE indices,                                                               &
    262279        ONLY:  nbgp
    … …  
    264281    USE kinds
    265282
    266     USE bulk_cloud_model_mod,                                                  &
    267         ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
    268 
    269283    USE multi_agent_system_mod,                                                &
    270284        ONLY:  agents_active, mas_last_actions
    271285
    272286    USE netcdf_data_input_mod,                                                 &
    273         ONLY:  netcdf_data_input_inquire_file, netcdf_data_input_init,         &
     287        ONLY:  chem_emis, chem_emis_att, netcdf_data_input_chemistry_data,     &
     288               netcdf_data_input_inquire_file, netcdf_data_input_init,         &
    274289               netcdf_data_input_surface_data, netcdf_data_input_topo
    275290
    … …  
    410425!-- --> Needs to be moved!! What is the dependency about?
    411426! IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
    412     IF ( air_chemistry )  THEN
     427    IF ( air_chemistry  .AND.  do_emis )  THEN
     428
     429       IF ( do_emis ) THEN
     430          CALL netcdf_data_input_chemistry_data(chem_emis_att,chem_emis)
     431       ENDIF
     432
    413433       CALL chem_init
    414        CALL photolysis_init   ! probably also required for restart
     434!       CALL photolysis_init   ! probably also required for restart
     435
     436       CALL init_date_and_time     !initialize the time of chemistry emissions
     437
    415438    ENDIF
    416439! END IF
    … …  
    421444!
    422445!-- Initialize all necessary variables
     446    CALL calc_date_and_time !this is required for chemistry emissions
     447
    423448    CALL init_3d_model
    424449

  • palm/trunk/SOURCE/parin.f90

    r3294 r3298  
    2525! -----------------
    2626! $Id$
     27! Introduced reading of date_init to inipar.(Russo)
     28! Add extra profiles for chemistry (basit)
     29!
     30! 3294 2018-10-01 02:37:10Z raasch
    2731! changes concerning modularization of ocean option
    2832!
    … …  
    421425! ------------
    422426!> This subroutine reads variables controling the run from the NAMELIST files
     427!>
     428!> @todo: Revise max_pr_cs (profiles for chemistry)
    423429!------------------------------------------------------------------------------!
    424430 SUBROUTINE parin
    … …  
    443449
    444450    USE date_and_time_mod,                                                     &
    445         ONLY:  day_of_year_init, time_utc_init
     451        ONLY:  date_init, day_of_year_init, time_utc_init
    446452
    447453    USE dvrp_variables,                                                        &
    … …  
    548554             cycle_mg, damp_level_1d,                                          &
    549555             data_output_during_spinup,                                        &
     556             date_init,                                                        &
    550557             day_of_year_init,                                                 &
    551558             dissipation_1d,                                                   &
    … …  
    620627             cycle_mg, damp_level_1d,                                          &
    621628             data_output_during_spinup,                                        &
     629             date_init,                                                        &
    622630             day_of_year_init,                                                 &
    623631             dissipation_1d,                                                   &
    … …  
    10661074                       ref_state(0:nz+1), sa_init(0:nz+1), ug(0:nz+1),         &
    10671075                       u_init(0:nz+1), v_init(0:nz+1), vg(0:nz+1),             &
    1068                        hom(0:nz+1,2,pr_palm+max_pr_user,0:statistic_regions),  &
    1069                        hom_sum(0:nz+1,pr_palm+max_pr_user,0:statistic_regions) )
     1076                       hom(0:nz+1,2,pr_palm+max_pr_user+max_pr_cs,0:statistic_regions),  &
     1077                       hom_sum(0:nz+1,pr_palm+max_pr_user+max_pr_cs,0:statistic_regions) )
    10701078
    10711079             hom = 0.0_wp

  • palm/trunk/SOURCE/prognostic_equations.f90

    r3294 r3298  
    2525! -----------------
    2626! $Id$
     27! Code added for decycling chemistry (basit)
     28!
     29! 3294 2018-10-01 02:37:10Z raasch
    2730! changes concerning modularization of ocean option
    2831!
    … …  
    335338
    336339    USE chemistry_model_mod,                                                   &
    337         ONLY:  chem_integrate, chem_prognostic_equations,                      &
    338                chem_species, nspec, nvar, spc_names
     340        ONLY:  chem_integrate, chem_species,  chem_prognostic_equations,       &
     341               nspec, nvar, spc_names, chem_boundary_conds
    339342           
    340343    USE chem_modules,                                                          &
    … …  
    343346    USE chem_photolysis_mod,                                                   &
    344347        ONLY:  photolysis_control
     348
     349    USE chem_modules,                                                          &
     350        ONLY:  call_chem_at_all_substeps, chem_gasphase_on, cs_name
    345351
    346352    USE control_parameters,                                                    &
    … …  
    456462
    457463    LOGICAL      ::  loop_start          !<
    458     INTEGER      ::  n, lsp              !< lsp running index for chem spcs
     464    INTEGER(iwp) ::  n
     465    INTEGER(iwp) :: lsp
     466    INTEGER(iwp) :: lsp_usr              !< lsp running index for chem spcs
    459467
    460468
    … …  
    468476!-- concentrations of chemical species                                   
    469477    IF ( air_chemistry )  THEN
     478       lsp_usr = 1
    470479!
    471480!--    If required, calculate photolysis frequencies -
    … …  
    492501!--    Loop over chemical species       
    493502       CALL cpu_log( log_point_s(84), 'chemistry exch-horiz ', 'start' )
    494        DO  n = 1, nspec
    495           CALL exchange_horiz( chem_species(n)%conc, nbgp )     
     503       DO  lsp = 1, nspec
     504          CALL exchange_horiz( chem_species(lsp)%conc, nbgp )   
     505          lsp_usr = 1 
     506          DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
     507             IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) )  THEN
     508
     509                CALL chem_boundary_conds( chem_species(lsp)%conc_p,                                 &
     510                                          chem_species(lsp)%conc_pr_init )
     511             
     512             ENDIF
     513             lsp_usr = lsp_usr +1
     514          ENDDO
     515         
    496516       ENDDO
    497517       CALL cpu_log( log_point_s(84), 'chemistry exch-horiz ', 'stop' )

  • palm/trunk/SOURCE/read_restart_data_mod.f90

    r3294 r3298  
    2525! -----------------
    2626! $Id$
     27! Added chemistry profiles for restart run (basit)
     28!
     29! 3294 2018-10-01 02:37:10Z raasch
    2730! changes concerning modularization of ocean option
    2831!
    … …  
    7578! ------------
    7679!> Reads restart data from restart-file(s) (binary format).
     80!>
     81!> @todo: Revise max_pr_cs (profiles for chemistry)
    7782!------------------------------------------------------------------------------!
    7883 MODULE read_restart_data_mod
    … …  
    8085
    8186    USE control_parameters
     87
     88    USE chem_modules,                                                                              &
     89       ONLY: max_pr_cs
     90
    8291
    8392    IMPLICIT NONE
    … …  
    266275                    v_init(0:nz+1), pt_init(0:nz+1), q_init(0:nz+1),           &
    267276                    ref_state(0:nz+1), s_init(0:nz+1), sa_init(0:nz+1),        &
    268                     hom(0:nz+1,2,pr_palm+max_pr_user,0:statistic_regions),     &
    269                     hom_sum(0:nz+1,pr_palm+max_pr_user,0:statistic_regions) )
     277                    hom(0:nz+1,2,pr_palm+max_pr_user+max_pr_cs,0:statistic_regions),     &
     278                    hom_sum(0:nz+1,pr_palm+max_pr_user+max_pr_cs,0:statistic_regions) )
    270279       ENDIF
    271280

  • palm/trunk/SOURCE/surface_mod.f90

    r3294 r3298  
    2626! -----------------
    2727! $Id$
     28! - Minor formatting/clean-up (kanani)
     29! - Removed initialization of cssws with surface_csflux values. Only the value 0
     30!   is now passed to all the chemical species of the selected mechanism.
     31!   The real inizialization is conducted in init_3d_model_mod.f90 (Russo)
     32! - Added calls and variables for inizialization of chemistry emission
     33!   fluxes (Russo)
     34! - Setting of cssws to zero moved (forkel)
     35! - Set cssws to zero before setting values (forkel)
     36! - Distinguish flux conversion factors for gases and PM (forkel)
     37! - Some questions/todo's on conversion factors (forkel)
     38!
     39! 3294 2018-10-01 02:37:10Z raasch
    2840! changes concerning modularization of ocean option
    2941!
    … …  
    20712083
    20722084             TYPE( surf_type ) :: surf          !< respective surface type
     2085
    20732086!
    20742087!--          Store indices of respective surface element
    … …  
    21732186
    21742187             DO  lsp = 1, nvar
    2175                 IF ( air_chemistry )  surf%css(lsp,num_h) = 0.0_wp
     2188                IF ( air_chemistry )  surf%css(lsp,num_h)   = 0.0_wp
     2189!
     2190!--             Ensure that fluxes of compounds which are not specified in
     2191!--             namelist are all zero --> kanani: revise
     2192                IF ( air_chemistry )  surf%cssws(lsp,num_h) = 0.0_wp
    21762193             ENDDO
    21772194!
    … …  
    22202237                   ELSE
    22212238                      surf%qsws(num_h) = wall_humidityflux(5) *                &
    2222                                              heatflux_input_conversion(k)
     2239                                                waterflux_input_conversion(k)
    22232240                   ENDIF
    22242241                ENDIF
    … …  
    22272244                   IF ( upward_facing )  THEN
    22282245                      IF ( constant_scalarflux )  THEN
    2229                          surf%ssws(num_h) = surface_scalarflux * rho_air_zw(k-1)
     2246                         surf%ssws(num_h) = surface_scalarflux  * rho_air_zw(k-1)
    22302247
    22312248                         IF ( k-1 /= 0 )                                       &
    … …  
    33863403                                  nxr_on_file, nynf, nync, nyn_on_file, nysf,  &
    33873404                                  nysc, nys_on_file, found )
     3405
    33883406
    33893407       IMPLICIT NONE

  • palm/trunk/SOURCE/time_integration.f90

    r3294 r3298  
    2525! -----------------
    2626! $Id$
     27! - Formatting, clean-up, comments (kanani)
     28! - Added CALL to chem_emissions_setup (Russo)
     29! - Code added for decycling chemistry (basit)
     30!
     31! 3294 2018-10-01 02:37:10Z raasch
    2732! changes concerning modularization of ocean option
    2833!
    … …  
    380385        ONLY:  calc_mean_profile
    381386
     387    USE chem_emissions_mod,                                                    &
     388        ONLY:  chem_emissions_setup
     389
     390    USE chem_modules,                                                          &
     391        ONLY:  bc_cs_t_val, call_chem_at_all_substeps, cs_name,                &
     392               constant_csflux, do_emis, nspec, nspec_out
     393
    382394    USE chemistry_model_mod,                                                   &
    383         ONLY:  chem_emissions, chem_species
    384 
    385     USE chem_modules,                                                          &
    386         ONLY:  nspec
     395        ONLY:  chem_boundary_conds, chem_species
    387396
    388397    USE control_parameters,                                                    &
    … …  
    426435        ONLY:  cpu_log, log_point, log_point_s
    427436
     437    USE date_and_time_mod,                                                     &
     438        ONLY:  calc_date_and_time, hour_call_emis, hour_of_year
     439
    428440    USE flight_mod,                                                            &
    429441        ONLY:  flight_measurement
    … …  
    450462               lsf_nesting_offline, lsf_nesting_offline_mass_conservation
    451463
    452     USE netcdf_data_input_mod,                                                 &
    453         ONLY:  nest_offl, netcdf_data_input_lsf
    454 
    455464    USE multi_agent_system_mod,                                                &
    456465        ONLY:  agents_active, multi_agent_system
     466
     467    USE netcdf_data_input_mod,                                                 &
     468        ONLY:  chem_emis, chem_emis_att, nest_offl, netcdf_data_input_lsf
    457469
    458470    USE ocean_mod,                                                             &
    … …  
    525537    IMPLICIT NONE
    526538
    527     CHARACTER (LEN=9) ::  time_to_string          !<
    528 
    529     INTEGER(iwp)      ::  n  !< loop counter for chemistry species
     539    CHARACTER (LEN=9) ::  time_to_string   !<
     540
     541    INTEGER(iwp)      ::  lsp       !<
     542    INTEGER(iwp)      ::  lsp_usr   !<
     543    INTEGER(iwp)      ::  n         !< loop counter for chemistry species
    530544
    531545    REAL(wp) ::  dt_3d_old  !< temporary storage of timestep to be used for
    … …  
    611625           bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
    612626           bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
     627           IF ( air_chemistry )  THEN
     628              DO  lsp = 1, nspec
     629                 bc_cs_t_val = (  chem_species(lsp)%conc_pr_init(nzt+1)       &
     630                                - chem_species(lsp)%conc_pr_init(nzt) )       &
     631                               / dzu(nzt+1)
     632              ENDDO
     633           ENDIF
    613634       ENDIF
    614635!
    … …  
    783804          IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )
    784805          IF ( air_chemistry )  THEN
    785              DO  n = 1, nspec     
    786                 CALL exchange_horiz( chem_species(n)%conc_p, nbgp )
     806             DO  lsp = 1, nspec
     807                CALL exchange_horiz( chem_species(lsp)%conc_p, nbgp )
     808!
     809!--             kanani: Push chem_boundary_conds after CALL boundary_conds
     810                lsp_usr = 1
     811                DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
     812                   IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) )  THEN
     813                      CALL chem_boundary_conds( chem_species(lsp)%conc_p,              &
     814                                                chem_species(lsp)%conc_pr_init )
     815                   ENDIF
     816                   lsp_usr = lsp_usr + 1
     817                ENDDO
    787818             ENDDO
    788819          ENDIF
    … …  
    11221153!
    11231154!--    If required, consider chemical emissions
    1124 !--    (todo (FK): Implement hourly call of emissions, using time_utc from
    1125 !--                data_and_time_mod.f90;
    1126 !--                move the CALL to appropriate location)
    1127        IF ( air_chemistry ) THEN
    1128           CALL chem_emissions
     1155       IF ( air_chemistry  .AND.  do_emis )  THEN
     1156!
     1157!--       Update the time --> kanani: revise location of this CALL
     1158          CALL calc_date_and_time
     1159!
     1160!--       Call emission routine only once an hour
     1161          IF (hour_of_year  .GT.  hour_call_emis )  THEN
     1162             CALL chem_emissions_setup( chem_emis_att, chem_emis, nspec_out )
     1163             hour_call_emis = hour_of_year
     1164          ENDIF
    11291165       ENDIF
    11301166!
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |