Home


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

Changeset 3600


Ignore:
Timestamp:
Dec 4, 2018 1:49:07 PM (6 years ago)
Author:
banzhafs
Message:

chemistry_model_mod code update to comply PALM coding rules

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/chemistry_model_mod.f90

    r3586 r3600  
    2121!
    2222! Current revisions:
    23 ! ------------------
    24 ! 
    25 ! 
     23! -----------------
     24!
     25!
    2626! Former revisions:
    2727! -----------------
    2828! $Id$
     29! Code update to comply PALM coding rules           
     30! Bug fix in par_dir_diff subroutine                 
     31! Small fixes (corrected 'conastant', added 'Unused')
     32!
     33! 3586 2018-11-30 13:20:29Z dom_dwd_user
    2934! Changed character lenth of name in species_def and photols_def to 15
    30 !
    31 !
     35!
    3236! 3570 2018-11-27 17:44:21Z kanani
    3337! resler:
    3438! Break lines at 132 characters
    35 ! 
     39!
    3640! 3543 2018-11-20 17:06:15Z suehring
    3741! Remove tabs
    38 ! 
     42!
    3943! 3542 2018-11-20 17:04:13Z suehring
    4044! working precision added to make code Fortran 2008 conform
    41 ! 
     45!
    4246! 3458 2018-10-30 14:51:23Z kanani
    4347! from chemistry branch r3443, banzhafs, basit:
    4448! replace surf_lsm_h%qv1(m) by q(k,j,i) for mixing ratio in chem_depo
     49!
    4550! bug fix in chem_depo: allow different surface fractions for one
    4651! surface element and set lai to zero for non vegetated surfaces
    … …  
    212217 MODULE chemistry_model_mod
    213218
    214     USE kinds,              ONLY: wp, iwp
    215     USE indices,            ONLY: nz, nzb,nzt,nysg,nyng,nxlg,nxrg,                                 &
    216                                  nys,nyn,nx,nxl,nxr,ny,wall_flags_0
    217     USE pegrid,             ONLY: myid, threads_per_task
    218 
    219    USE bulk_cloud_model_mod,                                               &
    220         ONLY:  bulk_cloud_model
    221 
    222     USE control_parameters, ONLY: bc_lr_cyc, bc_ns_cyc, dt_3d, humidity,                           &
    223                                  initializing_actions, message_string,                             &
    224                                  omega, tsc, intermediate_timestep_count,                          &
    225                                  intermediate_timestep_count_max, max_pr_user,                     &
    226                                  timestep_scheme, use_prescribed_profile_data, ws_scheme_sca 
    227    USE arrays_3d,          ONLY: exner, hyp, pt, q, ql, rdf_sc, tend, zu
    228     USE chem_gasphase_mod,  ONLY: nspec, spc_names, nkppctrl, nmaxfixsteps,                        &
    229                                  t_steps, chem_gasphase_integrate, vl_dim,                         &
    230                                  nvar, nreact,  atol, rtol, nphot, phot_names
    231     USE cpulog,             ONLY: cpu_log, log_point
     219    USE kinds,                                                                                     &
     220         ONLY:  iwp, wp
     221
     222    USE indices,                                                                                   &
     223         ONLY:  nz, nzb, nzt, nysg, nyng, nxlg, nxrg, nys, nyn, nx, nxl, nxr, ny, wall_flags_0
     224
     225    USE pegrid,                                                                                    &
     226         ONLY: myid, threads_per_task
     227
     228    USE bulk_cloud_model_mod,                                                                      &
     229         ONLY:  bulk_cloud_model
     230
     231    USE control_parameters,                                                                        &
     232         ONLY:  bc_lr_cyc, bc_ns_cyc, dt_3d, humidity, initializing_actions, message_string,        &
     233         omega, tsc, intermediate_timestep_count, intermediate_timestep_count_max,           &
     234         max_pr_user, timestep_scheme, use_prescribed_profile_data, ws_scheme_sca         
     235
     236    USE arrays_3d,                                                                                 &
     237         ONLY:  exner, hyp, pt, q, ql, rdf_sc, tend, zu
     238
     239    USE chem_gasphase_mod,                                                                         &
     240         ONLY:  nspec, spc_names, nkppctrl, nmaxfixsteps, t_steps, chem_gasphase_integrate,         &
     241         vl_dim, nvar, nreact,  atol, rtol, nphot, phot_names
     242
     243    USE cpulog,                                                                                    &
     244         ONLY:  cpu_log, log_point
    232245
    233246    USE chem_modules
    234  
     247
    235248    USE statistics
    236249
    237     IMPLICIT   NONE
     250    IMPLICIT NONE
    238251    PRIVATE
    239252    SAVE
    240253
    241     LOGICAL ::  nest_chemistry = .TRUE. !< flag for nesting mode of chemical species, independent on parent or not
    242 !
    243 !- Define chemical variables
     254    LOGICAL ::  nest_chemistry = .TRUE.  !< flag for nesting mode of chemical species, independent on parent or not
     255    !
     256    !- Define chemical variables
    244257    TYPE   species_def
    245        CHARACTER(LEN=15)                                  :: name
    246        CHARACTER(LEN=16)                                  :: unit
    247        REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: conc
    248        REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: conc_av
    249        REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: conc_p
    250        REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: tconc_m
    251        REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:)           :: cssws_av           
    252        REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:)           :: flux_s_cs
    253        REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:)           :: diss_s_cs
    254        REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:)         :: flux_l_cs
    255        REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:)         :: diss_l_cs
    256        REAL(kind=wp),ALLOCATABLE,DIMENSION(:)             :: conc_pr_init
     258       CHARACTER(LEN=15)                            :: name
     259       CHARACTER(LEN=15)                            :: unit
     260       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     :: conc
     261       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     :: conc_av
     262       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     :: conc_p
     263       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     :: tconc_m
     264       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   :: cssws_av           
     265       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   :: flux_s_cs
     266       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:)   :: diss_s_cs
     267       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) :: flux_l_cs
     268       REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) :: diss_l_cs
     269       REAL(kind=wp), ALLOCATABLE, DIMENSION(:)     :: conc_pr_init
    257270    END TYPE species_def
    258271
    259272
    260273    TYPE   photols_def                                                           
    261        CHARACTER(LEN=15)                                  :: name
    262        CHARACTER(LEN=16)                                  :: unit
    263        REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: freq
     274       CHARACTER(LEN=15)                            :: name
     275       CHARACTER(LEN=15)                            :: unit
     276       REAL(kind=wp), POINTER, DIMENSION(:,:,:)     :: freq
    264277    END TYPE photols_def
    265278
    … …  
    269282
    270283
    271     TYPE(species_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC     :: chem_species
    272     TYPE(photols_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC     :: phot_frequen 
    273 
    274     REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: spec_conc_1
    275     REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: spec_conc_2
    276     REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: spec_conc_3
    277     REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: spec_conc_av       
    278     REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: freq_1
    279 
    280     INTEGER,DIMENSION(nkppctrl)                           :: icntrl                            ! Fine tuning kpp
    281     REAL(kind=wp),DIMENSION(nkppctrl)                     :: rcntrl                            ! Fine tuning kpp
    282     LOGICAL :: decycle_chem_lr    = .FALSE.
    283     LOGICAL :: decycle_chem_ns    = .FALSE.
     284    TYPE(species_def), ALLOCATABLE, DIMENSION(:), TARGET, PUBLIC :: chem_species
     285    TYPE(photols_def), ALLOCATABLE, DIMENSION(:), TARGET, PUBLIC :: phot_frequen 
     286
     287    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_1
     288    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_2
     289    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_3
     290    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: spec_conc_av       
     291    REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:,:), TARGET :: freq_1
     292
     293    INTEGER, DIMENSION(nkppctrl)                           ::  icntrl                            !< Fine tuning kpp
     294    REAL(kind=wp), DIMENSION(nkppctrl)                     ::  rcntrl                            !< Fine tuning kpp
     295    LOGICAL ::  decycle_chem_lr    = .FALSE.
     296    LOGICAL ::  decycle_chem_ns    = .FALSE.
    284297    CHARACTER (LEN=20), DIMENSION(4) ::  decycle_method = &
    285                   (/'dirichlet','dirichlet','dirichlet','dirichlet'/)
    286                                  !< Decycling method at horizontal boundaries,
    287                                  !< 1=left, 2=right, 3=south, 4=north
    288                                  !< dirichlet = initial size distribution and
    289                                  !< chemical composition set for the ghost and
    290                                  !< first three layers
    291                                  !< neumann = zero gradient
    292 
    293     REAL(kind=wp), PUBLIC :: cs_time_step = 0._wp
    294     CHARACTER(10), PUBLIC :: photolysis_scheme
    295                                          ! 'constant',
    296                                          ! 'simple' (Simple parameterisation from MCM, Saunders et al., 2003, Atmos. Chem. Phys., 3, 161-180
    297                                          ! 'fastj'  (Wild et al., 2000, J. Atmos. Chem., 37, 245-282) STILL NOT IMPLEMENTED
    298 
    299    
    300     !-----------------------------------------------------------------------
    301     ! Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010)
    302     !
    303     INTEGER(iwp), PARAMETER :: nlu_dep = 15                   ! Number of DEPAC landuse classes (lu's)
    304     INTEGER(iwp), PARAMETER :: ncmp = 10                      ! Number of DEPAC gas components
    305     !------------------------------------------------------------------------
    306     ! DEPAC landuse classes as defined in LOTOS-EUROS model v2.1                             
    307     !
    308     INTEGER(iwp)  ::  ilu_grass              = 1       
    309     INTEGER(iwp)  ::  ilu_arable             = 2       
    310     INTEGER(iwp)  ::  ilu_permanent_crops    = 3         
    311     INTEGER(iwp)  ::  ilu_coniferous_forest  = 4         
    312     INTEGER(iwp)  ::  ilu_deciduous_forest   = 5         
    313     INTEGER(iwp)  ::  ilu_water_sea          = 6       
    314     INTEGER(iwp)  ::  ilu_urban              = 7       
    315     INTEGER(iwp)  ::  ilu_other              = 8 
    316     INTEGER(iwp)  ::  ilu_desert             = 9 
    317     INTEGER(iwp)  ::  ilu_ice                = 10
    318     INTEGER(iwp)  ::  ilu_savanna            = 11
    319     INTEGER(iwp)  ::  ilu_tropical_forest    = 12
    320     INTEGER(iwp)  ::  ilu_water_inland       = 13
    321     INTEGER(iwp)  ::  ilu_mediterrean_scrub  = 14
    322     INTEGER(iwp)  ::  ilu_semi_natural_veg   = 15
    323     !
    324     !--------------------------------------------------------------------------
    325     ! NH3/SO2 ratio regimes:
    326     INTEGER, PARAMETER  ::  iratns_low  = 1       ! low ratio NH3/SO2
    327     INTEGER, PARAMETER  ::  iratns_high = 2       ! high ratio NH3/SO2
    328     INTEGER, PARAMETER  ::  iratns_very_low = 3   ! very low ratio NH3/SO2
    329     ! Default:
    330     INTEGER, PARAMETER  ::  iratns_default = iratns_low
    331     !
    332     ! Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) naar W m-2
    333     REAL(wp), DIMENSION(nlu_dep), PARAMETER :: alpha   = &
    334       (/0.009,0.009, 0.009,0.006,0.006, -999., -999.,0.009,-999.,-999.,0.009,0.006,-999.,0.009,0.008/)*4.57
    335     !
    336     ! Set temperatures per land use for F_temp
    337     REAL(wp), DIMENSION(nlu_dep), PARAMETER :: Tmin = &
    338       (/12.0, 12.0,  12.0,  0.0,  0.0, -999., -999., 12.0, -999., -999., 12.0,  0.0, -999., 12.0,  8.0/)
    339     REAL(wp), DIMENSION(nlu_dep), PARAMETER :: Topt = &
    340       (/26.0, 26.0,  26.0, 18.0, 20.0, -999., -999., 26.0, -999., -999., 26.0, 20.0, -999., 26.0, 24.0 /)
    341     REAL(wp), DIMENSION(nlu_dep), PARAMETER :: Tmax = &
    342       (/40.0, 40.0,  40.0, 36.0, 35.0, -999., -999., 40.0, -999., -999., 40.0, 35.0, -999., 40.0, 39.0 /)
    343     !
    344     ! Set F_min:
    345     REAL(wp), DIMENSION(nlu_dep), PARAMETER :: F_min = &
    346       (/0.01, 0.01,  0.01, 0.1,  0.1,   -999., -999.,0.01, -999.,-999.,0.01,0.1,-999.,0.01, 0.04/)
    347    
    348     ! Set maximal conductance (m/s)
    349     ! (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from  mmol O3/m2/s to m/s
    350     ! Could be refined to a function of T and P. in Jones
    351     REAL(wp), DIMENSION(nlu_dep), PARAMETER :: g_max = &
    352       (/270., 300.,  300., 140., 150.,  -999., -999.,270., -999.,-999.,270., 150.,-999.,300., 422./)/41000
    353     !
    354     ! Set max, min for vapour pressure deficit vpd;
    355     REAL(wp), DIMENSION(nlu_dep), PARAMETER :: vpd_max = &
    356       (/1.3,  0.9,   0.9,  0.5,  1.0,   -999., -999.,1.3,  -999.,-999.,1.3,1.0,  -999.,0.9, 2.8/)
    357     REAL(wp), DIMENSION(nlu_dep), PARAMETER :: vpd_min = &
    358       (/3.0,  2.8,   2.8,  3.0,  3.25,  -999., -999.,3.0,  -999.,-999.,3.0,3.25, -999.,2.8, 4.5/)
    359     !
    360     ! 
    361     !------------------------------------------------------------------------
    362 
    363    
     298         (/'dirichlet','dirichlet','dirichlet','dirichlet'/)
     299                              !< Decycling method at horizontal boundaries,
     300                              !< 1=left, 2=right, 3=south, 4=north
     301                              !< dirichlet = initial size distribution and
     302                              !< chemical composition set for the ghost and
     303                              !< first three layers
     304                              !< neumann = zero gradient
     305
     306    REAL(kind=wp), PUBLIC ::  cs_time_step = 0._wp
     307    CHARACTER(10), PUBLIC ::  photolysis_scheme
     308                              !< 'constant',
     309                              !< 'simple' (Simple parameterisation from MCM, Saunders et al., 2003, Atmos. Chem. Phys., 3, 161-180
     310                              !< 'fastj'  (Wild et al., 2000, J. Atmos. Chem., 37, 245-282) STILL NOT IMPLEMENTED
     311
     312
     313    !
     314    !-- Parameter needed for Deposition calculation using DEPAC model (van Zanten et al., 2010)
     315    !
     316    INTEGER(iwp), PARAMETER ::  nlu_dep = 15                   !< Number of DEPAC landuse classes (lu's)
     317    INTEGER(iwp), PARAMETER ::  ncmp = 10                      !< Number of DEPAC gas components
     318    INTEGER(iwp), PARAMETER ::  nposp = 69                     !< Number of possible species for deposition
     319    !
     320    !-- DEPAC landuse classes as defined in LOTOS-EUROS model v2.1                             
     321    INTEGER(iwp) ::  ilu_grass              = 1       
     322    INTEGER(iwp) ::  ilu_arable             = 2       
     323    INTEGER(iwp) ::  ilu_permanent_crops    = 3         
     324    INTEGER(iwp) ::  ilu_coniferous_forest  = 4         
     325    INTEGER(iwp) ::  ilu_deciduous_forest   = 5         
     326    INTEGER(iwp) ::  ilu_water_sea          = 6       
     327    INTEGER(iwp) ::  ilu_urban              = 7       
     328    INTEGER(iwp) ::  ilu_other              = 8 
     329    INTEGER(iwp) ::  ilu_desert             = 9 
     330    INTEGER(iwp) ::  ilu_ice                = 10
     331    INTEGER(iwp) ::  ilu_savanna            = 11
     332    INTEGER(iwp) ::  ilu_tropical_forest    = 12
     333    INTEGER(iwp) ::  ilu_water_inland       = 13
     334    INTEGER(iwp) ::  ilu_mediterrean_scrub  = 14
     335    INTEGER(iwp) ::  ilu_semi_natural_veg   = 15
     336
     337    !
     338    !-- NH3/SO2 ratio regimes:
     339    INTEGER(iwp), PARAMETER ::  iratns_low      = 1       !< low ratio NH3/SO2
     340    INTEGER(iwp), PARAMETER ::  iratns_high     = 2       !< high ratio NH3/SO2
     341    INTEGER(iwp), PARAMETER ::  iratns_very_low = 3       !< very low ratio NH3/SO2
     342    !
     343    !-- Default:
     344    INTEGER, PARAMETER ::  iratns_default = iratns_low
     345    !
     346    !-- Set alpha for f_light (4.57 is conversion factor from 1./(mumol m-2 s-1) to W m-2
     347    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  alpha   =(/ 0.009, 0.009, 0.009, 0.006, 0.006, -999., -999., 0.009, -999.,      &
     348         -999., 0.009, 0.006, -999., 0.009, 0.008/)*4.57
     349    !
     350    !-- Set temperatures per land use for f_temp
     351    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmin = (/ 12.0, 12.0,  12.0,  0.0,  0.0, -999., -999., 12.0, -999., -999.,      &
     352         12.0,  0.0, -999., 12.0,  8.0/)
     353    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  topt = (/ 26.0, 26.0,  26.0, 18.0, 20.0, -999., -999., 26.0, -999., -999.,      &
     354         26.0, 20.0, -999., 26.0, 24.0 /)
     355    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  tmax = (/ 40.0, 40.0,  40.0, 36.0, 35.0, -999., -999., 40.0, -999., -999.,      &
     356         40.0, 35.0, -999., 40.0, 39.0 /)
     357    !
     358    !-- Set f_min:
     359    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  f_min = (/ 0.01, 0.01, 0.01, 0.1, 0.1, -999., -999., 0.01, -999., -999., 0.01,  &
     360         0.1, -999., 0.01, 0.04/)
     361
     362    !
     363    !-- Set maximal conductance (m/s)
     364    !-- (R T/P) = 1/41000 mmol/m3 is given for 20 deg C to go from  mmol O3/m2/s to m/s
     365    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  g_max = (/ 270., 300., 300., 140., 150., -999., -999., 270., -999., -999.,      &
     366         270., 150., -999., 300., 422./)/41000
     367    !
     368    !-- Set max, min for vapour pressure deficit vpd
     369    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_max = (/1.3, 0.9, 0.9, 0.5, 1.0, -999., -999., 1.3, -999., -999., 1.3,      &
     370         1.0, -999., 0.9, 2.8/)
     371    REAL(wp), DIMENSION(nlu_dep), PARAMETER ::  vpd_min = (/3.0, 2.8, 2.8, 3.0, 3.25, -999., -999., 3.0, -999., -999., 3.0,     &
     372         3.25, -999., 2.8, 4.5/)
     373
     374
    364375    PUBLIC nest_chemistry
    365376    PUBLIC nspec
    … …  
    369380    PUBLIC spec_conc_2 
    370381
    371 !-  Interface section
     382    !   
     383    !-- Interface section
    372384    INTERFACE chem_3d_data_averaging
    373385       MODULE PROCEDURE chem_3d_data_averaging
    … …  
    516528    END INTERFACE get_rb_cell
    517529
    518    
     530
    519531
    520532    PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_check_data_output, &
    521            chem_check_data_output_pr, chem_check_parameters,                    &
    522            chem_data_output_2d, chem_data_output_3d, chem_data_output_mask,     &
    523            chem_define_netcdf_grid, chem_header, chem_init,                     &
    524            chem_init_profiles, chem_integrate, chem_parin,                      &
    525            chem_prognostic_equations, chem_rrd_local,                           &
    526            chem_statistics, chem_swap_timelevel, chem_wrd_local, chem_depo
     533         chem_check_data_output_pr, chem_check_parameters,                    &
     534         chem_data_output_2d, chem_data_output_3d, chem_data_output_mask,     &
     535         chem_define_netcdf_grid, chem_header, chem_init,                     &
     536         chem_init_profiles, chem_integrate, chem_parin,                      &
     537         chem_prognostic_equations, chem_rrd_local,                           &
     538         chem_statistics, chem_swap_timelevel, chem_wrd_local, chem_depo
    527539
    528540
    … …  
    530542 CONTAINS
    531543
    532 !
    533 !------------------------------------------------------------------------------!
    534 !
    535 ! Description:
    536 ! ------------
    537 !> Subroutine for averaging 3D data of chemical species. Due to the fact that
    538 !> the averaged chem arrays are allocated in chem_init, no if-query concerning
    539 !> the allocation is required (in any mode). Attention: If you just specify an
    540 !> averaged output quantity in the _p3dr file during restarts the first output
    541 !> includes the time between the beginning of the restart run and the first
    542 !> output time (not necessarily the whole averaging_interval you have
    543 !> specified in your _p3d/_p3dr file )
    544 !------------------------------------------------------------------------------!
    545 
    546     SUBROUTINE chem_3d_data_averaging ( mode, variable )
    547 
    548        USE control_parameters
    549 
    550        USE indices
    551 
    552        USE kinds
    553 
    554        USE surface_mod,                                                         &
    555           ONLY: surf_def_h, surf_lsm_h, surf_usm_h   
    556  
    557        IMPLICIT NONE
    558  
    559        CHARACTER (LEN=*) ::  mode    !<
    560        CHARACTER (LEN=*) :: variable !<
    561  
    562        LOGICAL      ::  match_def !< flag indicating natural-type surface
    563        LOGICAL      ::  match_lsm !< flag indicating natural-type surface
    564        LOGICAL      ::  match_usm !< flag indicating urban-type surface
    565 
    566        INTEGER(iwp) ::  i                  !< grid index x direction
    567        INTEGER(iwp) ::  j                  !< grid index y direction
    568        INTEGER(iwp) ::  k                  !< grid index z direction
    569        INTEGER(iwp) ::  m                  !< running index surface type
    570        INTEGER(iwp) :: lsp                 !< running index for chem spcs
    571 
    572 
    573        IF ( mode == 'allocate' )  THEN
    574           DO lsp = 1, nspec
    575              IF ( TRIM(variable(1:3)) == 'kc_' .AND. &
    576                   TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
    577                 chem_species(lsp)%conc_av = 0.0_wp
    578              ENDIF
    579           ENDDO
    580 
    581        ELSEIF ( mode == 'sum' )  THEN
    582 
    583           DO lsp = 1, nspec
    584              IF ( TRIM(variable(1:3)) == 'kc_' .AND. &
    585                   TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
    586                 DO  i = nxlg, nxrg
    587                    DO  j = nysg, nyng
    588                       DO  k = nzb, nzt+1
    589                           chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) + &
    590                                                              chem_species(lsp)%conc(k,j,i)
    591                       ENDDO
     544 !
     545 !------------------------------------------------------------------------------!
     546 !
     547 ! Description:
     548 ! ------------
     549 !> Subroutine for averaging 3D data of chemical species. Due to the fact that
     550 !> the averaged chem arrays are allocated in chem_init, no if-query concerning
     551 !> the allocation is required (in any mode). Attention: If you just specify an
     552 !> averaged output quantity in the _p3dr file during restarts the first output
     553 !> includes the time between the beginning of the restart run and the first
     554 !> output time (not necessarily the whole averaging_interval you have
     555 !> specified in your _p3d/_p3dr file )
     556 !------------------------------------------------------------------------------!
     557
     558 SUBROUTINE chem_3d_data_averaging( mode, variable )
     559
     560    USE control_parameters
     561
     562    USE indices
     563
     564    USE kinds
     565
     566    USE surface_mod,                                                         &
     567         ONLY:  surf_def_h, surf_lsm_h, surf_usm_h   
     568
     569    IMPLICIT NONE
     570
     571    CHARACTER (LEN=*) ::  mode    !<
     572    CHARACTER (LEN=*) :: variable !<
     573
     574    LOGICAL      ::  match_def !< flag indicating natural-type surface
     575    LOGICAL      ::  match_lsm !< flag indicating natural-type surface
     576    LOGICAL      ::  match_usm !< flag indicating urban-type surface
     577
     578    INTEGER(iwp) ::  i                  !< grid index x direction
     579    INTEGER(iwp) ::  j                  !< grid index y direction
     580    INTEGER(iwp) ::  k                  !< grid index z direction
     581    INTEGER(iwp) ::  m                  !< running index surface type
     582    INTEGER(iwp) :: lsp                 !< running index for chem spcs
     583
     584
     585    IF ( mode == 'allocate' )  THEN
     586       DO lsp = 1, nspec
     587          IF ( TRIM(variable(1:3)) == 'kc_' .AND. &
     588               TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
     589             chem_species(lsp)%conc_av = 0.0_wp
     590          ENDIF
     591       ENDDO
     592
     593    ELSEIF ( mode == 'sum' )  THEN
     594
     595       DO lsp = 1, nspec
     596          IF ( TRIM(variable(1:3)) == 'kc_' .AND. &
     597               TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
     598             DO  i = nxlg, nxrg
     599                DO  j = nysg, nyng
     600                   DO  k = nzb, nzt+1
     601                      chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) + &
     602                           chem_species(lsp)%conc(k,j,i)
    592603                   ENDDO
    593604                ENDDO
    594              ELSEIF ( TRIM(variable(4:)) == TRIM('cssws*') )        THEN
    595                 DO  i = nxl, nxr
    596                    DO  j = nys, nyn
    597                       match_def = surf_def_h(0)%start_index(j,i) <=               &
    598                                   surf_def_h(0)%end_index(j,i)
    599                       match_lsm = surf_lsm_h%start_index(j,i) <=                  &
    600                                   surf_lsm_h%end_index(j,i)
    601                       match_usm = surf_usm_h%start_index(j,i) <=                  &
    602                                   surf_usm_h%end_index(j,i)
    603 
    604                       IF ( match_def )  THEN
    605                          m = surf_def_h(0)%end_index(j,i)
    606                          chem_species(lsp)%cssws_av(j,i) =                        &
    607                                                 chem_species(lsp)%cssws_av(j,i) + &
    608                                                 surf_def_h(0)%cssws(lsp,m)
    609                       ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
    610                          m = surf_lsm_h%end_index(j,i)
    611                          chem_species(lsp)%cssws_av(j,i) =                        &
    612                                                 chem_species(lsp)%cssws_av(j,i) + &
    613                                                 surf_lsm_h%cssws(lsp,m)
    614                       ELSEIF ( match_usm )  THEN
    615                          m = surf_usm_h%end_index(j,i)
    616                          chem_species(lsp)%cssws_av(j,i) =                        &
    617                                                 chem_species(lsp)%cssws_av(j,i) + &
    618                                                 surf_usm_h%cssws(lsp,m)
    619                       ENDIF
     605             ENDDO
     606          ELSEIF ( TRIM(variable(4:)) == TRIM('cssws*') )        THEN
     607             DO  i = nxl, nxr
     608                DO  j = nys, nyn
     609                   match_def = surf_def_h(0)%start_index(j,i) <=               &
     610                        surf_def_h(0)%end_index(j,i)
     611                   match_lsm = surf_lsm_h%start_index(j,i) <=                  &
     612                        surf_lsm_h%end_index(j,i)
     613                   match_usm = surf_usm_h%start_index(j,i) <=                  &
     614                        surf_usm_h%end_index(j,i)
     615
     616                   IF ( match_def )  THEN
     617                      m = surf_def_h(0)%end_index(j,i)
     618                      chem_species(lsp)%cssws_av(j,i) =                        &
     619                           chem_species(lsp)%cssws_av(j,i) + &
     620                           surf_def_h(0)%cssws(lsp,m)
     621                   ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
     622                      m = surf_lsm_h%end_index(j,i)
     623                      chem_species(lsp)%cssws_av(j,i) =                        &
     624                           chem_species(lsp)%cssws_av(j,i) + &
     625                           surf_lsm_h%cssws(lsp,m)
     626                   ELSEIF ( match_usm )  THEN
     627                      m = surf_usm_h%end_index(j,i)
     628                      chem_species(lsp)%cssws_av(j,i) =                        &
     629                           chem_species(lsp)%cssws_av(j,i) + &
     630                           surf_usm_h%cssws(lsp,m)
     631                   ENDIF
     632                ENDDO
     633             ENDDO
     634          ENDIF
     635       ENDDO
     636
     637    ELSEIF ( mode == 'average' )  THEN
     638       DO lsp = 1, nspec
     639          IF ( TRIM(variable(1:3)) == 'kc_' .AND. &
     640               TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
     641             DO  i = nxlg, nxrg
     642                DO  j = nysg, nyng
     643                   DO  k = nzb, nzt+1
     644                      chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    620645                   ENDDO
    621646                ENDDO
    622              ENDIF
    623           ENDDO
    624  
    625        ELSEIF ( mode == 'average' )  THEN
    626           DO lsp = 1, nspec
    627              IF ( TRIM(variable(1:3)) == 'kc_' .AND. &
    628                   TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
    629                 DO  i = nxlg, nxrg
    630                    DO  j = nysg, nyng
    631                       DO  k = nzb, nzt+1
    632                           chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    633                       ENDDO
    634                    ENDDO
     647             ENDDO
     648
     649          ELSEIF (TRIM(variable(4:)) == TRIM('cssws*') )            THEN
     650             DO i = nxlg, nxrg
     651                DO  j = nysg, nyng
     652                   chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp )
    635653                ENDDO
    636 
    637              ELSEIF (TRIM(variable(4:)) == TRIM('cssws*') )            THEN
    638                 DO i = nxlg, nxrg
    639                    DO  j = nysg, nyng
    640                         chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp )
    641                    ENDDO
    642                 ENDDO
    643                    CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp )                 
    644              ENDIF
    645           ENDDO
    646          
    647        ENDIF     
    648 
    649     END SUBROUTINE chem_3d_data_averaging
     654             ENDDO
     655             CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp )                 
     656          ENDIF
     657       ENDDO
     658
     659    ENDIF
     660
     661 END SUBROUTINE chem_3d_data_averaging
    650662
    651663!   
    … …  
    657669!------------------------------------------------------------------------------!
    658670
    659     SUBROUTINE chem_boundary_conds( mode )                                           
    660                                                                                      
    661        USE control_parameters,                                                    & 
    662           ONLY:  air_chemistry, bc_radiation_l, bc_radiation_n, bc_radiation_r,  &
    663                  bc_radiation_s               
    664        USE indices,                                                               & 
    665           ONLY:  nxl, nxr,  nxlg, nxrg, nyng, nysg, nzt                             
    666                                                                                      
    667 
    668        USE arrays_3d,                                                             &     
    669           ONLY:  dzu                                               
    670        USE surface_mod,                                                           &
    671           ONLY:  bc_h                                                           
    672 
    673        CHARACTER (len=*), INTENT(IN) :: mode
    674        INTEGER(iwp) ::  i                                                            !< grid index x direction.
    675        INTEGER(iwp) ::  j                                                            !< grid index y direction.
    676        INTEGER(iwp) ::  k                                                            !< grid index z direction.
    677        INTEGER(iwp) ::  kb                                                           !< variable to set respective boundary value, depends on facing.
    678        INTEGER(iwp) ::  l                                                            !< running index boundary type, for up- and downward-facing walls.
    679        INTEGER(iwp) ::  m                                                            !< running index surface elements.
    680        INTEGER(iwp) ::  lsp                                                          !< running index for chem spcs.
    681        INTEGER(iwp) ::  lph                                                          !< running index for photolysis frequencies
    682 
    683 
    684        SELECT CASE ( TRIM( mode ) )       
    685           CASE ( 'init' )
    686 
    687              IF ( bc_cs_b == 'dirichlet' )    THEN
    688                 ibc_cs_b = 0
    689              ELSEIF ( bc_cs_b == 'neumann' )  THEN
    690                 ibc_cs_b = 1
    691              ELSE
    692                 message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"' 
    693                 CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 )     !<
    694              ENDIF                                                                 
     671 SUBROUTINE chem_boundary_conds( mode )                                           
     672
     673    USE control_parameters,                                                    & 
     674        ONLY:  air_chemistry, bc_radiation_l, bc_radiation_n, bc_radiation_r,  &
     675               bc_radiation_s               
     676    USE indices,                                                               & 
     677        ONLY:  nxl, nxr,  nxlg, nxrg, nyng, nysg, nzt                             
     678
     679    USE arrays_3d,                                                             &     
     680        ONLY:  dzu                                                
     681
     682    USE surface_mod,                                                           &
     683        ONLY:  bc_h                                                           
     684
     685    CHARACTER (len=*), INTENT(IN) :: mode
     686    INTEGER(iwp) ::  i                            !< grid index x direction.
     687    INTEGER(iwp) ::  j                            !< grid index y direction.
     688    INTEGER(iwp) ::  k                            !< grid index z direction.
     689    INTEGER(iwp) ::  kb                           !< variable to set respective boundary value, depends on facing.
     690    INTEGER(iwp) ::  l                            !< running index boundary type, for up- and downward-facing walls.
     691    INTEGER(iwp) ::  m                            !< running index surface elements.
     692    INTEGER(iwp) ::  lsp                          !< running index for chem spcs.
     693    INTEGER(iwp) ::  lph                          !< running index for photolysis frequencies
     694
     695
     696    SELECT CASE ( TRIM( mode ) )       
     697       CASE ( 'init' )
     698
     699          IF ( bc_cs_b == 'dirichlet' )    THEN
     700             ibc_cs_b = 0
     701          ELSEIF ( bc_cs_b == 'neumann' )  THEN
     702             ibc_cs_b = 1
     703          ELSE
     704             message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"' 
     705             CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 )     !<
     706          ENDIF                                                                 
    695707!
    696708!--          Set Integer flags and check for possible erroneous settings for top
    697709!--          boundary condition.
    698              IF ( bc_cs_t == 'dirichlet' )             THEN
    699                 ibc_cs_t = 0
    700              ELSEIF ( bc_cs_t == 'neumann' )           THEN
    701                 ibc_cs_t = 1
    702              ELSEIF ( bc_cs_t == 'initial_gradient' )  THEN
    703                 ibc_cs_t = 2
    704              ELSEIF ( bc_cs_t == 'nested' )            THEN         
    705                 ibc_cs_t = 3
    706              ELSE
    707                 message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"'     
    708                 CALL message( 'check_parameters', 'CM0430', 1, 2, 0, 6, 0 )
    709              ENDIF
    710 
    711          
    712           CASE ( 'set_bc_bottomtop' )                   
     710          IF ( bc_cs_t == 'dirichlet' )             THEN
     711             ibc_cs_t = 0
     712          ELSEIF ( bc_cs_t == 'neumann' )           THEN
     713             ibc_cs_t = 1
     714          ELSEIF ( bc_cs_t == 'initial_gradient' )  THEN
     715             ibc_cs_t = 2
     716          ELSEIF ( bc_cs_t == 'nested' )            THEN         
     717             ibc_cs_t = 3
     718          ELSE
     719             message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"'     
     720             CALL message( 'check_parameters', 'CM0430', 1, 2, 0, 6, 0 )
     721          ENDIF
     722
     723
     724       CASE ( 'set_bc_bottomtop' )                   
    713725!--          Bottom boundary condtions for chemical species     
    714              DO lsp = 1, nspec                                                     
    715                 IF ( ibc_cs_b == 0 ) THEN                   
    716                    DO l = 0, 1
     726          DO lsp = 1, nspec                                                     
     727             IF ( ibc_cs_b == 0 ) THEN                   
     728                DO l = 0, 1
    717729!--                Set index kb: For upward-facing surfaces (l=0), kb=-1, i.e.
    718730!--                the chem_species(nspec)%conc_p value at the topography top (k-1)
    … …  
    720732!--                value at the topography bottom (k+1) is set.
    721733
    722                       kb = MERGE( -1, 1, l == 0 )
    723                       !$OMP PARALLEL DO PRIVATE( i, j, k )
    724                       DO m = 1, bc_h(l)%ns
    725                          i = bc_h(l)%i(m)           
    726                          j = bc_h(l)%j(m)
    727                          k = bc_h(l)%k(m)
    728                          chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc(k+kb,j,i)
    729                       ENDDO                                       
    730                    ENDDO                                       
    731              
    732                 ELSEIF ( ibc_cs_b == 1 ) THEN
     734                   kb = MERGE( -1, 1, l == 0 )
     735                   !$OMP PARALLEL DO PRIVATE( i, j, k )
     736                   DO m = 1, bc_h(l)%ns
     737                      i = bc_h(l)%i(m)           
     738                      j = bc_h(l)%j(m)
     739                      k = bc_h(l)%k(m)
     740                      chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc(k+kb,j,i)
     741                   ENDDO                                       
     742                ENDDO                                       
     743
     744             ELSEIF ( ibc_cs_b == 1 ) THEN
    733745!--             in boundary_conds there is som extra loop over m here for passive tracer
    734                    DO l = 0, 1
    735                       kb = MERGE( -1, 1, l == 0 )
    736                       !$OMP PARALLEL DO PRIVATE( i, j, k )                                           
    737                       DO m = 1, bc_h(l)%ns
    738                          i = bc_h(l)%i(m)           
    739                          j = bc_h(l)%j(m)
    740                          k = bc_h(l)%k(m)
    741                          chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc_p(k,j,i)
    742 
    743                       ENDDO
     746                DO l = 0, 1
     747                   kb = MERGE( -1, 1, l == 0 )
     748                   !$OMP PARALLEL DO PRIVATE( i, j, k )                                           
     749                   DO m = 1, bc_h(l)%ns
     750                      i = bc_h(l)%i(m)           
     751                      j = bc_h(l)%j(m)
     752                      k = bc_h(l)%k(m)
     753                      chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc_p(k,j,i)
     754
    744755                   ENDDO
    745                 ENDIF
    746           ENDDO    ! end lsp loop 
    747 
    748 !--    Top boundary conditions for chemical species - Should this not be done for all species?
    749              IF ( ibc_cs_t == 0 )  THEN                   
    750                 DO lsp = 1, nspec
    751                    chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)       
    752                 ENDDO
    753              ELSEIF ( ibc_cs_t == 1 )  THEN
    754                 DO lsp = 1, nspec
    755                    chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
    756                 ENDDO
    757              ELSEIF ( ibc_cs_t == 2 )  THEN
    758                 DO lsp = 1, nspec
    759                    chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1)
    760756                ENDDO
    761757             ENDIF
    762 !
    763           CASE ( 'set_bc_lateral' )                       
     758       ENDDO    ! end lsp loop 
     759
     760!--    Top boundary conditions for chemical species - Should this not be done for all species?
     761          IF ( ibc_cs_t == 0 )  THEN                   
     762             DO lsp = 1, nspec
     763                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)       
     764             ENDDO
     765          ELSEIF ( ibc_cs_t == 1 )  THEN
     766             DO lsp = 1, nspec
     767                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
     768             ENDDO
     769          ELSEIF ( ibc_cs_t == 2 )  THEN
     770             DO lsp = 1, nspec
     771                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1)
     772             ENDDO
     773          ENDIF
     774!
     775       CASE ( 'set_bc_lateral' )                       
    764776!--             Lateral boundary conditions for chem species at inflow boundary
    765777!--             are automatically set when chem_species concentration is
    … …  
    770782!--             Lateral boundary conditions for chem species at outflow boundary
    771783
    772              IF ( bc_radiation_s )  THEN
    773                 DO lsp = 1, nspec
    774                    chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
    775                 ENDDO
    776              ELSEIF ( bc_radiation_n )  THEN
    777                 DO lsp = 1, nspec
    778                    chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
    779                 ENDDO
    780              ELSEIF ( bc_radiation_l )  THEN
    781                 DO lsp = 1, nspec
    782                    chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
    783                 ENDDO
    784              ELSEIF ( bc_radiation_r )  THEN
    785                 DO lsp = 1, nspec
    786                    chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
    787                 ENDDO
    788              ENDIF
    789            
    790        END SELECT
    791 
    792     END SUBROUTINE chem_boundary_conds
     784          IF ( bc_radiation_s )  THEN
     785             DO lsp = 1, nspec
     786                chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
     787             ENDDO
     788          ELSEIF ( bc_radiation_n )  THEN
     789             DO lsp = 1, nspec
     790                chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
     791             ENDDO
     792          ELSEIF ( bc_radiation_l )  THEN
     793             DO lsp = 1, nspec
     794                chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
     795             ENDDO
     796          ELSEIF ( bc_radiation_r )  THEN
     797             DO lsp = 1, nspec
     798                chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
     799             ENDDO
     800          ENDIF
     801
     802    END SELECT
     803
     804 END SUBROUTINE chem_boundary_conds
    793805
    794806!
    … …  
    799811!> x-direction
    800812!-----------------------------------------------------------------------------
    801     SUBROUTINE chem_boundary_conds_decycle (cs_3d, cs_pr_init )
    802        USE pegrid,                                                             &
    803                  ONLY: myid
    804 
    805        IMPLICIT NONE
    806        INTEGER(iwp) ::  boundary !<
    807        INTEGER(iwp) ::  ee !<
    808        INTEGER(iwp) ::  copied !<
    809        INTEGER(iwp) ::  i !<
    810        INTEGER(iwp) ::  j !<
    811        INTEGER(iwp) ::  k !<
    812        INTEGER(iwp) ::  ss !<
    813        REAL(wp), DIMENSION(nzb:nzt+1) ::  cs_pr_init
    814        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  cs_3d
    815        REAL(wp) ::  flag !< flag to mask topography grid points
    816 
    817        flag = 0.0_wp
    818 
    819        
    820 !--    Left and right boundaries
    821        IF ( decycle_chem_lr  .AND.  bc_lr_cyc )  THEN
    822 
    823           DO  boundary = 1, 2
    824 
    825              IF ( decycle_method(boundary) == 'dirichlet' )  THEN
     813 SUBROUTINE chem_boundary_conds_decycle( cs_3d, cs_pr_init )
     814
     815    USE pegrid,                                                             &
     816        ONLY:  myid
     817
     818    IMPLICIT NONE
     819    INTEGER(iwp) ::  boundary  !<
     820    INTEGER(iwp) ::  ee        !<
     821    INTEGER(iwp) ::  copied    !<
     822    INTEGER(iwp) ::  i         !<
     823    INTEGER(iwp) ::  j         !<
     824    INTEGER(iwp) ::  k         !<
     825    INTEGER(iwp) ::  ss        !<
     826    REAL(wp), DIMENSION(nzb:nzt+1) ::  cs_pr_init
     827    REAL(wp), DIMENSION(nzb:nzt+1, nysg:nyng, nxlg:nxrg) ::  cs_3d
     828    REAL(wp) ::  flag !< flag to mask topography grid points
     829
     830    flag = 0.0_wp
     831
     832
     833!-- Left and right boundaries
     834    IF ( decycle_chem_lr  .AND.  bc_lr_cyc )  THEN
     835
     836       DO  boundary = 1, 2
     837
     838          IF ( decycle_method(boundary) == 'dirichlet' )  THEN
    826839!   
    827 !--             Initial profile is copied to ghost and first three layers         
    828                 ss = 1
    829                 ee = 0
    830                 IF ( boundary == 1  .AND.  nxl == 0 )  THEN
    831                    ss = nxlg
    832                    ee = nxl+2
    833                 ELSEIF ( boundary == 2  .AND.  nxr == nx )  THEN
    834                    ss = nxr-2
    835                    ee = nxrg
    836                 ENDIF
    837 
    838                 DO  i = ss, ee
    839                    DO  j = nysg, nyng
    840                       DO  k = nzb+1, nzt
    841                          flag = MERGE( 1.0_wp, 0.0_wp,                            &
    842                                        BTEST( wall_flags_0(k,j,i), 0 ) )
    843                          cs_3d(k,j,i) = cs_pr_init(k) * flag
    844                       ENDDO
     840!--          Initial profile is copied to ghost and first three layers         
     841             ss = 1
     842             ee = 0
     843             IF ( boundary == 1  .AND.  nxl == 0 )  THEN
     844                ss = nxlg
     845                ee = nxl+2
     846             ELSEIF ( boundary == 2  .AND.  nxr == nx )  THEN
     847                ss = nxr-2
     848                ee = nxrg
     849             ENDIF
     850
     851             DO  i = ss, ee
     852                DO  j = nysg, nyng
     853                   DO  k = nzb+1, nzt
     854                      flag = MERGE( 1.0_wp, 0.0_wp,                            &
     855                                    BTEST( wall_flags_0(k,j,i), 0 ) )
     856                      cs_3d(k,j,i) = cs_pr_init(k) * flag
    845857                   ENDDO
    846858                ENDDO
    847 
    848            ELSEIF ( decycle_method(boundary) == 'neumann' )  THEN
    849 !
    850 !--             The value at the boundary is copied to the ghost layers to simulate
    851 !--             an outlet with zero gradient
    852                 ss = 1
    853                 ee = 0
    854                 IF ( boundary == 1  .AND.  nxl == 0 )  THEN
    855                    ss = nxlg
    856                    ee = nxl-1
    857                    copied = nxl
    858                 ELSEIF ( boundary == 2  .AND.  nxr == nx )  THEN
    859                    ss = nxr+1
    860                    ee = nxrg
    861                    copied = nxr
    862                 ENDIF
    863 
    864                  DO  i = ss, ee
    865                    DO  j = nysg, nyng
    866                       DO  k = nzb+1, nzt
    867                          flag = MERGE( 1.0_wp, 0.0_wp,                            &
    868                                        BTEST( wall_flags_0(k,j,i), 0 ) )
    869                         cs_3d(k,j,i) = cs_3d(k,j,copied) * flag
    870                       ENDDO
     859             ENDDO
     860
     861        ELSEIF ( decycle_method(boundary) == 'neumann' )  THEN
     862!
     863!--          The value at the boundary is copied to the ghost layers to simulate
     864!--          an outlet with zero gradient
     865             ss = 1
     866             ee = 0
     867             IF ( boundary == 1  .AND.  nxl == 0 )  THEN
     868                ss = nxlg
     869                ee = nxl-1
     870                copied = nxl
     871             ELSEIF ( boundary == 2  .AND.  nxr == nx )  THEN
     872                ss = nxr+1
     873                ee = nxrg
     874                copied = nxr
     875             ENDIF
     876
     877              DO  i = ss, ee
     878                DO  j = nysg, nyng
     879                   DO  k = nzb+1, nzt
     880                      flag = MERGE( 1.0_wp, 0.0_wp,                            &
     881                                    BTEST( wall_flags_0(k,j,i), 0 ) )
     882                     cs_3d(k,j,i) = cs_3d(k,j,copied) * flag
    871883                   ENDDO
    872884                ENDDO
    873 
    874              ELSE
    875                 WRITE(message_string,*)                                           &
    876                                     'unknown decycling method: decycle_method (', &
    877                         boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
    878                 CALL message( 'chem_boundary_conds_decycle', 'CM0431',           &
    879                               1, 2, 0, 6, 0 )
     885             ENDDO
     886
     887          ELSE
     888             WRITE(message_string,*)                                           &
     889                                 'unknown decycling method: decycle_method (', &
     890                     boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
     891             CALL message( 'chem_boundary_conds_decycle', 'CM0431',           &
     892                           1, 2, 0, 6, 0 )
     893          ENDIF
     894       ENDDO
     895    ENDIF
     896
     897
     898!-- South and north boundaries
     899    IF ( decycle_chem_ns  .AND.  bc_ns_cyc )  THEN
     900
     901       DO  boundary = 3, 4
     902
     903          IF ( decycle_method(boundary) == 'dirichlet' )  THEN
     904!   
     905!--          Initial profile is copied to ghost and first three layers         
     906             ss = 1
     907             ee = 0
     908             IF ( boundary == 3  .AND.  nys == 0 )  THEN
     909                ss = nysg
     910                ee = nys+2
     911             ELSEIF ( boundary == 4  .AND.  nyn == ny )  THEN
     912                ss = nyn-2
     913                ee = nyng
    880914             ENDIF
    881           ENDDO
    882        ENDIF
    883 
    884        
    885 !--    South and north boundaries
    886        IF ( decycle_chem_ns  .AND.  bc_ns_cyc )  THEN
    887 
    888           DO  boundary = 3, 4
    889 
    890              IF ( decycle_method(boundary) == 'dirichlet' )  THEN
    891 !   
    892 !--             Initial profile is copied to ghost and first three layers         
    893                 ss = 1
    894                 ee = 0
    895                 IF ( boundary == 3  .AND.  nys == 0 )  THEN
    896                    ss = nysg
    897                    ee = nys+2
    898                 ELSEIF ( boundary == 4  .AND.  nyn == ny )  THEN
    899                    ss = nyn-2
    900                    ee = nyng
    901                 ENDIF
    902 
    903                 DO  i = nxlg, nxrg
    904                    DO  j = ss, ee
    905                       DO  k = nzb+1, nzt
    906                          flag = MERGE( 1.0_wp, 0.0_wp,                            &
    907                                        BTEST( wall_flags_0(k,j,i), 0 ) )
    908                          cs_3d(k,j,i) = cs_pr_init(k) * flag
    909                       ENDDO
     915
     916             DO  i = nxlg, nxrg
     917                DO  j = ss, ee
     918                   DO  k = nzb+1, nzt
     919                      flag = MERGE( 1.0_wp, 0.0_wp,                            &
     920                                    BTEST( wall_flags_0(k,j,i), 0 ) )
     921                      cs_3d(k,j,i) = cs_pr_init(k) * flag
    910922                   ENDDO
    911923                ENDDO
    912        
    913        
    914            ELSEIF ( decycle_method(boundary) == 'neumann' )  THEN
    915 !
    916 !--             The value at the boundary is copied to the ghost layers to simulate
    917 !--             an outlet with zero gradient
    918                 ss = 1
    919                 ee = 0
    920                 IF ( boundary == 3  .AND.  nys == 0 )  THEN
    921                    ss = nysg
    922                    ee = nys-1
    923                    copied = nys
    924                 ELSEIF ( boundary == 4  .AND.  nyn == ny )  THEN
    925                    ss = nyn+1
    926                    ee = nyng
    927                    copied = nyn
    928                 ENDIF
    929 
    930                  DO  i = nxlg, nxrg
    931                    DO  j = ss, ee
    932                       DO  k = nzb+1, nzt
    933                          flag = MERGE( 1.0_wp, 0.0_wp,                            &
    934                                        BTEST( wall_flags_0(k,j,i), 0 ) )
    935                          cs_3d(k,j,i) = cs_3d(k,copied,i) * flag
    936                       ENDDO
     924             ENDDO
     925
     926
     927        ELSEIF ( decycle_method(boundary) == 'neumann' )  THEN
     928!
     929!--          The value at the boundary is copied to the ghost layers to simulate
     930!--          an outlet with zero gradient
     931             ss = 1
     932             ee = 0
     933             IF ( boundary == 3  .AND.  nys == 0 )  THEN
     934                ss = nysg
     935                ee = nys-1
     936                copied = nys
     937             ELSEIF ( boundary == 4  .AND.  nyn == ny )  THEN
     938                ss = nyn+1
     939                ee = nyng
     940                copied = nyn
     941             ENDIF
     942
     943              DO  i = nxlg, nxrg
     944                DO  j = ss, ee
     945                   DO  k = nzb+1, nzt
     946                      flag = MERGE( 1.0_wp, 0.0_wp,                            &
     947                                    BTEST( wall_flags_0(k,j,i), 0 ) )
     948                      cs_3d(k,j,i) = cs_3d(k,copied,i) * flag
    937949                   ENDDO
    938950                ENDDO
    939 
    940              ELSE
    941                 WRITE(message_string,*)                                           &
    942                                     'unknown decycling method: decycle_method (', &
    943                         boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
    944                 CALL message( 'chem_boundary_conds_decycle', 'CM0432',           &
    945                               1, 2, 0, 6, 0 )
    946              ENDIF
    947           ENDDO
    948        ENDIF
    949     END SUBROUTINE chem_boundary_conds_decycle
     951             ENDDO
     952
     953          ELSE
     954             WRITE(message_string,*)                                           &
     955                                 'unknown decycling method: decycle_method (', &
     956                     boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
     957             CALL message( 'chem_boundary_conds_decycle', 'CM0432',           &
     958                           1, 2, 0, 6, 0 )
     959          ENDIF
     960       ENDDO
     961    ENDIF
     962 END SUBROUTINE chem_boundary_conds_decycle
    950963   !
    951964!------------------------------------------------------------------------------!
    … …  
    956969!------------------------------------------------------------------------------!
    957970
    958     SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
    959 
    960 
    961        USE control_parameters,                                                 &
    962           ONLY: data_output, message_string
    963 
    964        IMPLICIT NONE
    965 
    966        CHARACTER (LEN=*) ::  unit     !<
    967        CHARACTER (LEN=*) ::  var      !<
    968 
    969        INTEGER(iwp) :: i, lsp
    970        INTEGER(iwp) :: ilen
    971        INTEGER(iwp) :: k
    972 
    973        CHARACTER(len=16)    ::  spec_name
    974 
    975        unit = 'illegal'
    976 
    977        spec_name = TRIM(var(4:))             !< var 1:3 is 'kc_' or 'em_'.
    978 
    979        IF ( TRIM(var(1:3)) == 'em_' )  THEN
    980           DO lsp=1,nspec
    981              IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
    982                 unit = 'mol m-2 s-1'
    983              ENDIF
    984              !-- It is possible to plant PM10 and PM25 into the gasphase chemistry code
    985              !-- as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
    986              !-- set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
    987              IF (spec_name(1:2) == 'PM')   THEN
    988                 unit = 'kg m-2 s-1'
    989              ENDIF
    990           ENDDO
    991 
    992        ELSE
    993 
    994           DO lsp=1,nspec
    995              IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
    996                 unit = 'ppm'
    997              ENDIF
    998 !            It is possible  to plant PM10 and PM25 into the gasphase chemistry code
    999 !            as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
    1000 !            set unit to kilograms per m**3 for PM10 and PM25 (PM2.5)
    1001              IF (spec_name(1:2) == 'PM')   THEN 
    1002                unit = 'kg m-3'
    1003              ENDIF
    1004           ENDDO
    1005 
    1006           DO lsp=1,nphot
    1007              IF (TRIM(spec_name) == TRIM(phot_frequen(lsp)%name))   THEN
    1008                 unit = 'sec-1'
    1009              ENDIF
    1010           ENDDO
    1011        ENDIF
    1012 
    1013 
    1014        RETURN
    1015     END SUBROUTINE chem_check_data_output
     971 SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
     972
     973
     974    USE control_parameters,                                                 &
     975        ONLY:  data_output, message_string
     976
     977    IMPLICIT NONE
     978
     979    CHARACTER (LEN=*) ::  unit     !<
     980    CHARACTER (LEN=*) ::  var      !<
     981
     982    INTEGER(iwp) ::  i
     983    INTEGER(iwp) ::  lsp
     984    INTEGER(iwp) ::  ilen
     985    INTEGER(iwp) ::  k
     986
     987    CHARACTER(len=16)    ::  spec_name
     988
     989    unit = 'illegal'
     990
     991    spec_name = TRIM(var(4:))             !< var 1:3 is 'kc_' or 'em_'.
     992
     993    IF ( TRIM(var(1:3)) == 'em_' )  THEN
     994       DO lsp=1,nspec
     995          IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
     996             unit = 'mol m-2 s-1'
     997          ENDIF
     998          !
     999          !-- It is possible to plant PM10 and PM25 into the gasphase chemistry code
     1000          !-- as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
     1001          !-- set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
     1002          IF (spec_name(1:2) == 'PM')   THEN
     1003             unit = 'kg m-2 s-1'
     1004          ENDIF
     1005       ENDDO
     1006
     1007    ELSE
     1008
     1009       DO lsp=1,nspec
     1010          IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
     1011             unit = 'ppm'
     1012          ENDIF
     1013!
     1014!--            It is possible  to plant PM10 and PM25 into the gasphase chemistry code
     1015!--            as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
     1016!--            set unit to kilograms per m**3 for PM10 and PM25 (PM2.5)
     1017          IF (spec_name(1:2) == 'PM')   THEN 
     1018            unit = 'kg m-3'
     1019          ENDIF
     1020       ENDDO
     1021
     1022       DO lsp=1,nphot
     1023          IF (TRIM(spec_name) == TRIM(phot_frequen(lsp)%name))   THEN
     1024             unit = 'sec-1'
     1025          ENDIF
     1026       ENDDO
     1027    ENDIF
     1028
     1029
     1030    RETURN
     1031 END SUBROUTINE chem_check_data_output
    10161032!
    10171033!------------------------------------------------------------------------------!
    … …  
    10221038!------------------------------------------------------------------------------!
    10231039
    1024     SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
    1025 
    1026        USE arrays_3d
    1027        USE control_parameters,                                                    &
    1028            ONLY: data_output_pr, message_string, air_chemistry
    1029        USE indices
    1030        USE profil_parameter
    1031        USE statistics
    1032 
    1033 
    1034        IMPLICIT NONE
    1035 
    1036        CHARACTER (LEN=*) ::  unit     !<
    1037        CHARACTER (LEN=*) ::  variable !<
    1038        CHARACTER (LEN=*) ::  dopr_unit
    1039        CHARACTER(len=16) ::  spec_name
    1040    
    1041        INTEGER(iwp) ::  var_count, lsp  !<
    1042        
    1043 
    1044        spec_name = TRIM(variable(4:))   
    1045 
    1046           IF (  .NOT.  air_chemistry )  THEN
    1047              message_string = 'data_output_pr = ' //                        &
    1048              TRIM( data_output_pr(var_count) ) // ' is not imp' // &
    1049                          'lemented for air_chemistry = .FALSE.'
    1050              CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 )             
    1051 
    1052           ELSE
    1053              DO lsp = 1, nspec
    1054                 IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN
    1055                    cs_pr_count = cs_pr_count+1
    1056                    cs_pr_index(cs_pr_count) = lsp
    1057                    dopr_index(var_count) = pr_palm+cs_pr_count 
    1058                    dopr_unit  = 'ppm'
    1059                    IF (spec_name(1:2) == 'PM')   THEN
    1060                         dopr_unit = 'kg m-3'
    1061                    ENDIF
    1062                       hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
    1063                       unit = dopr_unit 
     1040 SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
     1041
     1042    USE arrays_3d
     1043    USE control_parameters,                                                    &
     1044        ONLY:  air_chemistry, data_output_pr, message_string
     1045    USE indices
     1046    USE profil_parameter
     1047    USE statistics
     1048
     1049
     1050    IMPLICIT NONE
     1051
     1052    CHARACTER (LEN=*) ::  unit     !<
     1053    CHARACTER (LEN=*) ::  variable !<
     1054    CHARACTER (LEN=*) ::  dopr_unit
     1055    CHARACTER(len=16) ::  spec_name
     1056
     1057    INTEGER(iwp) ::  var_count, lsp  !<
     1058
     1059
     1060    spec_name = TRIM(variable(4:))   
     1061
     1062       IF (  .NOT.  air_chemistry )  THEN
     1063          message_string = 'data_output_pr = ' //                        &
     1064          TRIM( data_output_pr(var_count) ) // ' is not imp' // &
     1065                      'lemented for air_chemistry = .FALSE.'
     1066          CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 )             
     1067
     1068       ELSE
     1069          DO lsp = 1, nspec
     1070             IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN
     1071                cs_pr_count = cs_pr_count+1
     1072                cs_pr_index(cs_pr_count) = lsp
     1073                dopr_index(var_count) = pr_palm+cs_pr_count 
     1074                dopr_unit  = 'ppm'
     1075                IF (spec_name(1:2) == 'PM')   THEN
     1076                     dopr_unit = 'kg m-3'
    10641077                ENDIF
    1065              ENDDO
    1066           ENDIF
    1067 
    1068     END SUBROUTINE chem_check_data_output_pr
     1078                   hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
     1079                   unit = dopr_unit 
     1080             ENDIF
     1081          ENDDO
     1082       ENDIF
     1083
     1084 END SUBROUTINE chem_check_data_output_pr
    10691085
    10701086!
    … …  
    10741090!> Check parameters routine for chemistry_model_mod
    10751091!------------------------------------------------------------------------------!
    1076     SUBROUTINE chem_check_parameters
    1077 
    1078        IMPLICIT NONE
    1079 
    1080        LOGICAL                        ::  found
    1081        INTEGER (iwp)                  ::  lsp_usr      !< running index for user defined chem spcs
    1082        INTEGER (iwp)                  ::  lsp          !< running index for chem spcs.
    1083 
    1084 
    1085 !!--   check for chemical reactions status
    1086        IF ( chem_gasphase_on )  THEN
    1087           message_string = 'Chemical reactions: ON'
    1088           CALL message( 'chem_check_parameters', 'CM0421', 0, 0, 0, 6, 0 )
    1089        ELSEIF ( .not. (chem_gasphase_on) ) THEN
    1090           message_string = 'Chemical reactions: OFF'
    1091           CALL message( 'chem_check_parameters', 'CM0422', 0, 0, 0, 6, 0 )
     1092 SUBROUTINE chem_check_parameters
     1093
     1094    IMPLICIT NONE
     1095
     1096    LOGICAL       ::  found
     1097    INTEGER (iwp) ::  lsp_usr      !< running index for user defined chem spcs
     1098    INTEGER (iwp) ::  lsp          !< running index for chem spcs.
     1099
     1100
     1101!
     1102!-- check for chemical reactions status
     1103    IF ( chem_gasphase_on )  THEN
     1104       message_string = 'Chemical reactions: ON'
     1105       CALL message( 'chem_check_parameters', 'CM0421', 0, 0, 0, 6, 0 )
     1106    ELSEIF ( .not. (chem_gasphase_on) ) THEN
     1107       message_string = 'Chemical reactions: OFF'
     1108       CALL message( 'chem_check_parameters', 'CM0422', 0, 0, 0, 6, 0 )
     1109    ENDIF
     1110
     1111!-- check for chemistry time-step
     1112    IF ( call_chem_at_all_substeps )  THEN
     1113       message_string = 'Chemistry is calculated at all meteorology time-step'
     1114       CALL message( 'chem_check_parameters', 'CM0423', 0, 0, 0, 6, 0 )
     1115    ELSEIF ( .not. (call_chem_at_all_substeps) ) THEN
     1116       message_string = 'Sub-time-steps are skipped for chemistry time-steps'
     1117       CALL message( 'chem_check_parameters', 'CM0424', 0, 0, 0, 6, 0 )
     1118    ENDIF
     1119
     1120!-- check for photolysis scheme
     1121    IF ( (photolysis_scheme /= 'simple') .AND. (photolysis_scheme /= 'constant')  )  THEN
     1122       message_string = 'Incorrect photolysis scheme selected, please check spelling'
     1123       CALL message( 'chem_check_parameters', 'CM0425', 1, 2, 0, 6, 0 )
     1124    ENDIF
     1125
     1126!-- check for  decycling of chem species
     1127    IF ( (.not. any(decycle_method == 'neumann') ) .AND. (.not. any(decycle_method == 'dirichlet') ) )   THEN
     1128       message_string = 'Incorrect boundary conditions. Only neumann or '   &
     1129                // 'dirichlet &available for decycling chemical species '
     1130       CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 )
     1131    ENDIF
     1132
     1133!---------------------
     1134!-- chem_check_parameters is called before the array chem_species is allocated!
     1135!-- temporary switch of this part of the check
     1136    RETURN
     1137!---------------------
     1138
     1139!-- check for initial chem species input
     1140    lsp_usr = 1
     1141    lsp     = 1
     1142    DO WHILE ( cs_name (lsp_usr) /= 'novalue')
     1143       found = .FALSE.
     1144       DO lsp = 1, nvar
     1145          IF ( TRIM(cs_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN
     1146             found = .TRUE.
     1147             EXIT
     1148          ENDIF
     1149       ENDDO
     1150       IF ( .not. found ) THEN
     1151          message_string = 'Unused/incorrect input for initial surface value: ' // TRIM(cs_name(lsp_usr))
     1152          CALL message( 'chem_check_parameters', 'CM0427', 0, 1, 0, 6, 0 )
    10921153       ENDIF
    1093 
    1094 !--    check for chemistry time-step
    1095        IF ( call_chem_at_all_substeps )  THEN
    1096           message_string = 'Chemistry is calculated at all meteorology time-step'
    1097           CALL message( 'chem_check_parameters', 'CM0423', 0, 0, 0, 6, 0 )
    1098        ELSEIF ( .not. (call_chem_at_all_substeps) ) THEN
    1099           message_string = 'Sub-time-steps are skipped for chemistry time-steps'
    1100           CALL message( 'chem_check_parameters', 'CM0424', 0, 0, 0, 6, 0 )
     1154       lsp_usr = lsp_usr + 1
     1155    ENDDO
     1156
     1157!-- check for surface  emission flux chem species
     1158
     1159    lsp_usr = 1
     1160    lsp     = 1
     1161    DO WHILE ( surface_csflux_name (lsp_usr) /= 'novalue')
     1162       found = .FALSE.
     1163       DO lsp = 1, nvar
     1164          IF ( TRIM(surface_csflux_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN
     1165             found = .TRUE.
     1166             EXIT
     1167          ENDIF
     1168       ENDDO
     1169       IF ( .not. found ) THEN
     1170          message_string = 'Unused/incorrect input of chemical species for surface emission fluxes: '  &
     1171                            // TRIM(surface_csflux_name(lsp_usr))
     1172          CALL message( 'chem_check_parameters', 'CM0428', 0, 1, 0, 6, 0 )
    11011173       ENDIF
    1102 
    1103 !--    check for photolysis scheme
    1104        IF ( (photolysis_scheme /= 'simple') .AND. (photolysis_scheme /= 'constant')  )  THEN
    1105           message_string = 'Incorrect photolysis scheme selected, please check spelling'
    1106           CALL message( 'chem_check_parameters', 'CM0425', 1, 2, 0, 6, 0 )
    1107        ENDIF
    1108 
    1109 !--    check for  decycling of chem species
    1110        IF ( (.not. any(decycle_method == 'neumann') ) .AND. (.not. any(decycle_method == 'dirichlet') ) )   THEN
    1111           message_string = 'Incorrect boundary conditions. Only neumann or '   &
    1112                    // 'dirichlet &available for decycling chemical species '
    1113           CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 )
    1114        ENDIF
    1115 
    1116 !---------------------
    1117 !--    chem_check_parameters is called before the array chem_species is allocated!
    1118 !--    temporary switch of this part of the check
    1119        RETURN
    1120 !---------------------
    1121 
    1122 !--    check for initial chem species input
    1123        lsp_usr = 1
    1124        lsp     = 1
    1125        DO WHILE ( cs_name (lsp_usr) /= 'novalue')
    1126           found = .FALSE.
    1127           DO lsp = 1, nvar
    1128              IF ( TRIM(cs_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN
    1129                 found = .TRUE.
    1130                 EXIT
    1131              ENDIF
    1132           ENDDO
    1133           IF ( .not. found ) THEN
    1134              message_string = 'Incorrect input for initial surface vaue: ' // TRIM(cs_name(lsp_usr))
    1135              CALL message( 'chem_check_parameters', 'CM0427', 0, 1, 0, 6, 0 )
    1136           ENDIF
    1137           lsp_usr = lsp_usr + 1
    1138        ENDDO
    1139 
    1140 !--    check for surface  emission flux chem species
    1141 
    1142        lsp_usr = 1
    1143        lsp     = 1
    1144        DO WHILE ( surface_csflux_name (lsp_usr) /= 'novalue')
    1145           found = .FALSE.
    1146           DO lsp = 1, nvar
    1147              IF ( TRIM(surface_csflux_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN
    1148                 found = .TRUE.
    1149                 EXIT
    1150              ENDIF
    1151           ENDDO
    1152           IF ( .not. found ) THEN
    1153              message_string = 'Incorrect input of chemical species for surface emission fluxes: '  &
    1154                                // TRIM(surface_csflux_name(lsp_usr))
    1155              CALL message( 'chem_check_parameters', 'CM0428', 0, 1, 0, 6, 0 )
    1156           ENDIF
    1157           lsp_usr = lsp_usr + 1
    1158        ENDDO
    1159 
    1160     END SUBROUTINE chem_check_parameters
     1174       lsp_usr = lsp_usr + 1
     1175    ENDDO
     1176
     1177 END SUBROUTINE chem_check_parameters
    11611178
    11621179!
    … …  
    11691186!------------------------------------------------------------------------------!
    11701187   
    1171     SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf,   &
    1172                                      two_d, nzb_do, nzt_do, fill_value )
    1173    
    1174        USE indices
    1175 
    1176        USE kinds
    1177 
    1178        USE pegrid,             ONLY: myid, threads_per_task
    1179 
    1180        IMPLICIT NONE
    1181 
    1182        CHARACTER (LEN=*) ::  grid       !<
    1183        CHARACTER (LEN=*) ::  mode       !<
    1184        CHARACTER (LEN=*) ::  variable   !<
    1185        INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or time-averaged data
    1186        INTEGER(iwp) ::  nzb_do          !< lower limit of the domain (usually nzb)
    1187        INTEGER(iwp) ::  nzt_do          !< upper limit of the domain (usually nzt+1)
    1188        LOGICAL      ::  found           !<
    1189        LOGICAL      ::  two_d           !< flag parameter that indicates 2D variables (horizontal cross sections)
    1190        REAL(wp)     ::  fill_value
    1191        REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf !<
    1192 
    1193        !-- local variables.
    1194        CHARACTER(len=16)    ::  spec_name
    1195        INTEGER(iwp) ::  lsp
    1196        INTEGER(iwp) ::  i               !< grid index along x-direction
    1197        INTEGER(iwp) ::  j               !< grid index along y-direction
    1198        INTEGER(iwp) ::  k               !< grid index along z-direction
    1199        INTEGER(iwp) ::  m               !< running index surface elements
    1200        INTEGER(iwp) ::  char_len        !< length of a character string
    1201        found = .TRUE.
    1202        char_len  = LEN_TRIM(variable)
    1203 
    1204        spec_name = TRIM( variable(4:char_len-3) )
    1205 
    1206           DO lsp=1,nspec
    1207              IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)    .AND.                             &
    1208                    ( (variable(char_len-2:) == '_xy')               .OR.                              &
    1209                      (variable(char_len-2:) == '_xz')               .OR.                              &
    1210                      (variable(char_len-2:) == '_yz') ) )               THEN             
     1188 SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf,   &
     1189                                  two_d, nzb_do, nzt_do, fill_value )
     1190
     1191    USE indices
     1192
     1193    USE kinds
     1194
     1195    USE pegrid,                                                               &
     1196        ONLY:  myid, threads_per_task
     1197
     1198    IMPLICIT NONE
     1199
     1200    CHARACTER (LEN=*) ::  grid       !<
     1201    CHARACTER (LEN=*) ::  mode       !<
     1202    CHARACTER (LEN=*) ::  variable   !<
     1203    INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or time-averaged data
     1204    INTEGER(iwp) ::  nzb_do          !< lower limit of the domain (usually nzb)
     1205    INTEGER(iwp) ::  nzt_do          !< upper limit of the domain (usually nzt+1)
     1206    LOGICAL      ::  found           !<
     1207    LOGICAL      ::  two_d           !< flag parameter that indicates 2D variables (horizontal cross sections)
     1208    REAL(wp)     ::  fill_value
     1209    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf !<
     1210
     1211    !
     1212    !-- local variables.
     1213    CHARACTER(len=16)    ::  spec_name
     1214    INTEGER(iwp) ::  lsp
     1215    INTEGER(iwp) ::  i               !< grid index along x-direction
     1216    INTEGER(iwp) ::  j               !< grid index along y-direction
     1217    INTEGER(iwp) ::  k               !< grid index along z-direction
     1218    INTEGER(iwp) ::  m               !< running index surface elements
     1219    INTEGER(iwp) ::  char_len        !< length of a character string
     1220    found = .TRUE.
     1221    char_len  = LEN_TRIM(variable)
     1222
     1223    spec_name = TRIM( variable(4:char_len-3) )
     1224
     1225       DO lsp=1,nspec
     1226          IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)    .AND.                             &
     1227                ( (variable(char_len-2:) == '_xy')               .OR.                              &
     1228                  (variable(char_len-2:) == '_xz')               .OR.                              &
     1229                  (variable(char_len-2:) == '_yz') ) )               THEN             
    12111230!
    12121231!--   todo: remove or replace by "CALL message" mechanism (kanani)
    12131232!                    IF(myid == 0)  WRITE(6,*) 'Output of species ' // TRIM(variable)  //               &
    12141233!                                                             TRIM(chem_species(lsp)%name)       
    1215                 IF (av == 0) THEN
    1216                    DO  i = nxl, nxr
    1217                       DO  j = nys, nyn
    1218                          DO  k = nzb_do, nzt_do
    1219                               local_pf(i,j,k) = MERGE(                                                &
    1220                                                  chem_species(lsp)%conc(k,j,i),                       &
    1221                                                  REAL( fill_value, KIND = wp ),                       &
    1222                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
    1223 
    1224 
    1225                          ENDDO
     1234             IF (av == 0) THEN
     1235                DO  i = nxl, nxr
     1236                   DO  j = nys, nyn
     1237                      DO  k = nzb_do, nzt_do
     1238                           local_pf(i,j,k) = MERGE(                                                &
     1239                                              chem_species(lsp)%conc(k,j,i),                       &
     1240                                              REAL( fill_value, KIND = wp ),                       &
     1241                                              BTEST( wall_flags_0(k,j,i), 0 ) )
     1242
     1243
    12261244                      ENDDO
    12271245                   ENDDO
    1228              
    1229                 ELSE
    1230                    DO  i = nxl, nxr
    1231                       DO  j = nys, nyn
    1232                          DO  k = nzb_do, nzt_do
    1233                               local_pf(i,j,k) = MERGE(                                                &
    1234                                                  chem_species(lsp)%conc_av(k,j,i),                    &
    1235                                                  REAL( fill_value, KIND = wp ),                       &
    1236                                                  BTEST( wall_flags_0(k,j,i), 0 ) )
    1237                          ENDDO
     1246                ENDDO
     1247
     1248             ELSE
     1249                DO  i = nxl, nxr
     1250                   DO  j = nys, nyn
     1251                      DO  k = nzb_do, nzt_do
     1252                           local_pf(i,j,k) = MERGE(                                                &
     1253                                              chem_species(lsp)%conc_av(k,j,i),                    &
     1254                                              REAL( fill_value, KIND = wp ),                       &
     1255                                              BTEST( wall_flags_0(k,j,i), 0 ) )
    12381256                      ENDDO
    12391257                   ENDDO
    1240                 ENDIF
    1241                  grid = 'zu'           
     1258                ENDDO
    12421259             ENDIF
    1243           ENDDO
    1244 
    1245           RETURN
    1246      
    1247     END SUBROUTINE chem_data_output_2d     
     1260              grid = 'zu'           
     1261          ENDIF
     1262       ENDDO
     1263
     1264       RETURN
     1265
     1266 END SUBROUTINE chem_data_output_2d     
    12481267
    12491268!
    … …  
    12551274!------------------------------------------------------------------------------!
    12561275
    1257     SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
    1258 
    1259 
    1260        USE indices
    1261 
    1262        USE kinds
    1263 
    1264        USE surface_mod
    1265 
    1266 
    1267        IMPLICIT NONE
    1268 
    1269        CHARACTER (LEN=*)    ::  variable     !<
    1270        INTEGER(iwp)         ::  av           !<
    1271        INTEGER(iwp) ::  nzb_do               !< lower limit of the data output (usually 0)
    1272        INTEGER(iwp) ::  nzt_do               !< vertical upper limit of the data output (usually nz_do3d)
    1273 
    1274        LOGICAL      ::  found                !<
    1275 
    1276        REAL(wp)             ::  fill_value !<
    1277        REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf
    1278 
    1279 
    1280        !-- local variables
    1281        CHARACTER(len=16)    ::  spec_name
    1282        INTEGER(iwp)         ::  i, j, k
    1283        INTEGER(iwp)         ::  m, l    !< running indices for surfaces
    1284        INTEGER(iwp)         ::  lsp  !< running index for chem spcs
    1285 
    1286 
    1287        found = .FALSE.
    1288        IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) ) THEN
    1289           RETURN
    1290        ENDIF
    1291 
    1292        spec_name = TRIM(variable(4:))
    1293 
    1294        IF ( variable(1:3) == 'em_' ) THEN
    1295 
    1296          local_pf = 0.0_wp
    1297 
    1298          DO lsp = 1, nvar  !!! cssws - nvar species, chem_species - nspec species !!!
    1299             IF ( TRIM(spec_name) == TRIM(chem_species(lsp)%name) )   THEN
    1300                ! no average for now
    1301                DO m = 1, surf_usm_h%ns
    1302                   local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) = &
    1303                      local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) + surf_usm_h%cssws(lsp,m)
     1276 SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
     1277
     1278
     1279    USE indices
     1280
     1281    USE kinds
     1282
     1283    USE surface_mod
     1284
     1285
     1286    IMPLICIT NONE
     1287
     1288    CHARACTER (LEN=*)    ::  variable     !<
     1289    INTEGER(iwp)         ::  av           !<
     1290    INTEGER(iwp) ::  nzb_do               !< lower limit of the data output (usually 0)
     1291    INTEGER(iwp) ::  nzt_do               !< vertical upper limit of the data output (usually nz_do3d)
     1292
     1293    LOGICAL      ::  found                !<
     1294
     1295    REAL(wp)             ::  fill_value   !<
     1296    REAL(sp), DIMENSION(nxl:nxr, nys:nyn, nzb_do:nzt_do) ::  local_pf
     1297
     1298
     1299    !-- local variables
     1300    CHARACTER(len=16)    ::  spec_name
     1301    INTEGER(iwp)         ::  i
     1302    INTEGER(iwp)         ::  j
     1303    INTEGER(iwp)         ::  k
     1304    INTEGER(iwp)         ::  m       !< running indices for surfaces
     1305    INTEGER(iwp)         ::  l
     1306    INTEGER(iwp)         ::  lsp     !< running index for chem spcs
     1307
     1308
     1309    found = .FALSE.
     1310    IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) ) THEN
     1311       RETURN
     1312    ENDIF
     1313
     1314    spec_name = TRIM(variable(4:))
     1315
     1316    IF ( variable(1:3) == 'em_' ) THEN
     1317
     1318      local_pf = 0.0_wp
     1319
     1320      DO lsp = 1, nvar   !!! cssws - nvar species, chem_species - nspec species !!!
     1321         IF ( TRIM(spec_name) == TRIM(chem_species(lsp)%name) )   THEN
     1322            !
     1323            !-- no average for now
     1324            DO m = 1, surf_usm_h%ns
     1325               local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) = &
     1326                  local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) + surf_usm_h%cssws(lsp,m)
     1327            ENDDO
     1328            DO m = 1, surf_lsm_h%ns
     1329               local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) = &
     1330                  local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) + surf_lsm_h%cssws(lsp,m)
     1331            ENDDO
     1332            DO l = 0, 3
     1333               DO m = 1, surf_usm_v(l)%ns
     1334                  local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) = &
     1335                     local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) + surf_usm_v(l)%cssws(lsp,m)
    13041336               ENDDO
    1305                DO m = 1, surf_lsm_h%ns
    1306                   local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) = &
    1307                      local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) + surf_lsm_h%cssws(lsp,m)
     1337               DO m = 1, surf_lsm_v(l)%ns
     1338                  local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) = &
     1339                     local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) + surf_lsm_v(l)%cssws(lsp,m)
    13081340               ENDDO
    1309                DO l = 0, 3
    1310                   DO m = 1, surf_usm_v(l)%ns
    1311                      local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) = &
    1312                         local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) + surf_usm_v(l)%cssws(lsp,m)
    1313                   ENDDO
    1314                   DO m = 1, surf_lsm_v(l)%ns
    1315                      local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) = &
    1316                         local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) + surf_lsm_v(l)%cssws(lsp,m)
    1317                   ENDDO
    1318                ENDDO
    1319                found = .TRUE.
    1320             ENDIF
    1321          ENDDO
    1322        ELSE
    1323          DO lsp=1,nspec
    1324             IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
     1341            ENDDO
     1342            found = .TRUE.
     1343         ENDIF
     1344      ENDDO
     1345    ELSE
     1346      DO lsp=1,nspec
     1347         IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
    13251348!
    13261349!--   todo: remove or replace by "CALL message" mechanism (kanani)
    13271350!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM(variable)  //  &
    13281351!                                                           TRIM(chem_species(lsp)%name)       
    1329                IF (av == 0) THEN
    1330                   DO  i = nxl, nxr
    1331                      DO  j = nys, nyn
    1332                         DO  k = nzb_do, nzt_do
    1333                             local_pf(i,j,k) = MERGE(                             &
    1334                                                 chem_species(lsp)%conc(k,j,i),   &
    1335                                                 REAL( fill_value, KIND = wp ),   &
    1336                                                 BTEST( wall_flags_0(k,j,i), 0 ) )
    1337                         ENDDO
     1352            IF (av == 0) THEN
     1353               DO  i = nxl, nxr
     1354                  DO  j = nys, nyn
     1355                     DO  k = nzb_do, nzt_do
     1356                         local_pf(i,j,k) = MERGE(                             &
     1357                                             chem_species(lsp)%conc(k,j,i),   &
     1358                                             REAL( fill_value, KIND = wp ),   &
     1359                                             BTEST( wall_flags_0(k,j,i), 0 ) )
    13381360                     ENDDO
    13391361                  ENDDO
    1340 
    1341                ELSE
    1342                   DO  i = nxl, nxr
    1343                      DO  j = nys, nyn
    1344                         DO  k = nzb_do, nzt_do
    1345                             local_pf(i,j,k) = MERGE(                             &
    1346                                                 chem_species(lsp)%conc_av(k,j,i),&
    1347                                                 REAL( fill_value, KIND = wp ),   &
    1348                                                 BTEST( wall_flags_0(k,j,i), 0 ) )
    1349                         ENDDO
     1362               ENDDO
     1363
     1364            ELSE
     1365               DO  i = nxl, nxr
     1366                  DO  j = nys, nyn
     1367                     DO  k = nzb_do, nzt_do
     1368                         local_pf(i,j,k) = MERGE(                             &
     1369                                             chem_species(lsp)%conc_av(k,j,i),&
     1370                                             REAL( fill_value, KIND = wp ),   &
     1371                                             BTEST( wall_flags_0(k,j,i), 0 ) )
    13501372                     ENDDO
    13511373                  ENDDO
    1352                ENDIF
    1353                found = .TRUE.
     1374               ENDDO
    13541375            ENDIF
    1355          ENDDO
    1356        ENDIF
    1357 
    1358        RETURN
    1359 
    1360     END SUBROUTINE chem_data_output_3d
     1376            found = .TRUE.
     1377         ENDIF
     1378      ENDDO
     1379    ENDIF
     1380
     1381    RETURN
     1382
     1383 END SUBROUTINE chem_data_output_3d
    13611384!
    13621385!------------------------------------------------------------------------------!
    … …  
    13671390!------------------------------------------------------------------------------!
    13681391   
    1369     SUBROUTINE chem_data_output_mask( av, variable, found, local_pf )
    1370    
    1371        USE control_parameters
    1372        USE indices
    1373        USE kinds
    1374        USE pegrid,             ONLY: myid, threads_per_task
    1375        USE surface_mod,        ONLY: get_topography_top_index_ji
    1376 
    1377        IMPLICIT NONE
    1378 
    1379        CHARACTER(LEN=5) ::  grid        !< flag to distinquish between staggered grids
    1380 
    1381        CHARACTER (LEN=*)::  variable    !<
    1382        INTEGER(iwp)     ::  av          !< flag to control data output of instantaneous or time-averaged data
    1383        LOGICAL          ::  found       !<
    1384        REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
    1385                  local_pf   !<
    1386 
    1387 
    1388        !-- local variables.
    1389        CHARACTER(len=16)    ::  spec_name
    1390        INTEGER(iwp) ::  lsp
    1391        INTEGER(iwp) ::  i               !< grid index along x-direction
    1392        INTEGER(iwp) ::  j               !< grid index along y-direction
    1393        INTEGER(iwp) ::  k               !< grid index along z-direction
    1394        INTEGER(iwp) ::  topo_top_ind    !< k index of highest horizontal surface
    1395 
    1396        found = .TRUE.
    1397        grid  = 's'
    1398 
    1399        spec_name = TRIM( variable(4:) )
    1400 
    1401        DO lsp=1,nspec
    1402           IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name) )               THEN             
    1403 !
    1404 !--   todo: remove or replace by "CALL message" mechanism (kanani)
    1405 !                 IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM(variable)  // &
    1406 !                                                           TRIM(chem_species(lsp)%name)       
    1407              IF (av == 0) THEN
    1408                 IF ( .NOT. mask_surface(mid) )  THEN
    1409 
    1410                    DO  i = 1, mask_size_l(mid,1)
    1411                       DO  j = 1, mask_size_l(mid,2)
    1412                          DO  k = 1, mask_size(mid,3)
    1413                              local_pf(i,j,k) = chem_species(lsp)%conc(  &
    1414                                                   mask_k(mid,k),        &
    1415                                                   mask_j(mid,j),        &
    1416                                                   mask_i(mid,i)      )
    1417                          ENDDO
     1392 SUBROUTINE chem_data_output_mask( av, variable, found, local_pf )
     1393
     1394    USE control_parameters
     1395    USE indices
     1396    USE kinds
     1397    USE pegrid,                                                                       &
     1398        ONLY:  myid, threads_per_task
     1399    USE surface_mod,                                                                  &
     1400        ONLY:  get_topography_top_index_ji
     1401
     1402    IMPLICIT NONE
     1403
     1404    CHARACTER(LEN=5) ::  grid        !< flag to distinquish between staggered grids
     1405
     1406    CHARACTER (LEN=*)::  variable    !<
     1407    INTEGER(iwp)     ::  av          !< flag to control data output of instantaneous or time-averaged data
     1408    LOGICAL          ::  found       !<
     1409    REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
     1410              local_pf   !<
     1411
     1412
     1413    !-- local variables.
     1414    CHARACTER(len=16)    ::  spec_name
     1415    INTEGER(iwp) ::  lsp
     1416    INTEGER(iwp) ::  i               !< grid index along x-direction
     1417    INTEGER(iwp) ::  j               !< grid index along y-direction
     1418    INTEGER(iwp) ::  k               !< grid index along z-direction
     1419    INTEGER(iwp) ::  topo_top_ind    !< k index of highest horizontal surface
     1420
     1421    found = .TRUE.
     1422    grid  = 's'
     1423
     1424    spec_name = TRIM( variable(4:) )
     1425
     1426    DO lsp=1,nspec
     1427       IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name) )               THEN             
     1428!
     1429!-- todo: remove or replace by "CALL message" mechanism (kanani)
     1430!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM(variable)  // &
     1431!                                                        TRIM(chem_species(lsp)%name)       
     1432          IF (av == 0) THEN
     1433             IF ( .NOT. mask_surface(mid) )  THEN
     1434
     1435                DO  i = 1, mask_size_l(mid,1)
     1436                   DO  j = 1, mask_size_l(mid,2)
     1437                      DO  k = 1, mask_size(mid,3)
     1438                          local_pf(i,j,k) = chem_species(lsp)%conc(  &
     1439                                               mask_k(mid,k),        &
     1440                                               mask_j(mid,j),        &
     1441                                               mask_i(mid,i)      )
    14181442                      ENDDO
    14191443                   ENDDO
    1420 
    1421                 ELSE
     1444                ENDDO
     1445
     1446             ELSE
    14221447!             
    1423 !--                Terrain-following masked output
    1424                    DO  i = 1, mask_size_l(mid,1)
    1425                       DO  j = 1, mask_size_l(mid,2)
     1448!--             Terrain-following masked output
     1449                DO  i = 1, mask_size_l(mid,1)
     1450                   DO  j = 1, mask_size_l(mid,2)
    14261451!             
    1427 !--                      Get k index of highest horizontal surface
    1428                          topo_top_ind = get_topography_top_index_ji( &
    1429                                            mask_j(mid,j),  &
    1430                                            mask_i(mid,i),  &
    1431                                            grid                    )
     1452!--                   Get k index of highest horizontal surface
     1453                      topo_top_ind = get_topography_top_index_ji( &
     1454                                        mask_j(mid,j),  &
     1455                                        mask_i(mid,i),  &
     1456                                        grid                    )
    14321457!             
    1433 !--                      Save output array
    1434                          DO  k = 1, mask_size_l(mid,3)
    1435                             local_pf(i,j,k) = chem_species(lsp)%conc( &
    1436                                                  MIN( topo_top_ind+mask_k(mid,k), &
    1437                                                       nzt+1 ),        &
    1438                                                  mask_j(mid,j),       &
    1439                                                  mask_i(mid,i)      )
    1440                          ENDDO
     1458!--                   Save output array
     1459                      DO  k = 1, mask_size_l(mid,3)
     1460                         local_pf(i,j,k) = chem_species(lsp)%conc( &
     1461                                              MIN( topo_top_ind+mask_k(mid,k), &
     1462                                                   nzt+1 ),        &
     1463                                              mask_j(mid,j),       &
     1464                                              mask_i(mid,i)      )
    14411465                      ENDDO
    14421466                   ENDDO
    1443 
    1444                 ENDIF
    1445              ELSE
    1446                 IF ( .NOT. mask_surface(mid) )  THEN
    1447 
    1448                    DO  i = 1, mask_size_l(mid,1)
    1449                       DO  j = 1, mask_size_l(mid,2)
    1450                          DO  k =  1, mask_size_l(mid,3)
    1451                              local_pf(i,j,k) = chem_species(lsp)%conc_av(  &
    1452                                                   mask_k(mid,k),           &
    1453                                                   mask_j(mid,j),           &
    1454                                                   mask_i(mid,i)         )
    1455                          ENDDO
     1467                ENDDO
     1468
     1469             ENDIF
     1470          ELSE
     1471             IF ( .NOT. mask_surface(mid) )  THEN
     1472
     1473                DO  i = 1, mask_size_l(mid,1)
     1474                   DO  j = 1, mask_size_l(mid,2)
     1475                      DO  k =  1, mask_size_l(mid,3)
     1476                          local_pf(i,j,k) = chem_species(lsp)%conc_av(  &
     1477                                               mask_k(mid,k),           &
     1478                                               mask_j(mid,j),           &
     1479                                               mask_i(mid,i)         )
    14561480                      ENDDO
    14571481                   ENDDO
    1458 
    1459                 ELSE
     1482                ENDDO
     1483
     1484             ELSE
    14601485!             
    1461 !--                Terrain-following masked output
    1462                    DO  i = 1, mask_size_l(mid,1)
    1463                       DO  j = 1, mask_size_l(mid,2)
     1486!--             Terrain-following masked output
     1487                DO  i = 1, mask_size_l(mid,1)
     1488                   DO  j = 1, mask_size_l(mid,2)
    14641489!             
    1465 !--                      Get k index of highest horizontal surface
    1466                          topo_top_ind = get_topography_top_index_ji( &
    1467                                            mask_j(mid,j),  &
    1468                                            mask_i(mid,i),  &
    1469                                            grid                    )
     1490!--                   Get k index of highest horizontal surface
     1491                      topo_top_ind = get_topography_top_index_ji( &
     1492                                        mask_j(mid,j),  &
     1493                                        mask_i(mid,i),  &
     1494                                        grid                    )
    14701495!             
    1471 !--                      Save output array
    1472                          DO  k = 1, mask_size_l(mid,3)
    1473                             local_pf(i,j,k) = chem_species(lsp)%conc_av(  &
    1474                                                  MIN( topo_top_ind+mask_k(mid,k), &
    1475                                                       nzt+1 ),            &
    1476                                                  mask_j(mid,j),           &
    1477                                                  mask_i(mid,i)         )
    1478                          ENDDO
     1496!--                   Save output array
     1497                      DO  k = 1, mask_size_l(mid,3)
     1498                         local_pf(i,j,k) = chem_species(lsp)%conc_av(  &
     1499                                              MIN( topo_top_ind+mask_k(mid,k), &
     1500                                                   nzt+1 ),            &
     1501                                              mask_j(mid,j),           &
     1502                                              mask_i(mid,i)         )
    14791503                      ENDDO
    14801504                   ENDDO
    1481 
    1482                 ENDIF
    1483 
     1505                ENDDO
    14841506
    14851507             ENDIF
    1486              found = .FALSE.
     1508
     1509
    14871510          ENDIF
    1488        ENDDO
    1489 
    1490        RETURN
    1491      
    1492     END SUBROUTINE chem_data_output_mask     
     1511          found = .FALSE.
     1512       ENDIF
     1513    ENDDO
     1514
     1515    RETURN
     1516
     1517 END SUBROUTINE chem_data_output_mask     
    14931518
    14941519!
    … …  
    15001525!> It is called out from subroutine netcdf.
    15011526!------------------------------------------------------------------------------!
    1502     SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
    1503 
    1504        IMPLICIT NONE
    1505 
    1506        CHARACTER (LEN=*), INTENT(IN)  ::  var          !<
    1507        LOGICAL, INTENT(OUT)           ::  found        !<
    1508        CHARACTER (LEN=*), INTENT(OUT) ::  grid_x       !<
    1509        CHARACTER (LEN=*), INTENT(OUT) ::  grid_y       !<
    1510        CHARACTER (LEN=*), INTENT(OUT) ::  grid_z       !<
    1511 
    1512        found  = .TRUE.
    1513 
    1514        IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' )   THEN                   !< always the same grid for chemistry variables
    1515           grid_x = 'x'
    1516           grid_y = 'y'
    1517           grid_z = 'zu'                             
    1518        ELSE
    1519           found  = .FALSE.
    1520           grid_x = 'none'
    1521           grid_y = 'none'
    1522           grid_z = 'none'
    1523        ENDIF
    1524 
    1525 
    1526     END SUBROUTINE chem_define_netcdf_grid
     1527 SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
     1528
     1529    IMPLICIT NONE
     1530
     1531    CHARACTER (LEN=*), INTENT(IN)  ::  var          !<
     1532    LOGICAL, INTENT(OUT)           ::  found        !<
     1533    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x       !<
     1534    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y       !<
     1535    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z       !<
     1536
     1537    found  = .TRUE.
     1538
     1539    IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' )   THEN                   !< always the same grid for chemistry variables
     1540       grid_x = 'x'
     1541       grid_y = 'y'
     1542       grid_z = 'zu'                             
     1543    ELSE
     1544       found  = .FALSE.
     1545       grid_x = 'none'
     1546       grid_y = 'none'
     1547       grid_z = 'none'
     1548    ENDIF
     1549
     1550
     1551 END SUBROUTINE chem_define_netcdf_grid
    15271552!
    15281553!------------------------------------------------------------------------------!
    … …  
    15321557!> Subroutine defining header output for chemistry model
    15331558!------------------------------------------------------------------------------!
    1534     SUBROUTINE chem_header ( io )
    1535        
    1536        IMPLICIT NONE
    1537  
    1538        INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
    1539        INTEGER(iwp)             :: lsp            !< running index for chem spcs
    1540        INTEGER(iwp)             :: cs_fixed
    1541        CHARACTER (LEN=80)       :: docsflux_chr
    1542        CHARACTER (LEN=80)       :: docsinit_chr
    1543 
    1544 !
    1545 !--    Write chemistry model  header
    1546        WRITE( io, 1 )
    1547 
    1548 !--    Gasphase reaction status
    1549        IF ( chem_gasphase_on )  THEN
    1550           WRITE( io, 2 )
    1551        ELSE
    1552           WRITE( io, 3 )
     1559 SUBROUTINE chem_header( io )
     1560
     1561    IMPLICIT NONE
     1562
     1563    INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
     1564    INTEGER(iwp)             :: lsp            !< running index for chem spcs
     1565    INTEGER(iwp)             :: cs_fixed
     1566    CHARACTER (LEN=80)       :: docsflux_chr
     1567    CHARACTER (LEN=80)       :: docsinit_chr
     1568
     1569!
     1570!-- Write chemistry model  header
     1571    WRITE( io, 1 )
     1572
     1573!-- Gasphase reaction status
     1574    IF ( chem_gasphase_on )  THEN
     1575       WRITE( io, 2 )
     1576    ELSE
     1577       WRITE( io, 3 )
     1578    ENDIF
     1579
     1580!
     1581!-- Chemistry time-step
     1582    WRITE ( io, 4 ) cs_time_step
     1583
     1584!-- Emission mode info
     1585    IF ( mode_emis == "DEFAULT" )  THEN
     1586       WRITE( io, 5 )
     1587    ELSEIF ( mode_emis == "PARAMETERIZED" )  THEN
     1588       WRITE( io, 6 )
     1589    ELSEIF ( mode_emis == "PRE-PROCESSED" )  THEN
     1590       WRITE( io, 7 )
     1591    ENDIF
     1592
     1593!-- Photolysis scheme info
     1594    IF ( photolysis_scheme == "simple" )      THEN
     1595       WRITE( io, 8 )
     1596    ELSEIF (photolysis_scheme == "constant" ) THEN
     1597       WRITE( io, 9 )
     1598    ENDIF
     1599
     1600!-- Emission flux info
     1601    lsp = 1
     1602    docsflux_chr ='Chemical species for surface emission flux: '
     1603    DO WHILE ( surface_csflux_name(lsp) /= 'novalue' )
     1604       docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ',' 
     1605       IF ( LEN_TRIM( docsflux_chr ) >= 75 )  THEN
     1606        WRITE ( io, 10 )  docsflux_chr
     1607        docsflux_chr = '       '
    15531608       ENDIF
    1554 
    1555 !      Chemistry time-step
    1556        WRITE ( io, 4 ) cs_time_step
    1557 
    1558 !--    Emission mode info
    1559        IF ( mode_emis == "DEFAULT" )  THEN
    1560           WRITE( io, 5 )
    1561        ELSEIF ( mode_emis == "PARAMETERIZED" )  THEN
    1562           WRITE( io, 6 )
    1563        ELSEIF ( mode_emis == "PRE-PROCESSED" )  THEN
    1564           WRITE( io, 7 )
    1565        ENDIF
    1566 
    1567 !--    Photolysis scheme info
    1568        IF ( photolysis_scheme == "simple" )      THEN
    1569           WRITE( io, 8 )
    1570        ELSEIF (photolysis_scheme == "conastant" ) THEN
    1571           WRITE( io, 9 )
     1609       lsp = lsp + 1
     1610    ENDDO
     1611
     1612    IF ( docsflux_chr /= '' )  THEN
     1613       WRITE ( io, 10 )  docsflux_chr
     1614    ENDIF
     1615
     1616
     1617!-- initializatoin of Surface and profile chemical species
     1618
     1619    lsp = 1
     1620    docsinit_chr ='Chemical species for initial surface and profile emissions: '
     1621    DO WHILE ( cs_name(lsp) /= 'novalue' )
     1622       docsinit_chr = TRIM( docsinit_chr ) // ' ' // TRIM( cs_name(lsp) ) // ',' 
     1623       IF ( LEN_TRIM( docsinit_chr ) >= 75 )  THEN
     1624        WRITE ( io, 11 )  docsinit_chr
     1625        docsinit_chr = '       '
    15721626       ENDIF
    1573  
    1574  !--   Emission flux info
    1575        lsp = 1
    1576        docsflux_chr ='Chemical species for surface emission flux: '
    1577        DO WHILE ( surface_csflux_name(lsp) /= 'novalue' )
    1578           docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ',' 
    1579           IF ( LEN_TRIM( docsflux_chr ) >= 75 )  THEN
    1580            WRITE ( io, 10 )  docsflux_chr
    1581            docsflux_chr = '       '
    1582           ENDIF
    1583           lsp = lsp + 1
    1584        ENDDO
    1585  
    1586        IF ( docsflux_chr /= '' )  THEN
    1587           WRITE ( io, 10 )  docsflux_chr
    1588        ENDIF
    1589 
    1590 
    1591 !--    initializatoin of Surface and profile chemical species
    1592 
    1593        lsp = 1
    1594        docsinit_chr ='Chemical species for initial surface and profile emissions: '
    1595        DO WHILE ( cs_name(lsp) /= 'novalue' )
    1596           docsinit_chr = TRIM( docsinit_chr ) // ' ' // TRIM( cs_name(lsp) ) // ',' 
    1597           IF ( LEN_TRIM( docsinit_chr ) >= 75 )  THEN
    1598            WRITE ( io, 11 )  docsinit_chr
    1599            docsinit_chr = '       '
    1600           ENDIF
    1601           lsp = lsp + 1
    1602        ENDDO
    1603  
    1604        IF ( docsinit_chr /= '' )  THEN
    1605           WRITE ( io, 11 )  docsinit_chr
    1606        ENDIF
    1607 
    1608 !--    number of variable and fix chemical species and number of reactions
    1609        cs_fixed = nspec - nvar
    1610        WRITE ( io, * ) '   --> Chemical species, variable: ', nvar
    1611        WRITE ( io, * ) '   --> Chemical species, fixed   : ', cs_fixed
    1612        WRITE ( io, * ) '   --> Total number of reactions : ', nreact
     1627       lsp = lsp + 1
     1628    ENDDO
     1629
     1630    IF ( docsinit_chr /= '' )  THEN
     1631       WRITE ( io, 11 )  docsinit_chr
     1632    ENDIF
     1633
     1634!-- number of variable and fix chemical species and number of reactions
     1635    cs_fixed = nspec - nvar
     1636    WRITE ( io, * ) '   --> Chemical species, variable: ', nvar
     1637    WRITE ( io, * ) '   --> Chemical species, fixed   : ', cs_fixed
     1638    WRITE ( io, * ) '   --> Total number of reactions : ', nreact
    16131639
    16141640
    161516411   FORMAT (//' Chemistry model information:'/                                  &
    1616               ' ----------------------------'/)
     1642           ' ----------------------------'/)
    161716432   FORMAT ('    --> Chemical reactions are turned on')
    161816443   FORMAT ('    --> Chemical reactions are turned off')
    … …  
    16271653!
    16281654!
    1629     END SUBROUTINE chem_header
     1655 END SUBROUTINE chem_header
    16301656
    16311657!
    … …  
    16361662!> Subroutine initializating chemistry_model_mod
    16371663!------------------------------------------------------------------------------!
    1638     SUBROUTINE chem_init
    1639 
    1640        USE control_parameters,                                                  &
    1641           ONLY: message_string, io_blocks, io_group, turbulent_inflow
    1642        USE arrays_3d,                                                           &
    1643            ONLY: mean_inflow_profiles
    1644        USE pegrid
    1645 
    1646        IMPLICIT   none
    1647    !-- local variables
    1648        INTEGER                 :: i,j               !< running index for for horiz numerical grid points
    1649        INTEGER                 :: lsp               !< running index for chem spcs
    1650        INTEGER                 :: lpr_lev           !< running index for chem spcs profile level
     1664 SUBROUTINE chem_init
     1665
     1666    USE control_parameters,                                                  &
     1667        ONLY:  message_string, io_blocks, io_group, turbulent_inflow
     1668    USE arrays_3d,                                                           &
     1669        ONLY:  mean_inflow_profiles
     1670    USE pegrid
     1671
     1672    IMPLICIT NONE
     1673!-- local variables
     1674    INTEGER(iwp) ::  i                 !< running index for for horiz numerical grid points
     1675    INTEGER(iwp) ::  j                 !< running index for for horiz numerical grid points
     1676    INTEGER(iwp) ::  lsp               !< running index for chem spcs
     1677    INTEGER(iwp) ::  lpr_lev           !< running index for chem spcs profile level
    16511678!
    16521679!-- NOPOINTER version not implemented yet
    … …  
    16571684!
    16581685!-- Allocate memory for chemical species
    1659        ALLOCATE( chem_species(nspec) )
    1660        ALLOCATE( spec_conc_1 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
    1661        ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
    1662        ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
    1663        ALLOCATE( spec_conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
    1664        ALLOCATE( phot_frequen(nphot) )
    1665        ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) )
    1666        ALLOCATE( bc_cs_t_val(nspec) )
     1686    ALLOCATE( chem_species(nspec) )
     1687    ALLOCATE( spec_conc_1 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
     1688    ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
     1689    ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
     1690    ALLOCATE( spec_conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
     1691    ALLOCATE( phot_frequen(nphot) )
     1692    ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) )
     1693    ALLOCATE( bc_cs_t_val(nspec) )
    16671694!
    16681695!-- Initialize arrays
    1669        spec_conc_1 (:,:,:,:) = 0.0_wp
    1670        spec_conc_2 (:,:,:,:) = 0.0_wp
    1671        spec_conc_3 (:,:,:,:) = 0.0_wp
    1672        spec_conc_av(:,:,:,:) = 0.0_wp
    1673 
    1674 
    1675        DO lsp = 1, nspec
    1676           chem_species(lsp)%name    = spc_names(lsp)
    1677 
    1678           chem_species(lsp)%conc   (nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_1 (:,:,:,lsp)
    1679           chem_species(lsp)%conc_p (nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_2 (:,:,:,lsp)
    1680           chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_3 (:,:,:,lsp)
    1681           chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_av(:,:,:,lsp)     
    1682 
    1683           ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))                   
    1684           chem_species(lsp)%cssws_av    = 0.0_wp
    1685 !
    1686 !--       The following block can be useful when emission module is not applied. &
    1687 !--       if emission module is applied the following block will be overwritten.
    1688           ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))   
    1689           ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))   
    1690           ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
    1691           ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))   
    1692           chem_species(lsp)%flux_s_cs = 0.0_wp                                     
    1693           chem_species(lsp)%flux_l_cs = 0.0_wp                                     
    1694           chem_species(lsp)%diss_s_cs = 0.0_wp                                     
    1695           chem_species(lsp)%diss_l_cs = 0.0_wp                                     
     1696    spec_conc_1 (:,:,:,:) = 0.0_wp
     1697    spec_conc_2 (:,:,:,:) = 0.0_wp
     1698    spec_conc_3 (:,:,:,:) = 0.0_wp
     1699    spec_conc_av(:,:,:,:) = 0.0_wp
     1700
     1701
     1702    DO lsp = 1, nspec
     1703       chem_species(lsp)%name    = spc_names(lsp)
     1704
     1705       chem_species(lsp)%conc   (nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_1 (:,:,:,lsp)
     1706       chem_species(lsp)%conc_p (nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_2 (:,:,:,lsp)
     1707       chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_3 (:,:,:,lsp)
     1708       chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_av(:,:,:,lsp)     
     1709
     1710       ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))                   
     1711       chem_species(lsp)%cssws_av    = 0.0_wp
     1712!
     1713!--    The following block can be useful when emission module is not applied. &
     1714!--    if emission module is applied the following block will be overwritten.
     1715       ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))   
     1716       ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))   
     1717       ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))
     1718       ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))   
     1719       chem_species(lsp)%flux_s_cs = 0.0_wp                                     
     1720       chem_species(lsp)%flux_l_cs = 0.0_wp                                     
     1721       chem_species(lsp)%diss_s_cs = 0.0_wp                                     
     1722       chem_species(lsp)%diss_l_cs = 0.0_wp                                     
    16961723!
    16971724!--   Allocate memory for initial concentration profiles
    … …  
    17011728!--               We have to find another solution since chem_init should
    17021729!--               eventually be called from init_3d_model!!)
    1703           ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) )
    1704           chem_species(lsp)%conc_pr_init(:) = 0.0_wp
    1705 
    1706        ENDDO
    1707 
    1708 !
    1709 !--    Initial concentration of profiles is prescribed by parameters cs_profile
    1710 !--    and cs_heights in the namelist &chemistry_parameters
    1711        CALL chem_init_profiles     
    1712            
    1713            
     1730       ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) )
     1731       chem_species(lsp)%conc_pr_init(:) = 0.0_wp
     1732
     1733    ENDDO
     1734
     1735!
     1736!-- Initial concentration of profiles is prescribed by parameters cs_profile
     1737!-- and cs_heights in the namelist &chemistry_parameters
     1738    CALL chem_init_profiles     
     1739
     1740
    17141741!
    17151742!-- Initialize model variables
    17161743
    17171744
    1718        IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
    1719             TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
     1745    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
     1746         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
    17201747
    17211748
    17221749!--    First model run of a possible job queue.
    17231750!--    Initial profiles of the variables must be computed.
    1724           IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
    1725             CALL location_message( 'initializing with 1D chemistry model profiles', .FALSE. )
    1726 !
    1727 !--        Transfer initial profiles to the arrays of the 3D model
    1728              DO lsp = 1, nspec
    1729                 DO  i = nxlg, nxrg
    1730                    DO  j = nysg, nyng
    1731                       DO lpr_lev = 1, nz + 1
    1732                          chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)
    1733                       ENDDO
    1734                    ENDDO
    1735                 ENDDO   
    1736              ENDDO
    1737          
    1738           ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
    1739           THEN
    1740              CALL location_message( 'initializing with constant chemistry profiles', .FALSE. )
    1741 
    1742              DO lsp = 1, nspec
    1743                 DO  i = nxlg, nxrg
    1744                    DO  j = nysg, nyng
    1745                       chem_species(lsp)%conc(:,j,i) = chem_species(lsp)%conc_pr_init   
     1751       IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
     1752         CALL location_message( 'initializing with 1D chemistry model profiles', .FALSE. )
     1753!
     1754!--       Transfer initial profiles to the arrays of the 3D model
     1755          DO lsp = 1, nspec
     1756             DO  i = nxlg, nxrg
     1757                DO  j = nysg, nyng
     1758                   DO lpr_lev = 1, nz + 1
     1759                      chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)
    17461760                   ENDDO
    17471761                ENDDO
     1762             ENDDO   
     1763          ENDDO
     1764
     1765       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
     1766       THEN
     1767          CALL location_message( 'initializing with constant chemistry profiles', .FALSE. )
     1768
     1769          DO lsp = 1, nspec
     1770             DO  i = nxlg, nxrg
     1771                DO  j = nysg, nyng
     1772                   chem_species(lsp)%conc(:,j,i) = chem_species(lsp)%conc_pr_init   
     1773                ENDDO
    17481774             ENDDO
    1749 
    1750           ENDIF
    1751 
    1752 !
    1753 !--       If required, change the surface chem spcs at the start of the 3D run
    1754           IF ( cs_surface_initial_change(1) /= 0.0_wp ) THEN           
    1755              DO lsp = 1, nspec
    1756                 chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) +  &
    1757                                                   cs_surface_initial_change(lsp)
    1758              ENDDO
    1759           ENDIF
    1760 !
    1761 !--      Initiale old and new time levels.
    1762           DO lsp = 1, nvar
    1763              chem_species(lsp)%tconc_m = 0.0_wp                     
    1764              chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
    17651775          ENDDO
    17661776
    17671777       ENDIF
    17681778
    1769 
    1770 
    1771 !---   new code add above this line
    1772        DO lsp = 1, nphot
    1773           phot_frequen(lsp)%name = phot_names(lsp)
    1774 !
    1775 !--   todo: remove or replace by "CALL message" mechanism (kanani)
     1779!
     1780!--    If required, change the surface chem spcs at the start of the 3D run
     1781       IF ( cs_surface_initial_change(1) /= 0.0_wp ) THEN           
     1782          DO lsp = 1, nspec
     1783             chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) +  &
     1784                                               cs_surface_initial_change(lsp)
     1785          ENDDO
     1786       ENDIF
     1787!
     1788!--    Initiale old and new time levels.
     1789       DO lsp = 1, nvar
     1790          chem_species(lsp)%tconc_m = 0.0_wp                     
     1791          chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
     1792       ENDDO
     1793
     1794    ENDIF
     1795
     1796
     1797
     1798!--- new code add above this line
     1799    DO lsp = 1, nphot
     1800       phot_frequen(lsp)%name = phot_names(lsp)
     1801!
     1802!-- todo: remove or replace by "CALL message" mechanism (kanani)
    17761803!         IF( myid == 0 )  THEN
    1777 !            WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,trim(phot_names(lsp))
     1804!            WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,TRIM(phot_names(lsp))
    17781805!         ENDIF
    1779              phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => freq_1(:,:,:,lsp)
    1780        ENDDO
    1781 
    1782        RETURN
     1806          phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => freq_1(:,:,:,lsp)
     1807    ENDDO
     1808
     1809    RETURN
    17831810
    17841811 END SUBROUTINE chem_init
    … …  
    17931820!> analogue to parameters u_profile, v_profile and uv_heights)
    17941821!------------------------------------------------------------------------------!
    1795     SUBROUTINE chem_init_profiles              !< SUBROUTINE is called from chem_init in case of
    1796                                                !< TRIM( initializing_actions ) /= 'read_restart_data'
    1797                                                !< We still need to see what has to be done in case of restart run
    1798        USE chem_modules
    1799 
    1800        IMPLICIT NONE
    1801 
    1802    !-- Local variables
    1803        INTEGER ::  lsp        !< running index for number of species in derived data type species_def
    1804        INTEGER ::  lsp_usr     !< running index for number of species (user defined)  in cs_names, cs_profiles etc
    1805        INTEGER ::  lpr_lev    !< running index for profile level for each chem spcs.
    1806        INTEGER ::  npr_lev    !< the next available profile lev
    1807 
    1808 !-----------------
    1809 !-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles
    1810 !-- and heights. If parameter "cs_profile" is not prescribed then initial surface values
    1811 !-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and
    1812 !-- "cs_heights" are prescribed, their values will!override the constant profile given by
    1813 !-- "cs_surface".
    1814 !
    1815        IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
    1816           lsp_usr = 1
    1817           DO  WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )   !'novalue' is the default
    1818              DO  lsp = 1, nspec                                !
    1819 !--             create initial profile (conc_pr_init) for each chemical species
    1820                 IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN   !
    1821                    IF ( cs_profile(lsp_usr,1) == 9999999.9_wp ) THEN
    1822 !--                set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
    1823                       DO lpr_lev = 0, nzt+1
    1824                          chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
    1825                       ENDDO
    1826                    ELSE
    1827                       IF ( cs_heights(1,1) /= 0.0_wp )  THEN
    1828                          message_string = 'The surface value of cs_heights must be 0.0'
    1829                          CALL message( 'chem_check_parameters', 'CM0434', 1, 2, 0, 6, 0 )
     1822 SUBROUTINE chem_init_profiles              !< SUBROUTINE is called from chem_init in case of
     1823   !< TRIM( initializing_actions ) /= 'read_restart_data'
     1824   !< We still need to see what has to be done in case of restart run
     1825    USE chem_modules
     1826
     1827    IMPLICIT NONE
     1828
     1829    !-- Local variables
     1830    INTEGER ::  lsp        !< running index for number of species in derived data type species_def
     1831    INTEGER ::  lsp_usr     !< running index for number of species (user defined)  in cs_names, cs_profiles etc
     1832    INTEGER ::  lpr_lev    !< running index for profile level for each chem spcs.
     1833    INTEGER ::  npr_lev    !< the next available profile lev
     1834
     1835    !-----------------
     1836    !-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles
     1837    !-- and heights. If parameter "cs_profile" is not prescribed then initial surface values
     1838    !-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and
     1839    !-- "cs_heights" are prescribed, their values will!override the constant profile given by
     1840    !-- "cs_surface".
     1841    !
     1842    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
     1843       lsp_usr = 1
     1844       DO  WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )   !'novalue' is the default
     1845          DO  lsp = 1, nspec                                !
     1846             !-- create initial profile (conc_pr_init) for each chemical species
     1847             IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN   !
     1848                IF ( cs_profile(lsp_usr,1) == 9999999.9_wp ) THEN
     1849                   !-- set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
     1850                   DO lpr_lev = 0, nzt+1
     1851                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
     1852                   ENDDO
     1853                ELSE
     1854                   IF ( cs_heights(1,1) /= 0.0_wp )  THEN
     1855                      message_string = 'The surface value of cs_heights must be 0.0'
     1856                      CALL message( 'chem_check_parameters', 'CM0434', 1, 2, 0, 6, 0 )
     1857                   ENDIF
     1858
     1859                   use_prescribed_profile_data = .TRUE.
     1860
     1861                   npr_lev = 1
     1862                   !                     chem_species(lsp)%conc_pr_init(0) = 0.0_wp
     1863                   DO  lpr_lev = 1, nz+1
     1864                      IF ( npr_lev < 100 )  THEN
     1865                         DO  WHILE ( cs_heights(lsp_usr, npr_lev+1) <= zu(lpr_lev) )
     1866                            npr_lev = npr_lev + 1
     1867                            IF ( npr_lev == 100 )  THEN
     1868                               message_string = 'number of chem spcs exceeding the limit'
     1869                               CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )               
     1870                               EXIT
     1871                            ENDIF
     1872                         ENDDO
    18301873                      ENDIF
    1831      
    1832                       use_prescribed_profile_data = .TRUE.
    1833    
    1834                       npr_lev = 1
    1835 !                     chem_species(lsp)%conc_pr_init(0) = 0.0_wp
    1836                       DO  lpr_lev = 1, nz+1
    1837                          IF ( npr_lev < 100 )  THEN
    1838                             DO  WHILE ( cs_heights(lsp_usr, npr_lev+1) <= zu(lpr_lev) )
    1839                                npr_lev = npr_lev + 1
    1840                                IF ( npr_lev == 100 )  THEN
    1841                                    message_string = 'number of chem spcs exceeding the limit'
    1842                                    CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )               
    1843                                    EXIT
    1844                                ENDIF
    1845                             ENDDO
    1846                          ENDIF
    1847                          IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr, npr_lev + 1) /= 9999999.9_wp )  THEN
    1848                             chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +                          &
    1849                                                     ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                          &
    1850                                                     ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) *  &
    1851                                                     ( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
    1852                          ELSE
    1853                                chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev)
    1854                          ENDIF
    1855                       ENDDO
    1856                    ENDIF
    1857 !-- If a profile is prescribed explicity using cs_profiles and cs_heights, then 
    1858 !-- chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
    1859 !-- on the cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:).
     1874                      IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr, npr_lev + 1) /= 9999999.9_wp )  THEN
     1875                         chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +       &
     1876                              ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                          &
     1877                              ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) *  &
     1878                              ( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
     1879                      ELSE
     1880                         chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev)
     1881                      ENDIF
     1882                   ENDDO
    18601883                ENDIF
     1884                !-- If a profile is prescribed explicity using cs_profiles and cs_heights, then 
     1885                !-- chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
     1886                !-- on the cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:).
     1887             ENDIF
     1888          ENDDO
     1889          lsp_usr = lsp_usr + 1
     1890       ENDDO
     1891    ENDIF
     1892
     1893 END SUBROUTINE chem_init_profiles
     1894 !
     1895 !------------------------------------------------------------------------------!
     1896 !
     1897 ! Description:
     1898 ! ------------
     1899 !> Subroutine to integrate chemical species in the given chemical mechanism
     1900 !------------------------------------------------------------------------------!
     1901
     1902 SUBROUTINE chem_integrate_ij( i, j )
     1903
     1904    USE cpulog,                                                              &
     1905         ONLY:  cpu_log, log_point
     1906    USE statistics,                                                          &
     1907         ONLY:  weight_pres
     1908    USE control_parameters,                                                 &
     1909         ONLY:  dt_3d, intermediate_timestep_count, simulated_time
     1910
     1911    IMPLICIT NONE
     1912    INTEGER,INTENT(IN)       :: i
     1913    INTEGER,INTENT(IN)       :: j
     1914
     1915    !--   local variables
     1916    INTEGER(iwp) ::  lsp                                                     !< running index for chem spcs.
     1917    INTEGER(iwp) ::  lph                                                     !< running index for photolysis frequencies
     1918    INTEGER      ::  k
     1919    INTEGER      ::  m
     1920    INTEGER      ::  istatf
     1921    INTEGER, DIMENSION(20)    :: istatus
     1922    REAL(kind=wp), DIMENSION(nzb+1:nzt,nspec)                :: tmp_conc
     1923    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_temp
     1924    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_qvap
     1925    REAL(kind=wp), DIMENSION(nzb+1:nzt,nphot)                :: tmp_phot
     1926    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact
     1927    REAL(kind=wp), DIMENSION(nzb+1:nzt)                      :: tmp_fact_i    !< conversion factor between
     1928                                                                              !<    molecules cm^{-3} and ppm
     1929
     1930    INTEGER,DIMENSION(nzb+1:nzt)                            :: nacc          !< Number of accepted steps
     1931    INTEGER,DIMENSION(nzb+1:nzt)                            :: nrej          !< Number of rejected steps
     1932
     1933    REAL(wp)                         ::  conv                                !< conversion factor
     1934    REAL(wp), PARAMETER              ::  ppm2fr  = 1.0e-6_wp                 !< Conversion factor ppm to fraction
     1935!    REAL(wp), PARAMETER              ::  xm_air  = 28.96_wp                  !< Mole mass of dry air
     1936!    REAL(wp), PARAMETER              ::  xm_h2o  = 18.01528_wp               !< Mole mass of water vapor
     1937    REAL(wp), PARAMETER              ::  pref_i  = 1._wp / 100000.0_wp       !< inverse reference pressure (1/Pa)
     1938    REAL(wp), PARAMETER              ::  t_std   = 273.15_wp                 !< standard pressure (Pa)
     1939    REAL(wp), PARAMETER              ::  p_std   = 101325.0_wp               !< standard pressure (Pa)
     1940    REAL(wp), PARAMETER              ::  vmolcm  = 22.414e3_wp               !< Mole volume (22.414 l) in cm^{-3}
     1941    REAL(wp), PARAMETER              ::  xna     = 6.022e23_wp               !< Avogadro number (molecules/mol)
     1942
     1943    REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local
     1944
     1945
     1946    REAL(kind=wp)  :: dt_chem                                             
     1947
     1948    CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'start' )
     1949    !<     set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry
     1950    IF (chem_gasphase_on) THEN
     1951       nacc = 0
     1952       nrej = 0
     1953
     1954       tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt)
     1955       !
     1956       !--       ppm to molecules/cm**3
     1957       !--       tmp_fact = 1.e-6_wp*6.022e23_wp/(22.414_wp*1000._wp) * 273.15_wp * hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp ) 
     1958       conv = ppm2fr * xna / vmolcm
     1959       tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std)
     1960       tmp_fact_i = 1.0_wp/tmp_fact
     1961
     1962       IF ( humidity ) THEN
     1963          IF ( bulk_cloud_model )  THEN
     1964             tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) * xm_air/xm_h2o * tmp_fact(:)
     1965          ELSE
     1966             tmp_qvap(:) = q(nzb+1:nzt,j,i) * xm_air/xm_h2o * tmp_fact(:)
     1967          ENDIF
     1968       ELSE
     1969          tmp_qvap(:) = 0.01 * tmp_fact(:)                          !< Constant value for q if water vapor is not computed
     1970       ENDIF
     1971
     1972       DO lsp = 1,nspec
     1973          tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:)
     1974       ENDDO
     1975
     1976       DO lph = 1,nphot
     1977          tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)               
     1978       ENDDO
     1979       !
     1980       !--   todo: remove (kanani)
     1981       !           IF(myid == 0 .AND. chem_debug0 ) THEN
     1982       !              IF (i == 10 .AND. j == 10) WRITE(0,*) 'begin chemics step ',dt_3d
     1983       !           ENDIF
     1984
     1985       !--       Compute length of time step
     1986       IF ( call_chem_at_all_substeps )  THEN
     1987          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
     1988       ELSE
     1989          dt_chem = dt_3d
     1990       ENDIF
     1991
     1992       cs_time_step = dt_chem
     1993
     1994       CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'start' )
     1995
     1996       IF(maxval(rcntrl) > 0.0)   THEN    ! Only if rcntrl is set
     1997          IF( simulated_time <= 2*dt_3d)  THEN
     1998             rcntrl_local = 0
     1999             !
     2000             !--   todo: remove (kanani)
     2001             !                  WRITE(9,'(a,2f10.3)') 'USE Default rcntrl in the first steps ',simulated_time,dt_3d
     2002          ELSE
     2003             rcntrl_local = rcntrl
     2004          ENDIF
     2005       ELSE
     2006          rcntrl_local = 0
     2007       END IF
     2008
     2009       CALL chem_gasphase_integrate ( dt_chem, tmp_conc, tmp_temp, tmp_qvap, tmp_fact, tmp_phot, &
     2010            icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus )
     2011
     2012       CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'stop' )
     2013
     2014       DO lsp = 1,nspec
     2015          chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:)
     2016       ENDDO
     2017
     2018
     2019    ENDIF
     2020    CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'stop' )
     2021
     2022    RETURN
     2023 END SUBROUTINE chem_integrate_ij
     2024 !
     2025 !------------------------------------------------------------------------------!
     2026 !
     2027 ! Description:
     2028 ! ------------
     2029 !> Subroutine defining parin for &chemistry_parameters for chemistry model
     2030 !------------------------------------------------------------------------------!
     2031 SUBROUTINE chem_parin
     2032
     2033    USE chem_modules
     2034    USE control_parameters
     2035
     2036    USE kinds
     2037    USE pegrid
     2038    USE statistics
     2039
     2040    IMPLICIT NONE
     2041
     2042    CHARACTER (LEN=80) ::  line                        !< dummy string that contains the current line of the parameter file
     2043    CHARACTER (LEN=3)  ::  cs_prefix
     2044
     2045    REAL(wp), DIMENSION(nmaxfixsteps) ::   my_steps    !< List of fixed timesteps   my_step(1) = 0.0 automatic stepping
     2046    INTEGER(iwp) ::  i                                 !<
     2047    INTEGER(iwp) ::  j                                 !<
     2048    INTEGER(iwp) ::  max_pr_cs_tmp                     !<
     2049
     2050
     2051    NAMELIST /chemistry_parameters/  bc_cs_b,                          &
     2052         bc_cs_t,                          &
     2053         call_chem_at_all_substeps,        &
     2054         chem_debug0,                      &
     2055         chem_debug1,                      &
     2056         chem_debug2,                      &
     2057         chem_gasphase_on,                 &
     2058         cs_heights,                       &
     2059         cs_name,                          &
     2060         cs_profile,                       &
     2061         cs_surface,                       &
     2062         decycle_chem_lr,                  &
     2063         decycle_chem_ns,                  &           
     2064         decycle_method,                   &
     2065         do_depo,                          &
     2066         emiss_factor_main,                &
     2067         emiss_factor_side,                &                     
     2068         icntrl,                           &
     2069         main_street_id,                   &
     2070         max_street_id,                    &
     2071         my_steps,                         &
     2072         nest_chemistry,                   &
     2073         rcntrl,                           &
     2074         side_street_id,                   &
     2075         photolysis_scheme,                &
     2076         wall_csflux,                      &
     2077         cs_vertical_gradient,             &
     2078         top_csflux,                       &
     2079         surface_csflux,                   &
     2080         surface_csflux_name,              &
     2081         cs_surface_initial_change,        &
     2082         cs_vertical_gradient_level,       &
     2083         !                                       namelist parameters for emissions
     2084         mode_emis,                        &
     2085         time_fac_type,                    &
     2086         daytype_mdh,                      &
     2087         do_emis                             
     2088
     2089    !-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
     2090    !-- so this way we could prescribe a specific flux value for each species
     2091    !>  chemistry_parameters for initial profiles
     2092    !>  cs_names = 'O3', 'NO2', 'NO', ...   to set initial profiles)
     2093    !>  cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3)
     2094    !>  cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.)
     2095    !>  cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, .....  (chem spcs conc at height lvls chem_heights(1,:)) etc.
     2096    !>  If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)"
     2097    !>  then write these cs_surface values to chem_species(lsp)%conc_pr_init(:)
     2098
     2099    !
     2100    !--   Read chem namelist   
     2101
     2102    INTEGER             :: ier
     2103    CHARACTER(LEN=64)   :: text
     2104    CHARACTER(LEN=8)    :: solver_type
     2105
     2106    icntrl    = 0
     2107    rcntrl    = 0.0_wp
     2108    my_steps  = 0.0_wp
     2109    photolysis_scheme = 'simple'
     2110    atol = 1.0_wp
     2111    rtol = 0.01_wp
     2112    !
     2113    !--   Try to find chemistry package
     2114    REWIND ( 11 )
     2115    line = ' '
     2116    DO   WHILE ( INDEX( line, '&chemistry_parameters' ) == 0 )
     2117       READ ( 11, '(A)', END=20 )  line
     2118    ENDDO
     2119    BACKSPACE ( 11 )
     2120    !
     2121    !--   Read chemistry namelist
     2122    READ ( 11, chemistry_parameters, ERR = 10, END = 20 )     
     2123    !
     2124    !--   Enable chemistry model
     2125    air_chemistry = .TRUE.                   
     2126    GOTO 20
     2127
     2128 10 BACKSPACE( 11 )
     2129    READ( 11 , '(A)') line
     2130    CALL parin_fail_message( 'chemistry_parameters', line )
     2131
     2132 20 CONTINUE
     2133
     2134    !
     2135    !--    check for  emission mode for chem species
     2136    IF ( (mode_emis /= 'PARAMETERIZED')  .AND. ( mode_emis /= 'DEFAULT' ) .AND. ( mode_emis /= 'PRE-PROCESSED'  ) )  THEN
     2137       message_string = 'Incorrect mode_emiss  option select. Please check spelling'
     2138       CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
     2139    ENDIF
     2140
     2141    t_steps = my_steps         
     2142
     2143    !--    Determine the number of user-defined profiles and append them to the
     2144    !--    standard data output (data_output_pr)
     2145    max_pr_cs_tmp = 0
     2146    i = 1
     2147
     2148    DO  WHILE ( data_output_pr(i)  /= ' '  .AND.  i <= 100 )
     2149       IF ( TRIM( data_output_pr(i)(1:3)) == 'kc_' ) THEN
     2150          max_pr_cs_tmp = max_pr_cs_tmp+1
     2151       ENDIF
     2152       i = i +1
     2153    ENDDO
     2154
     2155    IF ( max_pr_cs_tmp > 0 ) THEN
     2156       cs_pr_namelist_found = .TRUE.
     2157       max_pr_cs = max_pr_cs_tmp
     2158    ENDIF
     2159
     2160    !      Set Solver Type
     2161    IF(icntrl(3) == 0)   THEN
     2162       solver_type = 'rodas3'           !Default
     2163    ELSE IF(icntrl(3) == 1)   THEN
     2164       solver_type = 'ros2'
     2165    ELSE IF(icntrl(3) == 2)   THEN
     2166       solver_type = 'ros3'
     2167    ELSE IF(icntrl(3) == 3)   THEN
     2168       solver_type = 'ro4'
     2169    ELSE IF(icntrl(3) == 4)   THEN
     2170       solver_type = 'rodas3'
     2171    ELSE IF(icntrl(3) == 5)   THEN
     2172       solver_type = 'rodas4'
     2173    ELSE IF(icntrl(3) == 6)   THEN
     2174       solver_type = 'Rang3'
     2175    ELSE
     2176       message_string = 'illegal solver type'
     2177       CALL message( 'chem_parin', 'PA0506', 1, 2, 0, 6, 0 )
     2178    END IF
     2179
     2180    !
     2181    !--   todo: remove or replace by "CALL message" mechanism (kanani)
     2182    !       write(text,*) 'gas_phase chemistry: solver_type = ',TRIM(solver_type)
     2183    !kk    Has to be changed to right calling sequence
     2184    !kk       CALL location_message( TRIM(text), .FALSE. )
     2185    !        IF(myid == 0)   THEN
     2186    !           write(9,*) ' '
     2187    !           write(9,*) 'kpp setup '
     2188    !           write(9,*) ' '
     2189    !           write(9,*) '    gas_phase chemistry: solver_type = ',TRIM(solver_type)
     2190    !           write(9,*) ' '
     2191    !           write(9,*) '    Hstart  = ',rcntrl(3)
     2192    !           write(9,*) '    FacMin  = ',rcntrl(4)
     2193    !           write(9,*) '    FacMax  = ',rcntrl(5)
     2194    !           write(9,*) ' '
     2195    !           IF(vl_dim > 1)    THEN
     2196    !              write(9,*) '    Vector mode                   vektor length = ',vl_dim
     2197    !           ELSE
     2198    !              write(9,*) '    Scalar mode'
     2199    !           ENDIF
     2200    !           write(9,*) ' '
     2201    !        END IF
     2202
     2203    RETURN
     2204
     2205 END SUBROUTINE chem_parin
     2206
     2207 !
     2208 !------------------------------------------------------------------------------!
     2209 !
     2210 ! Description:
     2211 ! ------------
     2212 !> Subroutine calculating prognostic equations for chemical species
     2213 !> (vector-optimized).
     2214 !> Routine is called separately for each chemical species over a loop from
     2215 !> prognostic_equations.
     2216 !------------------------------------------------------------------------------!
     2217 SUBROUTINE chem_prognostic_equations( cs_scalar_p, cs_scalar, tcs_scalar_m,   &
     2218      pr_init_cs, ilsp )
     2219
     2220    USE advec_s_pw_mod,                                                        &
     2221         ONLY:  advec_s_pw
     2222    USE advec_s_up_mod,                                                        &
     2223         ONLY:  advec_s_up
     2224    USE advec_ws,                                                              &
     2225         ONLY:  advec_s_ws
     2226    USE diffusion_s_mod,                                                       &
     2227         ONLY:  diffusion_s
     2228    USE indices,                                                               &
     2229         ONLY:  nxl, nxr, nyn, nys, wall_flags_0
     2230    USE pegrid
     2231    USE surface_mod,                                                           &
     2232         ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
     2233         surf_usm_v
     2234
     2235    IMPLICIT NONE
     2236
     2237    INTEGER ::  i   !< running index
     2238    INTEGER ::  j   !< running index
     2239    INTEGER ::  k   !< running index
     2240
     2241    INTEGER(iwp),INTENT(IN) ::  ilsp          !<
     2242
     2243    REAL(wp), DIMENSION(0:nz+1) ::  pr_init_cs   !<
     2244
     2245    REAL(wp), DIMENSION(:,:,:), POINTER ::  cs_scalar      !<
     2246    REAL(wp), DIMENSION(:,:,:), POINTER ::  cs_scalar_p    !<
     2247    REAL(wp), DIMENSION(:,:,:), POINTER ::  tcs_scalar_m   !<
     2248
     2249
     2250    !
     2251    !-- Tendency terms for chemical species
     2252    tend = 0.0_wp
     2253    !   
     2254    !-- Advection terms
     2255    IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2256       IF ( ws_scheme_sca )  THEN
     2257          CALL advec_s_ws( cs_scalar, 'kc' )
     2258       ELSE
     2259          CALL advec_s_pw( cs_scalar )
     2260       ENDIF
     2261    ELSE
     2262       CALL advec_s_up( cs_scalar )
     2263    ENDIF
     2264    !
     2265    !-- Diffusion terms  (the last three arguments are zero)
     2266    CALL diffusion_s( cs_scalar,                                               &
     2267         surf_def_h(0)%cssws(ilsp,:),                             &
     2268         surf_def_h(1)%cssws(ilsp,:),                             &
     2269         surf_def_h(2)%cssws(ilsp,:),                             &
     2270         surf_lsm_h%cssws(ilsp,:),                                &
     2271         surf_usm_h%cssws(ilsp,:),                                &
     2272         surf_def_v(0)%cssws(ilsp,:),                             &
     2273         surf_def_v(1)%cssws(ilsp,:),                             &
     2274         surf_def_v(2)%cssws(ilsp,:),                             &
     2275         surf_def_v(3)%cssws(ilsp,:),                             &
     2276         surf_lsm_v(0)%cssws(ilsp,:),                             &
     2277         surf_lsm_v(1)%cssws(ilsp,:),                             &
     2278         surf_lsm_v(2)%cssws(ilsp,:),                             &
     2279         surf_lsm_v(3)%cssws(ilsp,:),                             &
     2280         surf_usm_v(0)%cssws(ilsp,:),                             &
     2281         surf_usm_v(1)%cssws(ilsp,:),                             &
     2282         surf_usm_v(2)%cssws(ilsp,:),                             &
     2283         surf_usm_v(3)%cssws(ilsp,:) )
     2284    !   
     2285    !-- Prognostic equation for chemical species
     2286    DO  i = nxl, nxr
     2287       DO  j = nys, nyn   
     2288          DO  k = nzb+1, nzt
     2289             cs_scalar_p(k,j,i) =   cs_scalar(k,j,i)                           &
     2290                  + ( dt_3d  *                                 &
     2291                  (   tsc(2) * tend(k,j,i)                 &
     2292                  + tsc(3) * tcs_scalar_m(k,j,i)         &
     2293                  )                                        &
     2294                  - tsc(5) * rdf_sc(k)                       &
     2295                  * ( cs_scalar(k,j,i) - pr_init_cs(k) )    &   
     2296                  )                                          &
     2297                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )       
     2298
     2299             IF ( cs_scalar_p(k,j,i) < 0.0_wp )  cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i)
     2300          ENDDO
     2301       ENDDO
     2302    ENDDO
     2303    !
     2304    !-- Calculate tendencies for the next Runge-Kutta step
     2305    IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2306       IF ( intermediate_timestep_count == 1 )  THEN
     2307          DO  i = nxl, nxr
     2308             DO  j = nys, nyn   
     2309                DO  k = nzb+1, nzt
     2310                   tcs_scalar_m(k,j,i) = tend(k,j,i)
     2311                ENDDO
    18612312             ENDDO
    1862              lsp_usr = lsp_usr + 1
     2313          ENDDO
     2314       ELSEIF ( intermediate_timestep_count < &
     2315            intermediate_timestep_count_max )  THEN
     2316          DO  i = nxl, nxr
     2317             DO  j = nys, nyn
     2318                DO  k = nzb+1, nzt
     2319                   tcs_scalar_m(k,j,i) = - 9.5625_wp * tend(k,j,i)             &
     2320                        + 5.3125_wp * tcs_scalar_m(k,j,i)
     2321                ENDDO
     2322             ENDDO
    18632323          ENDDO
    18642324       ENDIF
    1865 
    1866     END SUBROUTINE chem_init_profiles
    1867 !
    1868 !------------------------------------------------------------------------------!
    1869 !
    1870 ! Description:
    1871 ! ------------
    1872 !> Subroutine to integrate chemical species in the given chemical mechanism
    1873 !------------------------------------------------------------------------------!
    1874 
    1875     SUBROUTINE chem_integrate_ij (i, j)
    1876 
    1877        USE cpulog,                                                              &
    1878           ONLY: cpu_log, log_point
    1879        USE statistics,                                                          &
    1880           ONLY:  weight_pres
    1881         USE control_parameters,                                                 &
    1882           ONLY:  dt_3d, intermediate_timestep_count,simulated_time
    1883 
    1884        IMPLICIT   none
    1885        INTEGER,INTENT(IN)       :: i,j
    1886 
    1887  !--   local variables
    1888        INTEGER(iwp) ::  lsp                                                     !< running index for chem spcs.
    1889        INTEGER(iwp) ::  lph                                                     !< running index for photolysis frequencies
    1890        INTEGER                  :: k,m,istatf
    1891        INTEGER,DIMENSION(20)    :: istatus
    1892        REAL(kind=wp),DIMENSION(nzb+1:nzt,nspec)                :: tmp_conc
    1893        REAL(kind=wp),DIMENSION(nzb+1:nzt)                      :: tmp_temp
    1894        REAL(kind=wp),DIMENSION(nzb+1:nzt)                      :: tmp_qvap
    1895        REAL(kind=wp),DIMENSION(nzb+1:nzt,nphot)                :: tmp_phot
    1896        REAL(kind=wp),DIMENSION(nzb+1:nzt)                      :: tmp_fact
    1897        REAL(kind=wp),DIMENSION(nzb+1:nzt)                      :: tmp_fact_i    !< conversion factor between
    1898                                                                                 !<    molecules cm^{-3} and ppm
    1899 
    1900        INTEGER,DIMENSION(nzb+1:nzt)                            :: nacc          !< Number of accepted steps
    1901        INTEGER,DIMENSION(nzb+1:nzt)                            :: nrej          !< Number of rejected steps
    1902 
    1903        REAL(wp)                         ::  conv                                !< conversion factor
    1904        REAL(wp), PARAMETER              ::  ppm2fr  = 1.0e-6_wp                 !< Conversion factor ppm to fraction
    1905        REAL(wp), PARAMETER              ::  xm_air  = 28.96_wp                  !< Mole mass of dry air
    1906        REAL(wp), PARAMETER              ::  xm_h2o  = 18.01528_wp               !< Mole mass of water vapor
    1907        REAL(wp), PARAMETER              ::  pref_i  = 1._wp / 100000.0_wp       !< inverse reference pressure (1/Pa)
    1908        REAL(wp), PARAMETER              ::  t_std   = 273.15_wp                 !< standard pressure (Pa)
    1909        REAL(wp), PARAMETER              ::  p_std   = 101325.0_wp               !< standard pressure (Pa)
    1910        REAL(wp), PARAMETER              ::  vmolcm  = 22.414e3_wp               !< Mole volume (22.414 l) in cm^{-3}
    1911        REAL(wp), PARAMETER              ::  xna     = 6.022e23_wp               !< Avogadro number (molecules/mol)
    1912 
    1913        REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local
    1914 
    1915 
    1916        REAL(kind=wp)  :: dt_chem                                             
    1917 
    1918        CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'start' )
    1919 !<     set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry
    1920        IF (chem_gasphase_on) THEN
    1921           nacc = 0
    1922           nrej = 0
    1923 
    1924        tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt)
    1925 !         ppm to molecules/cm**3
    1926 !         tmp_fact = 1.e-6_wp*6.022e23_wp/(22.414_wp*1000._wp) * 273.15_wp * hyp(nzb+1:nzt)/( 101300.0_wp * tmp_temp ) 
    1927           conv = ppm2fr * xna / vmolcm
    1928           tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std)
    1929           tmp_fact_i = 1.0_wp/tmp_fact
    1930 
    1931           IF ( humidity ) THEN
    1932              IF ( bulk_cloud_model )  THEN
    1933                 tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) * xm_air/xm_h2o * tmp_fact(:)
    1934              ELSE
    1935                 tmp_qvap(:) = q(nzb+1:nzt,j,i) * xm_air/xm_h2o * tmp_fact(:)
    1936              ENDIF
    1937           ELSE
    1938              tmp_qvap(:) = 0.01 * tmp_fact(:)                          !< Constant value for q if water vapor is not computed
     2325    ENDIF
     2326
     2327 END SUBROUTINE chem_prognostic_equations
     2328
     2329 !------------------------------------------------------------------------------!
     2330 !
     2331 ! Description:
     2332 ! ------------
     2333 !> Subroutine calculating prognostic equations for chemical species
     2334 !> (cache-optimized).
     2335 !> Routine is called separately for each chemical species over a loop from
     2336 !> prognostic_equations.
     2337 !------------------------------------------------------------------------------!
     2338 SUBROUTINE chem_prognostic_equations_ij( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs,     &
     2339      i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs,                  &
     2340      flux_l_cs, diss_l_cs )
     2341    USE pegrid         
     2342    USE advec_ws,                                                                               &
     2343         ONLY:  advec_s_ws
     2344    USE advec_s_pw_mod,                                                                         &
     2345         ONLY:  advec_s_pw
     2346    USE advec_s_up_mod,                                                                         &
     2347         ONLY:  advec_s_up
     2348    USE diffusion_s_mod,                                                                        &
     2349         ONLY:  diffusion_s
     2350    USE indices,                                                                                &
     2351         ONLY:  wall_flags_0
     2352    USE surface_mod,                                                                            &
     2353         ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
     2354
     2355
     2356    IMPLICIT NONE
     2357
     2358    REAL(wp), DIMENSION(:,:,:), POINTER   :: cs_scalar_p, cs_scalar, tcs_scalar_m
     2359
     2360    INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn, ilsp
     2361    REAL(wp), DIMENSION(nzb+1:nzt, 0:threads_per_task-1)          :: flux_s_cs   !<
     2362    REAL(wp), DIMENSION(nzb+1:nzt, 0:threads_per_task-1)          :: diss_s_cs   !<
     2363    REAL(wp), DIMENSION(nzb+1:nzt, nys:nyn, 0:threads_per_task-1) :: flux_l_cs   !<
     2364    REAL(wp), DIMENSION(nzb+1:nzt, nys:nyn, 0:threads_per_task-1) :: diss_l_cs   !<
     2365    REAL(wp), DIMENSION(0:nz+1)                                 :: pr_init_cs  !<
     2366
     2367    !
     2368    !-- local variables
     2369
     2370    INTEGER :: k
     2371    !
     2372    !--    Tendency-terms for chem spcs.
     2373    tend(:,j,i) = 0.0_wp
     2374    !   
     2375    !-- Advection terms
     2376    IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2377       IF ( ws_scheme_sca )  THEN
     2378          CALL advec_s_ws( i, j, cs_scalar, 'kc', flux_s_cs, diss_s_cs,                &
     2379               flux_l_cs, diss_l_cs, i_omp_start, tn )
     2380       ELSE
     2381          CALL advec_s_pw( i, j, cs_scalar )
     2382       ENDIF
     2383    ELSE
     2384       CALL advec_s_up( i, j, cs_scalar )
     2385    ENDIF
     2386
     2387    !
     2388
     2389    !-- Diffusion terms (the last three arguments are zero)
     2390
     2391    CALL diffusion_s( i, j, cs_scalar,                                                 &
     2392         surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:),        &
     2393         surf_def_h(2)%cssws(ilsp,:),                                     &
     2394         surf_lsm_h%cssws(ilsp,:), surf_usm_h%cssws(ilsp,:),              &
     2395         surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:),        &
     2396         surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:),        &
     2397         surf_lsm_v(0)%cssws(ilsp,:), surf_lsm_v(1)%cssws(ilsp,:),        &
     2398         surf_lsm_v(2)%cssws(ilsp,:), surf_lsm_v(3)%cssws(ilsp,:),        &
     2399         surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:),        &
     2400         surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) )
     2401
     2402    !   
     2403    !-- Prognostic equation for chem spcs
     2404    DO k = nzb+1, nzt
     2405       cs_scalar_p(k,j,i) = cs_scalar(k,j,i) + ( dt_3d  *                       &
     2406            ( tsc(2) * tend(k,j,i) +         &
     2407            tsc(3) * tcs_scalar_m(k,j,i) ) &
     2408            - tsc(5) * rdf_sc(k)             &
     2409            * ( cs_scalar(k,j,i) - pr_init_cs(k) )    &   
     2410            )                                                  &
     2411            * MERGE( 1.0_wp, 0.0_wp,                  &     
     2412            BTEST( wall_flags_0(k,j,i), 0 )   &             
     2413            )       
     2414
     2415       IF ( cs_scalar_p(k,j,i) < 0.0_wp )  cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i)    !FKS6
     2416    ENDDO
     2417
     2418    !
     2419    !-- Calculate tendencies for the next Runge-Kutta step
     2420    IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2421       IF ( intermediate_timestep_count == 1 )  THEN
     2422          DO  k = nzb+1, nzt
     2423             tcs_scalar_m(k,j,i) = tend(k,j,i)
     2424          ENDDO
     2425       ELSEIF ( intermediate_timestep_count < &
     2426            intermediate_timestep_count_max )  THEN
     2427          DO  k = nzb+1, nzt
     2428             tcs_scalar_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
     2429                  5.3125_wp * tcs_scalar_m(k,j,i)
     2430          ENDDO
     2431       ENDIF
     2432    ENDIF
     2433
     2434 END SUBROUTINE chem_prognostic_equations_ij
     2435
     2436 !
     2437 !------------------------------------------------------------------------------!
     2438 !
     2439 ! Description:
     2440 ! ------------
     2441 !> Subroutine to read restart data of chemical species
     2442 !------------------------------------------------------------------------------!
     2443
     2444 SUBROUTINE chem_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,         &
     2445      nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc,  &
     2446      nys_on_file, tmp_3d, found )   
     2447
     2448    USE control_parameters
     2449
     2450    USE indices
     2451
     2452    USE pegrid
     2453
     2454    IMPLICIT NONE
     2455
     2456    CHARACTER (LEN=20) :: spc_name_av !<   
     2457
     2458    INTEGER(iwp) ::  i, lsp          !<
     2459    INTEGER(iwp) ::  k               !<
     2460    INTEGER(iwp) ::  nxlc            !<
     2461    INTEGER(iwp) ::  nxlf            !<
     2462    INTEGER(iwp) ::  nxl_on_file     !<   
     2463    INTEGER(iwp) ::  nxrc            !<
     2464    INTEGER(iwp) ::  nxrf            !<
     2465    INTEGER(iwp) ::  nxr_on_file     !<   
     2466    INTEGER(iwp) ::  nync            !<
     2467    INTEGER(iwp) ::  nynf            !<
     2468    INTEGER(iwp) ::  nyn_on_file     !<   
     2469    INTEGER(iwp) ::  nysc            !<
     2470    INTEGER(iwp) ::  nysf            !<
     2471    INTEGER(iwp) ::  nys_on_file     !<   
     2472
     2473    LOGICAL, INTENT(OUT) :: found
     2474
     2475    REAL(wp), DIMENSION(nzb:nzt+1, nys_on_file-nbgp:nyn_on_file+nbgp, nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !< 3D array to temp store data
     2476
     2477
     2478    found = .FALSE. 
     2479
     2480
     2481    IF ( ALLOCATED(chem_species) )  THEN
     2482
     2483       DO lsp = 1, nspec
     2484
     2485          !< for time-averaged chemical conc.
     2486          spc_name_av  =  TRIM(chem_species(lsp)%name)//'_av'
     2487
     2488          IF ( restart_string(1:length) == TRIM(chem_species(lsp)%name) )    &
     2489               THEN
     2490             !< read data into tmp_3d
     2491             IF ( k == 1 )  READ ( 13 )  tmp_3d 
     2492             !< fill ..%conc in the restart run   
     2493             chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp,                  &
     2494                  nxlc-nbgp:nxrc+nbgp) =                  &
     2495                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     2496             found = .TRUE.
     2497          ELSEIF (restart_string(1:length) == spc_name_av )  THEN
     2498             IF ( k == 1 )  READ ( 13 )  tmp_3d
     2499             chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp,               &
     2500                  nxlc-nbgp:nxrc+nbgp) =               &
     2501                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     2502             found = .TRUE.
    19392503          ENDIF
    19402504
    1941           DO lsp = 1,nspec
    1942              tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:)
    1943           ENDDO
    1944 
    1945           DO lph = 1,nphot
    1946              tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)               
    1947           ENDDO
    1948 !
    1949 !--   todo: remove (kanani)
    1950 !           IF(myid == 0 .AND. chem_debug0 ) THEN
    1951 !              IF (i == 10 .and. j == 10) WRITE(0,*) 'begin chemics step ',dt_3d
    1952 !           ENDIF
    1953 
    1954 !--       Compute length of time step
    1955           IF ( call_chem_at_all_substeps )  THEN
    1956              dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
    1957           ELSE
    1958              dt_chem = dt_3d
    1959           ENDIF
    1960 
    1961           cs_time_step = dt_chem
    1962 
    1963           CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'start' )
    1964 
    1965           IF(maxval(rcntrl) > 0.0)   THEN    ! Only if rcntrl is set
    1966              IF( simulated_time <= 2*dt_3d)  THEN
    1967                  rcntrl_local = 0
    1968 !
    1969 !--   todo: remove (kanani)
    1970 !                  WRITE(9,'(a,2f10.3)') 'USE Default rcntrl in the first steps ',simulated_time,dt_3d
    1971              ELSE
    1972                  rcntrl_local = rcntrl
    1973              ENDIF
    1974           ELSE
    1975              rcntrl_local = 0
    1976           END IF
    1977 
    1978           CALL chem_gasphase_integrate (dt_chem, tmp_conc, tmp_temp, tmp_qvap, tmp_fact, tmp_phot, &
    1979                              icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus)
    1980 
    1981           CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'stop' )
    1982 
    1983           DO lsp = 1,nspec
    1984              chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:)
    1985           ENDDO
    1986 
    1987 
    1988        ENDIF
    1989           CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'stop' )
    1990 
    1991        RETURN
    1992     END SUBROUTINE chem_integrate_ij
    1993 !
    1994 !------------------------------------------------------------------------------!
    1995 !
    1996 ! Description:
    1997 ! ------------
    1998 !> Subroutine defining parin for &chemistry_parameters for chemistry model
    1999 !------------------------------------------------------------------------------!
    2000     SUBROUTINE chem_parin
    2001    
    2002        USE chem_modules
    2003        USE control_parameters
    2004      
    2005        USE kinds
    2006        USE pegrid
    2007        USE statistics
    2008            
    2009        IMPLICIT   NONE
    2010 
    2011        CHARACTER (LEN=80) ::  line                        !< dummy string that contains the current line of the parameter file
    2012        CHARACTER (LEN=3)  ::  cs_prefix
    2013 
    2014        REAL(wp), DIMENSION(nmaxfixsteps) ::   my_steps    !< List of fixed timesteps   my_step(1) = 0.0 automatic stepping
    2015        INTEGER(iwp) ::  i                                 !<
    2016        INTEGER(iwp) ::  j                                 !<
    2017        INTEGER(iwp) ::  max_pr_cs_tmp                     !<
    2018 
    2019 
    2020        NAMELIST /chemistry_parameters/  bc_cs_b,                          &
    2021                                         bc_cs_t,                          &
    2022                                         call_chem_at_all_substeps,        &
    2023                                         chem_debug0,                      &
    2024                                         chem_debug1,                      &
    2025                                         chem_debug2,                      &
    2026                                         chem_gasphase_on,                 &
    2027                                         cs_heights,                       &
    2028                                         cs_name,                          &
    2029                                         cs_profile,                       &
    2030                                         cs_surface,                       &
    2031                                         decycle_chem_lr,                  &
    2032                                         decycle_chem_ns,                  &           
    2033                                         decycle_method,                   &
    2034                                         do_depo,                          &
    2035                                         emiss_factor_main,                &
    2036                                         emiss_factor_side,                &                     
    2037                                         icntrl,                           &
    2038                                         main_street_id,                   &
    2039                                         max_street_id,                    &
    2040                                         my_steps,                         &
    2041                                         nest_chemistry,                   &
    2042                                         rcntrl,                           &
    2043                                         side_street_id,                   &
    2044                                         photolysis_scheme,                &
    2045                                         wall_csflux,                      &
    2046                                         cs_vertical_gradient,             &
    2047                                         top_csflux,                       &
    2048                                         surface_csflux,                   &
    2049                                         surface_csflux_name,              &
    2050                                         cs_surface_initial_change,        &
    2051                                         cs_vertical_gradient_level,       &
    2052 !                                       namelist parameters for emissions
    2053                                         mode_emis,                        &
    2054                                         time_fac_type,                    &
    2055                                         daytype_mdh,                      &
    2056                                         do_emis                             
    2057                              
    2058 !-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
    2059 !-- so this way we could prescribe a specific flux value for each species
    2060 !>  chemistry_parameters for initial profiles
    2061 !>  cs_names = 'O3', 'NO2', 'NO', ...   to set initial profiles)
    2062 !>  cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3)
    2063 !>  cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.)
    2064 !>  cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, .....  (chem spcs conc at height lvls chem_heights(1,:)) etc.
    2065 !>  If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)"
    2066 !>  then write these cs_surface values to chem_species(lsp)%conc_pr_init(:)
    2067 
    2068 !
    2069 !--   Read chem namelist   
    2070 
    2071        INTEGER             :: ier
    2072        CHARACTER(LEN=64)   :: text
    2073        CHARACTER(LEN=8)    :: solver_type
    2074 
    2075        icntrl    = 0
    2076        rcntrl    = 0.0_wp
    2077        my_steps  = 0.0_wp
    2078        photolysis_scheme = 'simple'
    2079        atol = 1.0_wp
    2080        rtol = 0.01_wp
     2505       ENDDO
     2506
     2507    ENDIF
     2508
     2509
     2510 END SUBROUTINE chem_rrd_local
    20812511 !
    2082  !--   Try to find chemistry package
    2083        REWIND ( 11 )
    2084        line = ' '
    2085        DO   WHILE ( INDEX( line, '&chemistry_parameters' ) == 0 )
    2086           READ ( 11, '(A)', END=20 )  line
    2087        ENDDO
    2088        BACKSPACE ( 11 )
    2089  !
    2090  !--   Read chemistry namelist
    2091        READ ( 11, chemistry_parameters, ERR = 10, END = 20 )     
    2092  !
    2093  !--   Enable chemistry model
    2094        air_chemistry = .TRUE.                   
    2095        GOTO 20
    2096 
    2097  10    BACKSPACE( 11 )
    2098        READ( 11 , '(A)') line
    2099        CALL parin_fail_message( 'chemistry_parameters', line )
    2100 
    2101  20    CONTINUE
    2102 
    2103 !
    2104 !--    check for  emission mode for chem species
    2105        IF ( (mode_emis /= 'PARAMETERIZED')  .AND. ( mode_emis /= 'DEFAULT') .AND. (mode_emis /= 'PRE-PROCESSED'  ) )   THEN
    2106             message_string = 'Incorrect mode_emiss  option select. Please check spelling'
    2107             CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
    2108        ENDIF
    2109 
    2110        t_steps = my_steps         
    2111                                    
    2112 !--    Determine the number of user-defined profiles and append them to the
    2113 !--    standard data output (data_output_pr)
    2114        max_pr_cs_tmp = 0
    2115        i = 1
    2116 
    2117        DO  WHILE ( data_output_pr(i)  /= ' '  .AND.  i <= 100 )
    2118           IF ( TRIM( data_output_pr(i)(1:3)) == 'kc_' ) THEN
    2119              max_pr_cs_tmp = max_pr_cs_tmp+1
    2120           ENDIF
    2121           i = i +1
    2122        ENDDO
    2123 
    2124        IF ( max_pr_cs_tmp > 0 ) THEN
    2125           cs_pr_namelist_found = .TRUE.
    2126           max_pr_cs = max_pr_cs_tmp
    2127        ENDIF
    2128 
    2129 !      Set Solver Type
    2130        IF(icntrl(3) == 0)   THEN
    2131           solver_type = 'rodas3'           !Default
    2132        ELSE IF(icntrl(3) == 1)   THEN
    2133           solver_type = 'ros2'
    2134        ELSE IF(icntrl(3) == 2)   THEN
    2135           solver_type = 'ros3'
    2136        ELSE IF(icntrl(3) == 3)   THEN
    2137           solver_type = 'ro4'
    2138        ELSE IF(icntrl(3) == 4)   THEN
    2139           solver_type = 'rodas3'
    2140        ELSE IF(icntrl(3) == 5)   THEN
    2141           solver_type = 'rodas4'
    2142        ELSE IF(icntrl(3) == 6)   THEN
    2143           solver_type = 'Rang3'
    2144        ELSE
    2145            message_string = 'illegal solver type'
    2146            CALL message( 'chem_parin', 'PA0506', 1, 2, 0, 6, 0 )
    2147        END IF
    2148 
    2149 !
    2150 !--   todo: remove or replace by "CALL message" mechanism (kanani)
    2151 !       write(text,*) 'gas_phase chemistry: solver_type = ',trim(solver_type)
    2152 !kk    Has to be changed to right calling sequence
    2153 !kk       CALL location_message( trim(text), .FALSE. )
    2154 !        IF(myid == 0)   THEN
    2155 !           write(9,*) ' '
    2156 !           write(9,*) 'kpp setup '
    2157 !           write(9,*) ' '
    2158 !           write(9,*) '    gas_phase chemistry: solver_type = ',trim(solver_type)
    2159 !           write(9,*) ' '
    2160 !           write(9,*) '    Hstart  = ',rcntrl(3)
    2161 !           write(9,*) '    FacMin  = ',rcntrl(4)
    2162 !           write(9,*) '    FacMax  = ',rcntrl(5)
    2163 !           write(9,*) ' '
    2164 !           IF(vl_dim > 1)    THEN
    2165 !              write(9,*) '    Vector mode                   vektor length = ',vl_dim
    2166 !           ELSE
    2167 !              write(9,*) '    Scalar mode'
    2168 !           ENDIF
    2169 !           write(9,*) ' '
    2170 !        END IF
    2171 
    2172        RETURN
    2173 
    2174     END SUBROUTINE chem_parin
    2175 
    2176 !
    2177 !------------------------------------------------------------------------------!
    2178 !
    2179 ! Description:
    2180 ! ------------
    2181 !> Subroutine calculating prognostic equations for chemical species
    2182 !> (vector-optimized).
    2183 !> Routine is called separately for each chemical species over a loop from
    2184 !> prognostic_equations.
    2185 !------------------------------------------------------------------------------!
    2186     SUBROUTINE chem_prognostic_equations ( cs_scalar_p, cs_scalar, tcs_scalar_m,  &
    2187                                            pr_init_cs, ilsp )
    2188 
    2189        USE advec_s_pw_mod,                                                        &
    2190            ONLY:  advec_s_pw
    2191        USE advec_s_up_mod,                                                        &
    2192            ONLY:  advec_s_up
    2193        USE advec_ws,                                                              &
    2194            ONLY:  advec_s_ws
    2195        USE diffusion_s_mod,                                                       &
    2196            ONLY:  diffusion_s
    2197        USE indices,                                                               &
    2198            ONLY:  nxl, nxr, nyn, nys, wall_flags_0
    2199        USE pegrid
    2200        USE surface_mod,                                                           &
    2201            ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
    2202                   surf_usm_v
    2203 
    2204        IMPLICIT NONE
    2205 
    2206        INTEGER ::  i   !< running index
    2207        INTEGER ::  j   !< running index
    2208        INTEGER ::  k   !< running index
    2209 
    2210        INTEGER(iwp),INTENT(IN) ::  ilsp          !<
    2211 
    2212        REAL(wp), DIMENSION(0:nz+1) ::  pr_init_cs   !<
    2213 
    2214        REAL(wp), DIMENSION(:,:,:), POINTER ::  cs_scalar      !<
    2215        REAL(wp), DIMENSION(:,:,:), POINTER ::  cs_scalar_p    !<
    2216        REAL(wp), DIMENSION(:,:,:), POINTER ::  tcs_scalar_m   !<
    2217 
    2218 
    2219    !
    2220    !-- Tendency terms for chemical species
    2221        tend = 0.0_wp
    2222    !   
    2223    !-- Advection terms
    2224        IF ( timestep_scheme(1:5) == 'runge' )  THEN
    2225           IF ( ws_scheme_sca )  THEN
    2226              CALL advec_s_ws( cs_scalar, 'kc' )
    2227           ELSE
    2228              CALL advec_s_pw( cs_scalar )
    2229           ENDIF
    2230        ELSE
    2231             CALL advec_s_up( cs_scalar )
    2232        ENDIF
    2233    !
    2234    !-- Diffusion terms  (the last three arguments are zero)
    2235        CALL diffusion_s( cs_scalar,                                               &
    2236                          surf_def_h(0)%cssws(ilsp,:),                             &
    2237                          surf_def_h(1)%cssws(ilsp,:),                             &
    2238                          surf_def_h(2)%cssws(ilsp,:),                             &
    2239                          surf_lsm_h%cssws(ilsp,:),                                &
    2240                          surf_usm_h%cssws(ilsp,:),                                &
    2241                          surf_def_v(0)%cssws(ilsp,:),                             &
    2242                          surf_def_v(1)%cssws(ilsp,:),                             &
    2243                          surf_def_v(2)%cssws(ilsp,:),                             &
    2244                          surf_def_v(3)%cssws(ilsp,:),                             &
    2245                          surf_lsm_v(0)%cssws(ilsp,:),                             &
    2246                          surf_lsm_v(1)%cssws(ilsp,:),                             &
    2247                          surf_lsm_v(2)%cssws(ilsp,:),                             &
    2248                          surf_lsm_v(3)%cssws(ilsp,:),                             &
    2249                          surf_usm_v(0)%cssws(ilsp,:),                             &
    2250                          surf_usm_v(1)%cssws(ilsp,:),                             &
    2251                          surf_usm_v(2)%cssws(ilsp,:),                             &
    2252                          surf_usm_v(3)%cssws(ilsp,:) )
    2253    !   
    2254    !-- Prognostic equation for chemical species
    2255        DO  i = nxl, nxr
    2256           DO  j = nys, nyn   
    2257              DO  k = nzb+1, nzt
    2258                 cs_scalar_p(k,j,i) =   cs_scalar(k,j,i)                           &
    2259                                      + ( dt_3d  *                                 &
    2260                                          (   tsc(2) * tend(k,j,i)                 &
    2261                                            + tsc(3) * tcs_scalar_m(k,j,i)         &
    2262                                          )                                        &
    2263                                        - tsc(5) * rdf_sc(k)                       &
    2264                                                 * ( cs_scalar(k,j,i) - pr_init_cs(k) )    &   
    2265                                        )                                          &
    2266                                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )       
    2267 
    2268                 IF ( cs_scalar_p(k,j,i) < 0.0_wp )  cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i)
    2269              ENDDO
    2270           ENDDO
    2271        ENDDO
    2272    !
    2273    !-- Calculate tendencies for the next Runge-Kutta step
    2274        IF ( timestep_scheme(1:5) == 'runge' )  THEN
    2275           IF ( intermediate_timestep_count == 1 )  THEN
    2276              DO  i = nxl, nxr
    2277                 DO  j = nys, nyn   
    2278                    DO  k = nzb+1, nzt
    2279                       tcs_scalar_m(k,j,i) = tend(k,j,i)
    2280                    ENDDO
    2281                 ENDDO
    2282              ENDDO
    2283           ELSEIF ( intermediate_timestep_count < &
    2284              intermediate_timestep_count_max )  THEN
    2285              DO  i = nxl, nxr
    2286                 DO  j = nys, nyn
    2287                    DO  k = nzb+1, nzt
    2288                       tcs_scalar_m(k,j,i) = - 9.5625_wp * tend(k,j,i)             &
    2289                                             + 5.3125_wp * tcs_scalar_m(k,j,i)
    2290                    ENDDO
    2291                 ENDDO
    2292              ENDDO
    2293           ENDIF
    2294        ENDIF
    2295 
    2296     END SUBROUTINE chem_prognostic_equations
    2297 
    2298 !------------------------------------------------------------------------------!
    2299 !
    2300 ! Description:
    2301 ! ------------
    2302 !> Subroutine calculating prognostic equations for chemical species
    2303 !> (cache-optimized).
    2304 !> Routine is called separately for each chemical species over a loop from
    2305 !> prognostic_equations.
    2306 !------------------------------------------------------------------------------!
    2307     SUBROUTINE chem_prognostic_equations_ij ( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs,    &
    2308                                i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs,                  &
    2309                                flux_l_cs, diss_l_cs )
    2310        USE pegrid         
    2311        USE advec_ws,        ONLY:  advec_s_ws
    2312        USE advec_s_pw_mod,  ONLY:  advec_s_pw
    2313        USE advec_s_up_mod,  ONLY:  advec_s_up
    2314        USE diffusion_s_mod, ONLY:  diffusion_s
    2315        USE indices,         ONLY:  wall_flags_0
    2316        USE surface_mod,     ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
    2317                                    surf_usm_v
    2318 
    2319 
    2320        IMPLICIT NONE
    2321 
    2322        REAL(wp), DIMENSION(:,:,:), POINTER   :: cs_scalar_p, cs_scalar, tcs_scalar_m
    2323 
    2324        INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn, ilsp
    2325        REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1)         :: flux_s_cs   !<
    2326        REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1)         :: diss_s_cs   !<
    2327        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: flux_l_cs   !<
    2328        REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: diss_l_cs   !<
    2329        REAL(wp), DIMENSION(0:nz+1)                                 :: pr_init_cs  !<
    2330 
    2331 !-- local variables
    2332 
    2333        INTEGER :: k
    2334 !
    2335 !--    Tendency-terms for chem spcs.
    2336        tend(:,j,i) = 0.0_wp
    2337 !   
    2338 !-- Advection terms
    2339        IF ( timestep_scheme(1:5) == 'runge' )  THEN
    2340           IF ( ws_scheme_sca )  THEN
    2341              CALL advec_s_ws( i, j, cs_scalar, 'kc', flux_s_cs, diss_s_cs,                &
    2342                 flux_l_cs, diss_l_cs, i_omp_start, tn )
    2343           ELSE
    2344              CALL advec_s_pw( i, j, cs_scalar )
    2345           ENDIF
    2346        ELSE
    2347             CALL advec_s_up( i, j, cs_scalar )
    2348        ENDIF
    2349 
    2350 !
    2351 
    2352 !-- Diffusion terms (the last three arguments are zero)
    2353 
    2354          CALL diffusion_s( i, j, cs_scalar,                                                 &
    2355                            surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:),        &
    2356                            surf_def_h(2)%cssws(ilsp,:),                                     &
    2357                            surf_lsm_h%cssws(ilsp,:), surf_usm_h%cssws(ilsp,:),              &
    2358                            surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:),        &
    2359                            surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:),        &
    2360                            surf_lsm_v(0)%cssws(ilsp,:), surf_lsm_v(1)%cssws(ilsp,:),        &
    2361                            surf_lsm_v(2)%cssws(ilsp,:), surf_lsm_v(3)%cssws(ilsp,:),        &
    2362                            surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:),        &
    2363                            surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) )
    2364    
    2365 !   
    2366 !-- Prognostic equation for chem spcs
    2367        DO k = nzb+1, nzt
    2368           cs_scalar_p(k,j,i) = cs_scalar(k,j,i) + ( dt_3d  *                       &
    2369                                                   ( tsc(2) * tend(k,j,i) +         &
    2370                                                     tsc(3) * tcs_scalar_m(k,j,i) ) &
    2371                                                   - tsc(5) * rdf_sc(k)             &
    2372                                                            * ( cs_scalar(k,j,i) - pr_init_cs(k) )    &   
    2373                                                   )                                                  &
    2374                                                            * MERGE( 1.0_wp, 0.0_wp,                  &     
    2375                                                                    BTEST( wall_flags_0(k,j,i), 0 )   &             
    2376                                                                     )       
    2377 
    2378           IF ( cs_scalar_p(k,j,i) < 0.0_wp )  cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i)    !FKS6
    2379        ENDDO
    2380 
    2381 !
    2382 !-- Calculate tendencies for the next Runge-Kutta step
    2383        IF ( timestep_scheme(1:5) == 'runge' )  THEN
    2384           IF ( intermediate_timestep_count == 1 )  THEN
    2385              DO  k = nzb+1, nzt
    2386                 tcs_scalar_m(k,j,i) = tend(k,j,i)
    2387              ENDDO
    2388           ELSEIF ( intermediate_timestep_count < &
    2389              intermediate_timestep_count_max )  THEN
    2390              DO  k = nzb+1, nzt
    2391                 tcs_scalar_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
    2392                    5.3125_wp * tcs_scalar_m(k,j,i)
    2393              ENDDO
    2394           ENDIF
    2395        ENDIF
    2396 
    2397     END SUBROUTINE chem_prognostic_equations_ij
    2398 
    2399 !
    2400 !------------------------------------------------------------------------------!
    2401 !
    2402 ! Description:
    2403 ! ------------
    2404 !> Subroutine to read restart data of chemical species
    2405 !------------------------------------------------------------------------------!
    2406 
    2407     SUBROUTINE chem_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,         &
    2408                                nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc,  &
    2409                                nys_on_file, tmp_3d, found )   
    2410                                        
    2411        USE control_parameters
    2412                
    2413        USE indices
    2414        
    2415        USE pegrid
    2416 
    2417        IMPLICIT NONE
    2418 
    2419        CHARACTER (LEN=20) :: spc_name_av !<   
    2420          
    2421        INTEGER(iwp) ::  i, lsp          !<
    2422        INTEGER(iwp) ::  k               !<
    2423        INTEGER(iwp) ::  nxlc            !<
    2424        INTEGER(iwp) ::  nxlf            !<
    2425        INTEGER(iwp) ::  nxl_on_file     !<   
    2426        INTEGER(iwp) ::  nxrc            !<
    2427        INTEGER(iwp) ::  nxrf            !<
    2428        INTEGER(iwp) ::  nxr_on_file     !<   
    2429        INTEGER(iwp) ::  nync            !<
    2430        INTEGER(iwp) ::  nynf            !<
    2431        INTEGER(iwp) ::  nyn_on_file     !<   
    2432        INTEGER(iwp) ::  nysc            !<
    2433        INTEGER(iwp) ::  nysf            !<
    2434        INTEGER(iwp) ::  nys_on_file     !<   
    2435        
    2436        LOGICAL, INTENT(OUT)  :: found
    2437 
    2438        REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !< 3D array to temp store data
    2439 
    2440 
    2441        found = .FALSE. 
    2442 
    2443 
    2444           IF ( ALLOCATED(chem_species) )  THEN
    2445 
    2446              DO  lsp = 1, nspec
    2447 
    2448                  !< for time-averaged chemical conc.
    2449                 spc_name_av  =  TRIM(chem_species(lsp)%name)//'_av'
    2450 
    2451                 IF (restart_string(1:length) == TRIM(chem_species(lsp)%name) )    &
    2452                 THEN
    2453                    !< read data into tmp_3d
    2454                    IF ( k == 1 )  READ ( 13 )  tmp_3d 
    2455                    !< fill ..%conc in the restart run   
    2456                    chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp,                  &
    2457                                           nxlc-nbgp:nxrc+nbgp) =                  &
    2458                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    2459                    found = .TRUE.
    2460                 ELSEIF (restart_string(1:length) == spc_name_av )  THEN
    2461                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    2462                    chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp,               &
    2463                                              nxlc-nbgp:nxrc+nbgp) =               &
    2464                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    2465                    found = .TRUE.
    2466                 ENDIF
    2467 
    2468              ENDDO
    2469 
    2470           ENDIF
    2471 
    2472 
    2473     END SUBROUTINE chem_rrd_local
    2474 !
    2475 !-------------------------------------------------------------------------------!
    2476 !> Description:
    2477 !> Calculation of horizontally averaged profiles
    2478 !> This routine is called for every statistic region (sr) defined by the user,
    2479 !> but at least for the region "total domain" (sr=0).
    2480 !> quantities.
    2481 !------------------------------------------------------------------------------!
    2482     SUBROUTINE chem_statistics( mode, sr, tn )
    2483 
    2484 
    2485        USE arrays_3d
    2486        USE indices
    2487        USE kinds
    2488        USE pegrid
    2489        USE statistics
     2512 !-------------------------------------------------------------------------------!
     2513 !> Description:
     2514 !> Calculation of horizontally averaged profiles
     2515 !> This routine is called for every statistic region (sr) defined by the user,
     2516 !> but at least for the region "total domain" (sr=0).
     2517 !> quantities.
     2518 !------------------------------------------------------------------------------!
     2519 SUBROUTINE chem_statistics( mode, sr, tn )
     2520
     2521
     2522    USE arrays_3d
     2523    USE indices
     2524    USE kinds
     2525    USE pegrid
     2526    USE statistics
    24902527
    24912528    USE user
    … …  
    24952532    CHARACTER (LEN=*) ::  mode   !<
    24962533
    2497     INTEGER(iwp) :: i    !< running index on x-axis
    2498     INTEGER(iwp) :: j    !< running index on y-axis
    2499     INTEGER(iwp) :: k    !< vertical index counter
    2500     INTEGER(iwp) :: sr   !< statistical region
    2501     INTEGER(iwp) :: tn   !< thread number
    2502     INTEGER(iwp) :: n    !<
    2503     INTEGER(iwp) :: m    !<
    2504     INTEGER(iwp) :: lpr  !< running index chem spcs
    2505 !    REAL(wp),                                                                                      &
    2506 !    DIMENSION(dots_num_palm+1:dots_max) ::                                                         &
    2507 !          ts_value_l   !<
     2534    INTEGER(iwp) ::  i    !< running index on x-axis
     2535    INTEGER(iwp) ::  j    !< running index on y-axis
     2536    INTEGER(iwp) ::  k    !< vertical index counter
     2537    INTEGER(iwp) ::  sr   !< statistical region
     2538    INTEGER(iwp) ::  tn   !< thread number
     2539    INTEGER(iwp) ::  n    !<
     2540    INTEGER(iwp) ::  m    !<
     2541    INTEGER(iwp) ::  lpr  !< running index chem spcs
     2542    !    REAL(wp),                                                                                      &
     2543    !    DIMENSION(dots_num_palm+1:dots_max) ::                                                         &
     2544    !          ts_value_l   !<
    25082545
    25092546    IF ( mode == 'profiles' )  THEN
    2510 !
    2511 !--    Sample on how to calculate horizontally averaged profiles of user-
    2512 !--    defined quantities. Each quantity is identified by the index
    2513 !--    "pr_palm+#" where "#" is an integer starting from 1. These
    2514 !--    user-profile-numbers must also be assigned to the respective strings
    2515 !--    given by data_output_pr_cs in routine user_check_data_output_pr.
    2516 !--    hom(:,:,:,:) =  dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
    2517 !                       w*pt*), dim-4 = statistical region.
     2547       !
     2548       !--    Sample on how to calculate horizontally averaged profiles of user-
     2549       !--    defined quantities. Each quantity is identified by the index
     2550       !--    "pr_palm+#" where "#" is an integer starting from 1. These
     2551       !--    user-profile-numbers must also be assigned to the respective strings
     2552       !--    given by data_output_pr_cs in routine user_check_data_output_pr.
     2553       !--    hom(:,:,:,:) =  dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
     2554       !                       w*pt*), dim-4 = statistical region.
    25182555
    25192556       !$OMP DO
    25202557       DO  i = nxl, nxr
    25212558          DO  j = nys, nyn
    2522               DO  k = nzb, nzt+1
     2559             DO  k = nzb, nzt+1
    25232560                DO lpr = 1, cs_pr_count
    25242561
    25252562                   sums_l(k,pr_palm+max_pr_user+lpr,tn) = sums_l(k,pr_palm+max_pr_user+lpr,tn) +    &
    2526                                                  chem_species(cs_pr_index(lpr))%conc(k,j,i) *       &
    2527                                                  rmask(j,i,sr)  *                                   &
    2528                                                  MERGE( 1.0_wp, 0.0_wp,                             &
    2529                                                  BTEST( wall_flags_0(k,j,i), 22 ) )
    2530                 ENDDO                                                                         
     2563                        chem_species(cs_pr_index(lpr))%conc(k,j,i) *       &
     2564                        rmask(j,i,sr)  *                                   &
     2565                        MERGE( 1.0_wp, 0.0_wp,                             &
     2566                        BTEST( wall_flags_0(k,j,i), 22 ) )
     2567                ENDDO
    25312568             ENDDO
    25322569          ENDDO
    … …  
    25362573    ENDIF
    25372574
    2538     END SUBROUTINE chem_statistics
    2539 !
    2540 !------------------------------------------------------------------------------!
    2541 !
    2542 ! Description:
    2543 ! ------------
    2544 !> Subroutine for swapping of timelevels for chemical species
    2545 !> called out from subroutine swap_timelevel
    2546 !------------------------------------------------------------------------------!
    2547 
    2548     SUBROUTINE chem_swap_timelevel (level)
    2549 
    2550           IMPLICIT   none
    2551 
    2552           INTEGER,INTENT(IN)        :: level
    2553 !--       local variables
    2554           INTEGER                   :: lsp
    2555 
    2556 
    2557           IF ( level == 0 ) THEN
    2558              DO lsp=1, nvar                                       
    2559                 chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_1(:,:,:,lsp)
    2560                 chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_2(:,:,:,lsp)
    2561              ENDDO
     2575 END SUBROUTINE chem_statistics
     2576
     2577 !
     2578 !------------------------------------------------------------------------------!
     2579 !
     2580 ! Description:
     2581 ! ------------
     2582 !> Subroutine for swapping of timelevels for chemical species
     2583 !> called out from subroutine swap_timelevel
     2584 !------------------------------------------------------------------------------!
     2585
     2586 SUBROUTINE chem_swap_timelevel( level )
     2587
     2588    IMPLICIT NONE
     2589
     2590    INTEGER(iwp), INTENT(IN) ::  level
     2591    !
     2592    !-- local variables
     2593    INTEGER(iwp)             ::  lsp
     2594
     2595
     2596    IF ( level == 0 ) THEN
     2597       DO lsp=1, nvar                                       
     2598          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_1(:,:,:,lsp)
     2599          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_2(:,:,:,lsp)
     2600       ENDDO
     2601    ELSE
     2602       DO lsp=1, nvar                                       
     2603          chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_2(:,:,:,lsp)
     2604          chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_1(:,:,:,lsp)
     2605       ENDDO
     2606    ENDIF
     2607
     2608    RETURN
     2609 END SUBROUTINE chem_swap_timelevel
     2610 !
     2611 !------------------------------------------------------------------------------!
     2612 !
     2613 ! Description:
     2614 ! ------------
     2615 !> Subroutine to write restart data for chemistry model
     2616 !------------------------------------------------------------------------------!
     2617 SUBROUTINE chem_wrd_local
     2618
     2619    IMPLICIT NONE
     2620
     2621    INTEGER(iwp) ::  lsp  !<
     2622
     2623    DO  lsp = 1, nspec
     2624       CALL wrd_write_string( TRIM( chem_species(lsp)%name ) )
     2625       WRITE ( 14 )  chem_species(lsp)%conc
     2626       CALL wrd_write_string( TRIM( chem_species(lsp)%name )//'_av' )
     2627       WRITE ( 14 )  chem_species(lsp)%conc_av
     2628    ENDDO
     2629
     2630 END SUBROUTINE chem_wrd_local
     2631
     2632
     2633
     2634 !-------------------------------------------------------------------------------!
     2635 !
     2636 ! Description:
     2637 ! ------------
     2638 !> Subroutine to calculate the deposition of gases and PMs. For now deposition
     2639 !> only takes place on lsm and usm horizontal surfaces. Default surfaces are NOT
     2640 !> considered. The deposition of particles is derived following Zhang et al.,
     2641 !> 2001, gases are deposited using the DEPAC module (van Zanten et al., 2010).
     2642 !>     
     2643 !> @TODO: Consider deposition on vertical surfaces   
     2644 !> @TODO: Consider overlaying horizontal surfaces
     2645 !> @TODO: Consider resolved vegetation   
     2646 !> @TODO: Check error messages
     2647 !-------------------------------------------------------------------------------!
     2648
     2649 SUBROUTINE chem_depo( i, j )
     2650
     2651    USE control_parameters,                                                 &   
     2652         ONLY:  dt_3d, intermediate_timestep_count, latitude
     2653
     2654    USE arrays_3d,                                                          &
     2655         ONLY:  dzw, rho_air_zw
     2656
     2657    USE date_and_time_mod,                                                  &
     2658         ONLY:  day_of_year
     2659
     2660    USE surface_mod,                                                        &
     2661         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_lsm_h,        &
     2662         surf_type, surf_usm_h
     2663
     2664    USE radiation_model_mod,                                                &
     2665         ONLY:  zenith
     2666
     2667
     2668
     2669    IMPLICIT NONE
     2670
     2671    INTEGER(iwp), INTENT(IN) ::  i
     2672    INTEGER(iwp), INTENT(IN) ::  j
     2673
     2674    INTEGER(iwp) ::  k                             !< matching k to surface m at i,j
     2675    INTEGER(iwp) ::  lsp                           !< running index for chem spcs.
     2676    INTEGER(iwp) ::  lu                            !< running index for landuse classes
     2677    INTEGER(iwp) ::  lu_palm                       !< index of PALM LSM vegetation_type at current surface element
     2678    INTEGER(iwp) ::  lup_palm                      !< index of PALM LSM pavement_type at current surface element
     2679    INTEGER(iwp) ::  luw_palm                      !< index of PALM LSM water_type at current surface element
     2680    INTEGER(iwp) ::  luu_palm                      !< index of PALM USM walls/roofs at current surface element
     2681    INTEGER(iwp) ::  lug_palm                      !< index of PALM USM green walls/roofs at current surface element
     2682    INTEGER(iwp) ::  lud_palm                      !< index of PALM USM windows at current surface element
     2683    INTEGER(iwp) ::  lu_dep                        !< matching DEPAC LU to lu_palm
     2684    INTEGER(iwp) ::  lup_dep                       !< matching DEPAC LU to lup_palm
     2685    INTEGER(iwp) ::  luw_dep                       !< matching DEPAC LU to luw_palm
     2686    INTEGER(iwp) ::  luu_dep                       !< matching DEPAC LU to luu_palm
     2687    INTEGER(iwp) ::  lug_dep                       !< matching DEPAC LU to lug_palm
     2688    INTEGER(iwp) ::  lud_dep                       !< matching DEPAC LU to lud_palm
     2689    INTEGER(iwp) ::  m                             !< index for horizontal surfaces
     2690
     2691    INTEGER(iwp) ::  pspec                         !< running index
     2692    INTEGER(iwp) ::  i_pspec                       !< index for matching depac gas component
     2693
     2694    !
     2695    !-- Vegetation                                              !<Match to DEPAC
     2696    INTEGER(iwp) ::  ind_lu_user = 0                         !< --> ERROR 
     2697    INTEGER(iwp) ::  ind_lu_b_soil = 1                       !< --> ilu_desert
     2698    INTEGER(iwp) ::  ind_lu_mixed_crops = 2                  !< --> ilu_arable
     2699    INTEGER(iwp) ::  ind_lu_s_grass = 3                      !< --> ilu_grass
     2700    INTEGER(iwp) ::  ind_lu_ev_needle_trees = 4              !< --> ilu_coniferous_forest
     2701    INTEGER(iwp) ::  ind_lu_de_needle_trees = 5              !< --> ilu_coniferous_forest
     2702    INTEGER(iwp) ::  ind_lu_ev_broad_trees = 6               !< --> ilu_tropical_forest
     2703    INTEGER(iwp) ::  ind_lu_de_broad_trees = 7               !< --> ilu_deciduous_forest
     2704    INTEGER(iwp) ::  ind_lu_t_grass = 8                      !< --> ilu_grass
     2705    INTEGER(iwp) ::  ind_lu_desert = 9                       !< --> ilu_desert
     2706    INTEGER(iwp) ::  ind_lu_tundra = 10                      !< --> ilu_other
     2707    INTEGER(iwp) ::  ind_lu_irr_crops = 11                   !< --> ilu_arable
     2708    INTEGER(iwp) ::  ind_lu_semidesert = 12                  !< --> ilu_other
     2709    INTEGER(iwp) ::  ind_lu_ice = 13                         !< --> ilu_ice
     2710    INTEGER(iwp) ::  ind_lu_marsh = 14                       !< --> ilu_other
     2711    INTEGER(iwp) ::  ind_lu_ev_shrubs = 15                   !< --> ilu_mediterrean_scrub
     2712    INTEGER(iwp) ::  ind_lu_de_shrubs = 16                   !< --> ilu_mediterrean_scrub
     2713    INTEGER(iwp) ::  ind_lu_mixed_forest = 17                !< --> ilu_coniferous_forest (ave(decid+conif))
     2714    INTEGER(iwp) ::  ind_lu_intrup_forest = 18               !< --> ilu_other (ave(other+decid))
     2715
     2716    !
     2717    !-- Water
     2718    INTEGER(iwp) ::  ind_luw_user = 0                        !< --> ERROR
     2719    INTEGER(iwp) ::  ind_luw_lake = 1                        !< --> ilu_water_inland
     2720    INTEGER(iwp) ::  ind_luw_river = 2                       !< --> ilu_water_inland
     2721    INTEGER(iwp) ::  ind_luw_ocean = 3                       !< --> ilu_water_sea
     2722    INTEGER(iwp) ::  ind_luw_pond = 4                        !< --> ilu_water_inland
     2723    INTEGER(iwp) ::  ind_luw_fountain = 5                    !< --> ilu_water_inland
     2724
     2725    !
     2726    !-- Pavement
     2727    INTEGER(iwp) ::  ind_lup_user = 0                        !< --> ERROR
     2728    INTEGER(iwp) ::  ind_lup_asph_conc = 1                   !< --> ilu_desert
     2729    INTEGER(iwp) ::  ind_lup_asph = 2                        !< --> ilu_desert
     2730    INTEGER(iwp) ::  ind_lup_conc = 3                        !< --> ilu_desert
     2731    INTEGER(iwp) ::  ind_lup_sett = 4                        !< --> ilu_desert
     2732    INTEGER(iwp) ::  ind_lup_pav_stones = 5                  !< --> ilu_desert
     2733    INTEGER(iwp) ::  ind_lup_cobblest = 6                    !< --> ilu_desert
     2734    INTEGER(iwp) ::  ind_lup_metal = 7                       !< --> ilu_desert
     2735    INTEGER(iwp) ::  ind_lup_wood = 8                        !< --> ilu_desert
     2736    INTEGER(iwp) ::  ind_lup_gravel = 9                      !< --> ilu_desert
     2737    INTEGER(iwp) ::  ind_lup_f_gravel = 10                   !< --> ilu_desert
     2738    INTEGER(iwp) ::  ind_lup_pebblest = 11                   !< --> ilu_desert
     2739    INTEGER(iwp) ::  ind_lup_woodchips = 12                  !< --> ilu_desert
     2740    INTEGER(iwp) ::  ind_lup_tartan = 13                     !< --> ilu_desert
     2741    INTEGER(iwp) ::  ind_lup_art_turf = 14                   !< --> ilu_desert
     2742    INTEGER(iwp) ::  ind_lup_clay = 15                       !< --> ilu_desert
     2743
     2744
     2745    !
     2746    !-- Particle parameters according to the respective aerosol classes (PM25, PM10)
     2747    INTEGER(iwp) ::  ind_p_size = 1     !< index for partsize in particle_pars
     2748    INTEGER(iwp) ::  ind_p_dens = 2     !< index for rhopart in particle_pars
     2749    INTEGER(iwp) ::  ind_p_slip = 3     !< index for slipcor in particle_pars
     2750
     2751    INTEGER(iwp) ::  part_type
     2752
     2753    INTEGER(iwp) ::  nwet           !< wetness indicator dor DEPAC; nwet=0 -> dry; nwet=1 -> wet; nwet=9 -> snow
     2754
     2755    REAL(wp) ::  dt_chem
     2756    REAL(wp) ::  dh               
     2757    REAL(wp) ::  inv_dh
     2758    REAL(wp) ::  dt_dh
     2759
     2760    REAL(wp) ::  dens              !< density at layer k at i,j 
     2761    REAL(wp) ::  r_aero_lu         !< aerodynamic resistance (s/m) at current surface element
     2762    REAL(wp) ::  ustar_lu          !< ustar at current surface element
     2763    REAL(wp) ::  z0h_lu            !< roughness length for heat at current surface element
     2764    REAL(wp) ::  glrad             !< rad_sw_in at current surface element
     2765    REAL(wp) ::  ppm_to_ugm3       !< conversion factor
     2766    REAL(wp) ::  relh              !< relative humidity at current surface element
     2767    REAL(wp) ::  lai               !< leaf area index at current surface element
     2768    REAL(wp) ::  sai               !< surface area index at current surface element assumed to be lai + 1
     2769
     2770    REAL(wp) ::  slinnfac       
     2771    REAL(wp) ::  visc              !< Viscosity
     2772    REAL(wp) ::  vs                !< Sedimentation velocity
     2773    REAL(wp) ::  vd_lu             !< deposition velocity (m/s)
     2774    REAL(wp) ::  Rs                !< Sedimentaion resistance (s/m)
     2775    REAL(wp) ::  Rb                !< quasi-laminar boundary layer resistance (s/m)
     2776    REAL(wp) ::  Rc_tot            !< total canopy resistance Rc (s/m)
     2777
     2778    REAL(wp) ::  catm              !< gasconc. at i, j, k in ug/m3
     2779    REAL(wp) ::  diffc             !< diffusivity
     2780
     2781
     2782    REAL(wp), DIMENSION(nspec) ::  bud_now_lu       !< budget for LSM vegetation type at current surface element
     2783    REAL(wp), DIMENSION(nspec) ::  bud_now_lup      !< budget for LSM pavement type at current surface element
     2784    REAL(wp), DIMENSION(nspec) ::  bud_now_luw      !< budget for LSM water type at current surface element
     2785    REAL(wp), DIMENSION(nspec) ::  bud_now_luu      !< budget for USM walls/roofs at current surface element
     2786    REAL(wp), DIMENSION(nspec) ::  bud_now_lug      !< budget for USM green surfaces at current surface element
     2787    REAL(wp), DIMENSION(nspec) ::  bud_now_lud      !< budget for USM windows at current surface element
     2788    REAL(wp), DIMENSION(nspec) ::  bud_now          !< overall budget at current surface element
     2789    REAL(wp), DIMENSION(nspec) ::  cc               !< concentration i,j,k
     2790    REAL(wp), DIMENSION(nspec) ::  ccomp_tot        !< total compensation point (ug/m3)
     2791    !< For now kept to zero for all species!
     2792
     2793    REAL(wp) ::  ttemp          !< temperatur at i,j,k
     2794    REAL(wp) ::  ts             !< surface temperatur in degrees celsius
     2795    REAL(wp) ::  qv_tmp         !< surface mixing ratio at current surface element
     2796
     2797    !
     2798    !-- Particle parameters (PM10 (1), PM25 (2))
     2799    !-- partsize (diameter in m), rhopart (density in kg/m3), slipcor
     2800    !-- (slip correction factor dimensionless, Seinfeld and Pandis 2006, Table 9.3)
     2801    REAL(wp), DIMENSION(1:3,1:2), PARAMETER ::  particle_pars = RESHAPE( (/ &
     2802         8.0e-6_wp, 1.14e3_wp, 1.016_wp, &  !<  1
     2803         0.7e-6_wp, 1.14e3_wp, 1.082_wp &   !<  2
     2804         /), (/ 3, 2 /) )
     2805
     2806
     2807    LOGICAL ::  match_lsm     !< flag indicating natural-type surface
     2808    LOGICAL ::  match_usm     !< flag indicating urban-type surface
     2809
     2810
     2811    !
     2812    !-- List of names of possible tracers
     2813    CHARACTER(len=*), PARAMETER ::  pspecnames(nposp) = (/ &
     2814         'NO2           ', &    !< NO2
     2815         'NO            ', &    !< NO
     2816         'O3            ', &    !< O3
     2817         'CO            ', &    !< CO
     2818         'form          ', &    !< FORM
     2819         'ald           ', &    !< ALD
     2820         'pan           ', &    !< PAN
     2821         'mgly          ', &    !< MGLY
     2822         'par           ', &    !< PAR
     2823         'ole           ', &    !< OLE
     2824         'eth           ', &    !< ETH
     2825         'tol           ', &    !< TOL
     2826         'cres          ', &    !< CRES
     2827         'xyl           ', &    !< XYL
     2828         'SO4a_f        ', &    !< SO4a_f
     2829         'SO2           ', &    !< SO2
     2830         'HNO2          ', &    !< HNO2
     2831         'CH4           ', &    !< CH4
     2832         'NH3           ', &    !< NH3
     2833         'NO3           ', &    !< NO3
     2834         'OH            ', &    !< OH
     2835         'HO2           ', &    !< HO2
     2836         'N2O5          ', &    !< N2O5
     2837         'SO4a_c        ', &    !< SO4a_c
     2838         'NH4a_f        ', &    !< NH4a_f
     2839         'NO3a_f        ', &    !< NO3a_f
     2840         'NO3a_c        ', &    !< NO3a_c
     2841         'C2O3          ', &    !< C2O3
     2842         'XO2           ', &    !< XO2
     2843         'XO2N          ', &    !< XO2N
     2844         'cro           ', &    !< CRO
     2845         'HNO3          ', &    !< HNO3
     2846         'H2O2          ', &    !< H2O2
     2847         'iso           ', &    !< ISO
     2848         'ispd          ', &    !< ISPD
     2849         'to2           ', &    !< TO2
     2850         'open          ', &    !< OPEN
     2851         'terp          ', &    !< TERP
     2852         'ec_f          ', &    !< EC_f
     2853         'ec_c          ', &    !< EC_c
     2854         'pom_f         ', &    !< POM_f
     2855         'pom_c         ', &    !< POM_c
     2856         'ppm_f         ', &    !< PPM_f
     2857         'ppm_c         ', &    !< PPM_c
     2858         'na_ff         ', &    !< Na_ff
     2859         'na_f          ', &    !< Na_f
     2860         'na_c          ', &    !< Na_c
     2861         'na_cc         ', &    !< Na_cc
     2862         'na_ccc        ', &    !< Na_ccc
     2863         'dust_ff       ', &    !< dust_ff
     2864         'dust_f        ', &    !< dust_f
     2865         'dust_c        ', &    !< dust_c
     2866         'dust_cc       ', &    !< dust_cc
     2867         'dust_ccc      ', &    !< dust_ccc
     2868         'tpm10         ', &    !< tpm10
     2869         'tpm25         ', &    !< tpm25
     2870         'tss           ', &    !< tss
     2871         'tdust         ', &    !< tdust
     2872         'tc            ', &    !< tc
     2873         'tcg           ', &    !< tcg
     2874         'tsoa          ', &    !< tsoa
     2875         'tnmvoc        ', &    !< tnmvoc
     2876         'SOxa          ', &    !< SOxa
     2877         'NOya          ', &    !< NOya
     2878         'NHxa          ', &    !< NHxa
     2879         'NO2_obs       ', &    !< NO2_obs
     2880         'tpm10_biascorr', &    !< tpm10_biascorr
     2881         'tpm25_biascorr', &    !< tpm25_biascorr
     2882         'O3_biascorr   ' /)    !< o3_biascorr
     2883
     2884
     2885    !
     2886    !-- tracer mole mass:
     2887    REAL(wp), PARAMETER ::  specmolm(nposp) = (/ &
     2888         xm_O * 2 + xm_N, &                         !< NO2
     2889         xm_O + xm_N, &                             !< NO
     2890         xm_O * 3, &                                !< O3
     2891         xm_C + xm_O, &                             !< CO
     2892         xm_H * 2 + xm_C + xm_O, &                  !< FORM
     2893         xm_H * 3 + xm_C * 2 + xm_O, &              !< ALD
     2894         xm_H * 3 + xm_C * 2 + xm_O * 5 + xm_N, &   !< PAN
     2895         xm_H * 4 + xm_C * 3 + xm_O * 2, &          !< MGLY
     2896         xm_H * 3 + xm_C, &                         !< PAR
     2897         xm_H * 3 + xm_C * 2, &                     !< OLE
     2898         xm_H * 4 + xm_C * 2, &                     !< ETH
     2899         xm_H * 8 + xm_C * 7, &                     !< TOL
     2900         xm_H * 8 + xm_C * 7 + xm_O, &              !< CRES
     2901         xm_H * 10 + xm_C * 8, &                    !< XYL
     2902         xm_S + xm_O * 4, &                         !< SO4a_f
     2903         xm_S + xm_O * 2, &                         !< SO2
     2904         xm_H + xm_O * 2 + xm_N, &                  !< HNO2
     2905         xm_H * 4 + xm_C, &                         !< CH4
     2906         xm_H * 3 + xm_N, &                         !< NH3
     2907         xm_O * 3 + xm_N, &                         !< NO3
     2908         xm_H + xm_O, &                             !< OH
     2909         xm_H + xm_O * 2, &                         !< HO2
     2910         xm_O * 5 + xm_N * 2, &                     !< N2O5
     2911         xm_S + xm_O * 4, &                         !< SO4a_c
     2912         xm_H * 4 + xm_N, &                         !< NH4a_f
     2913         xm_O * 3 + xm_N, &                         !< NO3a_f
     2914         xm_O * 3 + xm_N, &                         !< NO3a_c
     2915         xm_C * 2 + xm_O * 3, &                     !< C2O3
     2916         xm_dummy, &                                !< XO2
     2917         xm_dummy, &                                !< XO2N
     2918         xm_dummy, &                                !< CRO
     2919         xm_H + xm_O * 3 + xm_N, &                  !< HNO3
     2920         xm_H * 2 + xm_O * 2, &                     !< H2O2
     2921         xm_H * 8 + xm_C * 5, &                     !< ISO
     2922         xm_dummy, &                                !< ISPD
     2923         xm_dummy, &                                !< TO2
     2924         xm_dummy, &                                !< OPEN
     2925         xm_H * 16 + xm_C * 10, &                   !< TERP
     2926         xm_dummy, &                                !< EC_f
     2927         xm_dummy, &                                !< EC_c
     2928         xm_dummy, &                                !< POM_f
     2929         xm_dummy, &                                !< POM_c
     2930         xm_dummy, &                                !< PPM_f
     2931         xm_dummy, &                                !< PPM_c
     2932         xm_Na, &                                   !< Na_ff
     2933         xm_Na, &                                   !< Na_f
     2934         xm_Na, &                                   !< Na_c
     2935         xm_Na, &                                   !< Na_cc
     2936         xm_Na, &                                   !< Na_ccc
     2937         xm_dummy, &                                !< dust_ff
     2938         xm_dummy, &                                !< dust_f
     2939         xm_dummy, &                                !< dust_c
     2940         xm_dummy, &                                !< dust_cc
     2941         xm_dummy, &                                !< dust_ccc
     2942         xm_dummy, &                                !< tpm10
     2943         xm_dummy, &                                !< tpm25
     2944         xm_dummy, &                                !< tss
     2945         xm_dummy, &                                !< tdust
     2946         xm_dummy, &                                !< tc
     2947         xm_dummy, &                                !< tcg
     2948         xm_dummy, &                                !< tsoa
     2949         xm_dummy, &                                !< tnmvoc
     2950         xm_dummy, &                                !< SOxa
     2951         xm_dummy, &                                !< NOya
     2952         xm_dummy, &                                !< NHxa
     2953         xm_O * 2 + xm_N, &                         !< NO2_obs
     2954         xm_dummy, &                                !< tpm10_biascorr
     2955         xm_dummy, &                                !< tpm25_biascorr
     2956         xm_O * 3 /)                                !< o3_biascorr
     2957
     2958
     2959    !
     2960    !-- Initialize surface element m
     2961    m = 0
     2962    k = 0
     2963    !
     2964    !-- LSM or USM surface present at i,j:
     2965    !-- Default surfaces are NOT considered for deposition
     2966    match_lsm = surf_lsm_h%start_index(j,i) <= surf_lsm_h%end_index(j,i)
     2967    match_usm = surf_usm_h%start_index(j,i) <= surf_usm_h%end_index(j,i)
     2968
     2969
     2970    !
     2971    !--For LSM surfaces
     2972
     2973    IF ( match_lsm )  THEN
     2974
     2975       !
     2976       !-- Get surface element information at i,j:
     2977       m = surf_lsm_h%start_index(j,i)
     2978       k = surf_lsm_h%k(m)
     2979
     2980       !
     2981       !-- Get needed variables for surface element m
     2982       ustar_lu  = surf_lsm_h%us(m)
     2983       z0h_lu    = surf_lsm_h%z0h(m)
     2984       r_aero_lu = surf_lsm_h%r_a(m)
     2985       glrad     = surf_lsm_h%rad_sw_in(m)
     2986       lai = surf_lsm_h%lai(m)
     2987       sai = lai + 1
     2988
     2989       !
     2990       !-- For small grid spacing neglect R_a
     2991       IF ( dzw(k) <= 1.0 )  THEN
     2992          r_aero_lu = 0.0_wp
     2993       ENDIF
     2994
     2995       !
     2996       !--Initialize lu's
     2997       lu_palm = 0
     2998       lu_dep = 0
     2999       lup_palm = 0
     3000       lup_dep = 0
     3001       luw_palm = 0
     3002       luw_dep = 0
     3003
     3004       !
     3005       !--Initialize budgets
     3006       bud_now_lu  = 0.0_wp
     3007       bud_now_lup = 0.0_wp
     3008       bud_now_luw = 0.0_wp
     3009
     3010
     3011       !
     3012       !-- Get land use for i,j and assign to DEPAC lu
     3013       IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
     3014          lu_palm = surf_lsm_h%vegetation_type(m)
     3015          IF ( lu_palm == ind_lu_user )  THEN
     3016             message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
     3017             CALL message( 'chem_depo', 'CM0451', 1, 2, 0, 6, 0 )
     3018          ELSEIF ( lu_palm == ind_lu_b_soil )  THEN
     3019             lu_dep = 9
     3020          ELSEIF ( lu_palm == ind_lu_mixed_crops )  THEN
     3021             lu_dep = 2
     3022          ELSEIF ( lu_palm == ind_lu_s_grass )  THEN
     3023             lu_dep = 1
     3024          ELSEIF ( lu_palm == ind_lu_ev_needle_trees )  THEN
     3025             lu_dep = 4
     3026          ELSEIF ( lu_palm == ind_lu_de_needle_trees )  THEN
     3027             lu_dep = 4
     3028          ELSEIF ( lu_palm == ind_lu_ev_broad_trees )  THEN
     3029             lu_dep = 12
     3030          ELSEIF ( lu_palm == ind_lu_de_broad_trees )  THEN
     3031             lu_dep = 5
     3032          ELSEIF ( lu_palm == ind_lu_t_grass )  THEN
     3033             lu_dep = 1
     3034          ELSEIF ( lu_palm == ind_lu_desert )  THEN
     3035             lu_dep = 9
     3036          ELSEIF ( lu_palm == ind_lu_tundra )  THEN
     3037             lu_dep = 8
     3038          ELSEIF ( lu_palm == ind_lu_irr_crops )  THEN
     3039             lu_dep = 2
     3040          ELSEIF ( lu_palm == ind_lu_semidesert )  THEN
     3041             lu_dep = 8
     3042          ELSEIF ( lu_palm == ind_lu_ice )  THEN
     3043             lu_dep = 10
     3044          ELSEIF ( lu_palm == ind_lu_marsh )  THEN
     3045             lu_dep = 8
     3046          ELSEIF ( lu_palm == ind_lu_ev_shrubs )  THEN
     3047             lu_dep = 14
     3048          ELSEIF ( lu_palm == ind_lu_de_shrubs )  THEN
     3049             lu_dep = 14
     3050          ELSEIF ( lu_palm == ind_lu_mixed_forest )  THEN
     3051             lu_dep = 4
     3052          ELSEIF ( lu_palm == ind_lu_intrup_forest )  THEN
     3053             lu_dep = 8     
     3054          ENDIF
     3055       ENDIF
     3056
     3057       IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 )  THEN
     3058          lup_palm = surf_lsm_h%pavement_type(m)
     3059          IF ( lup_palm == ind_lup_user )  THEN
     3060             message_string = 'No pavement type defined. Please define pavement type to enable deposition calculation'
     3061             CALL message( 'chem_depo', 'CM0452', 1, 2, 0, 6, 0 )
     3062          ELSEIF ( lup_palm == ind_lup_asph_conc )  THEN
     3063             lup_dep = 9
     3064          ELSEIF ( lup_palm == ind_lup_asph )  THEN
     3065             lup_dep = 9
     3066          ELSEIF ( lup_palm ==  ind_lup_conc )  THEN
     3067             lup_dep = 9
     3068          ELSEIF ( lup_palm ==  ind_lup_sett )  THEN
     3069             lup_dep = 9
     3070          ELSEIF ( lup_palm == ind_lup_pav_stones )  THEN
     3071             lup_dep = 9
     3072          ELSEIF ( lup_palm == ind_lup_cobblest )  THEN
     3073             lup_dep = 9       
     3074          ELSEIF ( lup_palm == ind_lup_metal )  THEN
     3075             lup_dep = 9
     3076          ELSEIF ( lup_palm == ind_lup_wood )  THEN
     3077             lup_dep = 9   
     3078          ELSEIF ( lup_palm == ind_lup_gravel )  THEN
     3079             lup_dep = 9
     3080          ELSEIF ( lup_palm == ind_lup_f_gravel )  THEN
     3081             lup_dep = 9
     3082          ELSEIF ( lup_palm == ind_lup_pebblest )  THEN
     3083             lup_dep = 9
     3084          ELSEIF ( lup_palm == ind_lup_woodchips )  THEN
     3085             lup_dep = 9
     3086          ELSEIF ( lup_palm == ind_lup_tartan )  THEN
     3087             lup_dep = 9
     3088          ELSEIF ( lup_palm == ind_lup_art_turf )  THEN
     3089             lup_dep = 9
     3090          ELSEIF ( lup_palm == ind_lup_clay )  THEN
     3091             lup_dep = 9
     3092          ENDIF
     3093       ENDIF
     3094
     3095       IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 )  THEN
     3096          luw_palm = surf_lsm_h%water_type(m)     
     3097          IF ( luw_palm == ind_luw_user )  THEN
     3098             message_string = 'No water type defined. Please define water type to enable deposition calculation'
     3099             CALL message( 'chem_depo', 'CM0453', 1, 2, 0, 6, 0 )
     3100          ELSEIF ( luw_palm ==  ind_luw_lake )  THEN
     3101             luw_dep = 13
     3102          ELSEIF ( luw_palm == ind_luw_river )  THEN
     3103             luw_dep = 13
     3104          ELSEIF ( luw_palm == ind_luw_ocean )  THEN
     3105             luw_dep = 6
     3106          ELSEIF ( luw_palm == ind_luw_pond )  THEN
     3107             luw_dep = 13
     3108          ELSEIF ( luw_palm == ind_luw_fountain )  THEN
     3109             luw_dep = 13
     3110          ENDIF
     3111       ENDIF
     3112
     3113
     3114       !
     3115       !-- Set wetness indicator to dry or wet for lsm vegetation or pavement
     3116       IF ( surf_lsm_h%c_liq(m) > 0 )  THEN
     3117          nwet = 1
     3118       ELSE
     3119          nwet = 0
     3120       ENDIF
     3121
     3122       !
     3123       !--Compute length of time step
     3124       IF ( call_chem_at_all_substeps )  THEN
     3125          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
     3126       ELSE
     3127          dt_chem = dt_3d
     3128       ENDIF
     3129
     3130
     3131       dh = dzw(k)
     3132       inv_dh = 1.0_wp / dh
     3133       dt_dh = dt_chem/dh
     3134
     3135       !
     3136       !-- Concentration at i,j,k
     3137       DO  lsp = 1, nspec
     3138          cc(lsp) = chem_species(lsp)%conc(k,j,i)
     3139       ENDDO
     3140
     3141
     3142       !
     3143       !-- Temperature at i,j,k
     3144       ttemp = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
     3145
     3146       ts       = ttemp - 273.15  !< in degrees celcius
     3147
     3148       !
     3149       !-- Viscosity of air
     3150       visc = 1.496e-6 * ttemp**1.5 / (ttemp+120.0)
     3151
     3152       !
     3153       !-- Air density at k
     3154       dens = rho_air_zw(k)
     3155
     3156       !
     3157       !-- Calculate relative humidity from specific humidity for DEPAC
     3158       qv_tmp = MAX(q(k,j,i),0.0_wp)
     3159       relh = relativehumidity_from_specifichumidity(qv_tmp, ttemp, hyp(k) )
     3160
     3161
     3162
     3163       !
     3164       !-- Check if surface fraction (vegetation, pavement or water) > 0 and calculate vd and budget
     3165       !-- for each surface fraction. Then derive overall budget taking into account the surface fractions.
     3166
     3167       !
     3168       !-- Vegetation
     3169       IF ( surf_lsm_h%frac(ind_veg_wall,m) > 0 )  THEN
     3170
     3171
     3172          slinnfac = 1.0_wp
     3173
     3174          !
     3175          !-- Get vd
     3176          DO  lsp = 1, nvar
     3177             !
     3178             !-- Initialize
     3179             vs = 0.0_wp
     3180             vd_lu = 0.0_wp
     3181             Rs = 0.0_wp
     3182             Rb = 0.0_wp
     3183             Rc_tot = 0.0_wp
     3184             IF ( spc_names(lsp) == 'PM10' )  THEN
     3185                part_type = 1
     3186                !
     3187                !-- Sedimentation velocity
     3188                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3189                     particle_pars(ind_p_size, part_type), &
     3190                     particle_pars(ind_p_slip, part_type), &
     3191                     visc)
     3192
     3193                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3194                     vs, &
     3195                     particle_pars(ind_p_size, part_type), &
     3196                     particle_pars(ind_p_slip, part_type), &
     3197                     nwet, ttemp, dens, visc, &
     3198                     lu_dep,  &
     3199                     r_aero_lu, ustar_lu )
     3200
     3201                bud_now_lu(lsp) = - cc(lsp) * &
     3202                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     3203
     3204
     3205             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
     3206                part_type = 2
     3207                !
     3208                !-- Sedimentation velocity
     3209                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3210                     particle_pars(ind_p_size, part_type), &
     3211                     particle_pars(ind_p_slip, part_type), &
     3212                     visc)
     3213
     3214
     3215                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3216                     vs, &
     3217                     particle_pars(ind_p_size, part_type), &
     3218                     particle_pars(ind_p_slip, part_type), &
     3219                     nwet, ttemp, dens, visc, &
     3220                     lu_dep , &
     3221                     r_aero_lu, ustar_lu )
     3222
     3223                bud_now_lu(lsp) = - cc(lsp) * &
     3224                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     3225
     3226
     3227             ELSE !< GASES
     3228
     3229                !
     3230                !-- Read spc_name of current species for gas parameter
     3231
     3232                IF ( ANY( pspecnames(:) == spc_names(lsp) ) )  THEN
     3233                   i_pspec = 0
     3234                   DO  pspec = 1, nposp
     3235                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
     3236                         i_pspec = pspec
     3237                      END IF
     3238                   ENDDO
     3239
     3240                ELSE
     3241                   !
     3242                   !-- For now species not deposited
     3243                   CYCLE
     3244                ENDIF
     3245
     3246                !
     3247                !-- Factor used for conversion from ppb to ug/m3 :
     3248                !-- ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     3249                !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     3250                !   c           1e-9              xm_tracer         1e9       /       xm_air            dens
     3251                !-- thus:
     3252                !   c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     3253                !-- Use density at k:
     3254
     3255                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
     3256
     3257                !
     3258                !-- Atmospheric concentration in DEPAC is requested in ug/m3:
     3259                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
     3260                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     3261
     3262                !
     3263                !-- Diffusivity for DEPAC relevant gases
     3264                !-- Use default value
     3265                diffc            = 0.11e-4
     3266                !
     3267                !-- overwrite with known coefficients of diffusivity from Massman (1998)
     3268                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
     3269                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     3270                IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
     3271                IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
     3272                IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
     3273                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
     3274                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
     3275
     3276
     3277                !
     3278                !-- Get quasi-laminar boundary layer resistance Rb:
     3279                CALL get_rb_cell( (lu_dep == ilu_water_sea) .OR. (lu_dep == ilu_water_inland), &
     3280                     z0h_lu, ustar_lu, diffc, &
     3281                     Rb )
     3282
     3283                !
     3284                !-- Get Rc_tot
     3285                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     3286                     relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
     3287                     r_aero_lu , Rb )
     3288
     3289
     3290                !
     3291                !-- Calculate budget
     3292                IF ( Rc_tot <= 0.0 )  THEN
     3293
     3294                   bud_now_lu(lsp) = 0.0_wp
     3295
     3296                ELSE
     3297
     3298                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
     3299
     3300                   bud_now_lu(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     3301                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     3302                ENDIF
     3303
     3304             ENDIF
     3305          ENDDO
     3306       ENDIF
     3307
     3308
     3309       !
     3310       !-- Pavement
     3311       IF ( surf_lsm_h%frac(ind_pav_green,m) > 0 )  THEN
     3312
     3313
     3314          !
     3315          !-- No vegetation on pavements:
     3316          lai = 0.0_wp
     3317          sai = 0.0_wp
     3318
     3319          slinnfac = 1.0_wp
     3320
     3321          !
     3322          !-- Get vd
     3323          DO  lsp = 1, nvar
     3324             !
     3325             !-- Initialize
     3326             vs = 0.0_wp
     3327             vd_lu = 0.0_wp
     3328             Rs = 0.0_wp
     3329             Rb = 0.0_wp
     3330             Rc_tot = 0.0_wp
     3331             IF ( spc_names(lsp) == 'PM10' )  THEN
     3332                part_type = 1
     3333                !
     3334                !-- Sedimentation velocity:
     3335                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3336                     particle_pars(ind_p_size, part_type), &
     3337                     particle_pars(ind_p_slip, part_type), &
     3338                     visc)
     3339
     3340                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3341                     vs, &
     3342                     particle_pars(ind_p_size, part_type), &
     3343                     particle_pars(ind_p_slip, part_type), &
     3344                     nwet, ttemp, dens, visc, &
     3345                     lup_dep,  &
     3346                     r_aero_lu, ustar_lu )
     3347
     3348                bud_now_lup(lsp) = - cc(lsp) * &
     3349                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     3350
     3351
     3352             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
     3353                part_type = 2
     3354                !
     3355                !-- Sedimentation velocity:
     3356                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3357                     particle_pars(ind_p_size, part_type), &
     3358                     particle_pars(ind_p_slip, part_type), &
     3359                     visc)
     3360
     3361
     3362                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3363                     vs, &
     3364                     particle_pars(ind_p_size, part_type), &
     3365                     particle_pars(ind_p_slip, part_type), &
     3366                     nwet, ttemp, dens, visc, &
     3367                     lup_dep, &
     3368                     r_aero_lu, ustar_lu )
     3369
     3370                bud_now_lup(lsp) = - cc(lsp) * &
     3371                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     3372
     3373
     3374             ELSE  !<GASES
     3375
     3376                !
     3377                !-- Read spc_name of current species for gas parameter
     3378
     3379                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
     3380                   i_pspec = 0
     3381                   DO  pspec = 1, nposp
     3382                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
     3383                         i_pspec = pspec
     3384                      END IF
     3385                   ENDDO
     3386
     3387                ELSE
     3388                   !
     3389                   !-- For now species not deposited
     3390                   CYCLE
     3391                ENDIF
     3392
     3393                !-- Factor used for conversion from ppb to ug/m3 :
     3394                !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     3395                !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     3396                !   c           1e-9               xm_tracer         1e9       /       xm_air            dens
     3397                !-- thus:
     3398                !   c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     3399                !-- Use density at lowest layer:
     3400
     3401                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
     3402
     3403                !
     3404                !-- Atmospheric concentration in DEPAC is requested in ug/m3:
     3405                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
     3406                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     3407
     3408                !
     3409                !-- Diffusivity for DEPAC relevant gases
     3410                !-- Use default value
     3411                diffc            = 0.11e-4
     3412                !
     3413                !-- overwrite with known coefficients of diffusivity from Massman (1998)
     3414                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
     3415                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     3416                IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
     3417                IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
     3418                IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
     3419                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
     3420                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
     3421
     3422
     3423                !
     3424                !-- Get quasi-laminar boundary layer resistance Rb:
     3425                CALL get_rb_cell( (lup_dep == ilu_water_sea) .OR. (lup_dep == ilu_water_inland), &
     3426                     z0h_lu, ustar_lu, diffc, &
     3427                     Rb )
     3428
     3429
     3430                !
     3431                !-- Get Rc_tot
     3432                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     3433                     relh, lai, sai, nwet, lup_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
     3434                     r_aero_lu , Rb )
     3435
     3436
     3437                !
     3438                !-- Calculate budget
     3439                IF ( Rc_tot <= 0.0 )  THEN
     3440
     3441                   bud_now_lup(lsp) = 0.0_wp
     3442
     3443                ELSE
     3444
     3445                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
     3446
     3447                   bud_now_lup(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     3448                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     3449                ENDIF
     3450
     3451
     3452             ENDIF
     3453          ENDDO
     3454       ENDIF
     3455
     3456
     3457       !
     3458       !-- Water
     3459       IF ( surf_lsm_h%frac(ind_wat_win,m) > 0 )  THEN
     3460
     3461
     3462          !
     3463          !-- No vegetation on water:
     3464          lai = 0.0_wp
     3465          sai = 0.0_wp
     3466
     3467          slinnfac = 1.0_wp
     3468
     3469          !
     3470          !-- Get vd
     3471          DO  lsp = 1, nvar
     3472             !
     3473             !-- Initialize
     3474             vs = 0.0_wp
     3475             vd_lu = 0.0_wp
     3476             Rs = 0.0_wp
     3477             Rb = 0.0_wp
     3478             Rc_tot = 0.0_wp
     3479             IF ( spc_names(lsp) == 'PM10' )  THEN
     3480                part_type = 1
     3481                !
     3482                !-- Sedimentation velocity:
     3483                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3484                     particle_pars(ind_p_size, part_type), &
     3485                     particle_pars(ind_p_slip, part_type), &
     3486                     visc)
     3487
     3488                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3489                     vs, &
     3490                     particle_pars(ind_p_size, part_type), &
     3491                     particle_pars(ind_p_slip, part_type), &
     3492                     nwet, ttemp, dens, visc, &
     3493                     luw_dep, &
     3494                     r_aero_lu, ustar_lu )
     3495
     3496                bud_now_luw(lsp) = - cc(lsp) * &
     3497                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     3498
     3499
     3500             ELSEIF ( spc_names(lsp) == 'PM25' )  THEN
     3501                part_type = 2
     3502                !
     3503                !-- Sedimentation velocity:
     3504                vs = slinnfac * sedimentation_velocity( particle_pars(ind_p_dens, part_type), &
     3505                     particle_pars(ind_p_size, part_type), &
     3506                     particle_pars(ind_p_slip, part_type), &
     3507                     visc)
     3508
     3509
     3510                CALL drydepo_aero_zhang_vd( vd_lu, Rs, &
     3511                     vs, &
     3512                     particle_pars(ind_p_size, part_type), &
     3513                     particle_pars(ind_p_slip, part_type), &
     3514                     nwet, ttemp, dens, visc, &
     3515                     luw_dep, &
     3516                     r_aero_lu, ustar_lu )
     3517
     3518                bud_now_luw(lsp) = - cc(lsp) * &
     3519                     (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     3520
     3521
     3522             ELSE  !<GASES
     3523
     3524                !
     3525                !-- Read spc_name of current species for gas PARAMETER
     3526
     3527                IF ( ANY(pspecnames(:) == spc_names(lsp) ) )  THEN
     3528                   i_pspec = 0
     3529                   DO  pspec = 1, nposp
     3530                      IF ( pspecnames(pspec) == spc_names(lsp) )  THEN
     3531                         i_pspec = pspec
     3532                      END IF
     3533                   ENDDO
     3534
     3535                ELSE
     3536                   !
     3537                   !-- For now species not deposited
     3538                   CYCLE
     3539                ENDIF
     3540
     3541                !-- Factor used for conversion from ppb to ug/m3 :
     3542                !   ppb (mole tr)/(mole air)/ppb (kg tr)/(mole tr) (ug tr)/(kg tr) &
     3543                !   (mole air)/(kg air) (kg air)/(m3 air) (kg air(ug/m3)/ppb/(kg/mole) = / (kg/mole)
     3544                !   c           1e-9               xm_tracer         1e9       /       xm_air            dens
     3545                !-- thus:
     3546                !   c_in_ppb * xm_tracer * [ dens / xm_air ] = c_in_ugm3
     3547                !-- Use density at lowest layer:
     3548
     3549                ppm_to_ugm3 =  (dens/xm_air) * 0.001_wp  !< (mole air)/m3
     3550
     3551                !
     3552                !-- Atmospheric concentration in DEPAC is requested in ug/m3:
     3553                !   ug/m3              ppm          (ug/m3)/ppm/(kg/mole)     kg/mole
     3554                catm = cc(lsp)         * ppm_to_ugm3 *   specmolm(i_pspec)  ! in ug/m3
     3555
     3556                !
     3557                !-- Diffusivity for DEPAC relevant gases
     3558                !-- Use default value
     3559                diffc            = 0.11e-4
     3560                !
     3561                !-- overwrite with known coefficients of diffusivity from Massman (1998)
     3562                IF ( spc_names(lsp) == 'NO2' ) diffc = 0.136e-4
     3563                IF ( spc_names(lsp) == 'NO'  ) diffc = 0.199e-4
     3564                IF ( spc_names(lsp) == 'O3'  ) diffc = 0.144e-4
     3565                IF ( spc_names(lsp) == 'CO'  ) diffc = 0.176e-4
     3566                IF ( spc_names(lsp) == 'SO2' ) diffc = 0.112e-4
     3567                IF ( spc_names(lsp) == 'CH4' ) diffc = 0.191e-4
     3568                IF ( spc_names(lsp) == 'NH3' ) diffc = 0.191e-4
     3569
     3570
     3571                !
     3572                !-- Get quasi-laminar boundary layer resistance Rb:
     3573                CALL get_rb_cell( (luw_dep == ilu_water_sea) .OR. (luw_dep == ilu_water_inland), &
     3574                     z0h_lu, ustar_lu, diffc, &
     3575                     Rb )
     3576
     3577                !
     3578                !-- Get Rc_tot
     3579                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
     3580                     relh, lai, sai, nwet, luw_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
     3581                     r_aero_lu , Rb )
     3582
     3583
     3584                ! Calculate budget
     3585                IF ( Rc_tot <= 0.0 )  THEN
     3586
     3587                   bud_now_luw(lsp) = 0.0_wp
     3588
     3589                ELSE
     3590
     3591                   vd_lu = 1.0_wp / (r_aero_lu + Rb + Rc_tot )
     3592
     3593                   bud_now_luw(lsp) = - (cc(lsp) - ccomp_tot(lsp)) * &
     3594                        (1.0_wp - exp(-vd_lu * dt_dh )) * dh
     3595                ENDIF
     3596
     3597             ENDIF
     3598          ENDDO
     3599       ENDIF
     3600
     3601
     3602       bud_now = 0.0_wp
     3603       !
     3604       !-- Calculate overall budget for surface m and adapt concentration
     3605       DO  lsp = 1, nspec
     3606
     3607
     3608          bud_now(lsp) = surf_lsm_h%frac(ind_veg_wall,m) * bud_now_lu(lsp) + &
     3609               surf_lsm_h%frac(ind_pav_green,m) * bud_now_lup(lsp) + &
     3610               surf_lsm_h%frac(ind_wat_win,m) * bud_now_luw(lsp)
     3611
     3612          !
     3613          !-- Compute new concentration:
     3614          cc(lsp) = cc(lsp) + bud_now(lsp) * inv_dh
     3615
     3616          chem_species(lsp)%conc(k,j,i) = MAX(0.0_wp, cc(lsp))
     3617
     3618       ENDDO
     3619
     3620    ENDIF
     3621
     3622
     3623
     3624
     3625    !
     3626    !-- For USM surfaces   
     3627
     3628    IF ( match_usm )  THEN
     3629
     3630       !
     3631       !-- Get surface element information at i,j:
     3632       m = surf_usm_h%start_index(j,i)
     3633       k = surf_usm_h%k(m)
     3634
     3635       !
     3636       !-- Get needed variables for surface element m
     3637       ustar_lu  = surf_usm_h%us(m)
     3638       z0h_lu    = surf_usm_h%z0h(m)
     3639       r_aero_lu = surf_usm_h%r_a(m)
     3640       glrad     = surf_usm_h%rad_sw_in(m)
     3641       lai = surf_usm_h%lai(m)
     3642       sai = lai + 1
     3643
     3644       !
     3645       !-- For small grid spacing neglect R_a
     3646       IF ( dzw(k) <= 1.0 )  THEN
     3647          r_aero_lu = 0.0_wp
     3648       ENDIF
     3649
     3650       !
     3651       !-- Initialize lu's
     3652       luu_palm = 0
     3653       luu_dep = 0
     3654       lug_palm = 0
     3655       lug_dep = 0
     3656       lud_palm = 0
     3657       lud_dep = 0
     3658
     3659       !
     3660       !-- Initialize budgets
     3661       bud_now_luu  = 0.0_wp
     3662       bud_now_lug = 0.0_wp
     3663       bud_now_lud = 0.0_wp
     3664
     3665
     3666       !
     3667       !-- Get land use for i,j and assign to DEPAC lu
     3668       IF ( surf_usm_h%frac(ind_pav_green,m) > 0 )  THEN
     3669          !
     3670          !-- For green urban surfaces (e.g. green roofs
     3671          !-- assume LU short grass
     3672          lug_palm = ind_lu_s_grass
     3673          IF ( lug_palm == ind_lu_user )  THEN
     3674             message_string = 'No vegetation type defined. Please define vegetation type to enable deposition calculation'
     3675             CALL message( 'chem_depo', 'CM0454', 1, 2, 0, 6, 0 )
     3676          ELSEIF ( lug_palm == ind_lu_b_soil )  THEN
     3677             lug_dep = 9
     3678          ELSEIF ( lug_palm == ind_lu_mixed_crops )  THEN
     3679             lug_dep = 2
     3680          ELSEIF ( lug_palm == ind_lu_s_grass )  THEN
     3681             lug_dep = 1
     3682          ELSEIF ( lug_palm == ind_lu_ev_needle_trees )  THEN
     3683             lug_dep = 4
     3684          ELSEIF ( lug_palm == ind_lu_de_needle_trees )  THEN
     3685             lug_dep = 4
     3686          ELSEIF ( lug_palm == ind_lu_ev_broad_trees )  THEN
     3687             lug_dep = 12
     3688          ELSEIF ( lug_palm == ind_lu_de_broad_trees )  THEN
     3689             lug_dep = 5
     3690          ELSEIF ( lug_palm == ind_lu_t_grass )  THEN
     3691             lug_dep = 1
     3692          ELSEIF ( lug_palm == ind_lu_desert )  THEN
     3693             lug_dep = 9
     3694          ELSEIF ( lug_palm == ind_lu_tundra  )  THEN
     3695             lug_dep = 8
     3696          ELSEIF ( lug_palm == ind