source: palm/trunk/SOURCE/chemistry_model_mod.f90 @ 3435

Last change on this file since 3435 was 3435, checked in by gronemeier, 4 years ago

new: terrain-following masked output; bugfixes: increase vertical dimension of gamma_w_green_sat by 1, add checks for masked output for chemistry_model and radiation_model, reordered calls to xxx_define_netcdf_grid in masked output part

  • Property svn:keywords set to Id
File size: 98.4 KB
Line 
1!> @file chemistry_model_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2017-2018 Leibniz Universitaet Hannover
18! Copyright 2017-2018 Karlsruhe Institute of Technology
19! Copyright 2017-2018 Freie Universitaet Berlin
20!------------------------------------------------------------------------------!
21!
22! Current revisions:
23! -----------------
24!
25!
26! Former revisions:
27! -----------------
28! $Id: chemistry_model_mod.f90 3435 2018-10-26 18:25:44Z gronemeier $
29! Add terrain-following masked output
30!
31! 3373 2018-10-18 15:25:56Z kanani
32! Remove MPI_Abort, replace by message
33!
34! 3318 2018-10-08 11:43:01Z sward
35! Fixed faulty syntax of message string
36!
37! 3298 2018-10-02 12:21:11Z kanani
38! Add remarks (kanani)
39! Merge with trunk, replaced cloud_physics by bulk_cloud_model         28.09.2018 forkel
40! Subroutines header and chem_check_parameters added                   25.09.2018 basit
41! Removed chem_emission routine now declared in chem_emissions.f90     30.07.2018 ERUSSO
42! Introduced emissions namelist parameters                             30.07.2018 ERUSSO
43!
44! Timestep steering added in subroutine chem_integrate_ij and
45! output of chosen solver in chem_parin added                          30.07.2018 ketelsen
46!
47! chem_check_data_output_pr: added unit for PM compounds               20.07.2018 forkel
48! replaced : by nzb+1:nzt for pt,q,ql (found by kk)                    18.07.2018 forkel
49! debugged restart run for chem species               06.07.2018 basit
50! reorganized subroutines in alphabetical order.      27.06.2018 basit
51! subroutine chem_parin updated for profile output    27.06.2018 basit
52! Added humidity arrays to USE section and tmp_qvap in chem_integrate  26.6.2018 forkel
53! Merged chemistry with with trunk (nzb_do and nzt_do in 3d output)    26.6.2018 forkel
54!
55! reorganized subroutines in alphabetical order.      basit 22.06.2018
56! subroutine chem_parin updated for profile output    basit 22.06.2018
57! subroutine chem_statistics added                 
58! subroutine chem_check_data_output_pr add            21.06.2018 basit
59! subroutine chem_data_output_mask added              20.05.2018 basit
60! subroutine chem_data_output_2d added                20.05.2018 basit
61! subroutine chem_statistics added                    04.06.2018 basit
62! subroutine chem_emissions: Set cssws to zero before setting values 20.03.2018 forkel
63! subroutine chem_emissions: Introduced different conversion factors
64! for PM and gaseous compounds                                    15.03.2018 forkel
65! subroutine chem_emissions updated to take variable number of chem_spcs and
66! emission factors.                                               13.03.2018 basit
67! chem_boundary_conds_decycle improved.                           05.03.2018 basit
68! chem_boundary_conds_decycle subroutine added                    21.02.2018 basit
69! chem_init_profiles subroutines re-activated after correction    21.02.2018 basit
70!
71!
72! 3293 2018-09-28 12:45:20Z forkel
73! Modularization of all bulk cloud physics code components
74!
75! 3248 2018-09-14 09:42:06Z sward
76! Minor formating changes
77!
78! 3246 2018-09-13 15:14:50Z sward
79! Added error handling for input namelist via parin_fail_message
80!
81! 3241 2018-09-12 15:02:00Z raasch
82! +nest_chemistry
83!
84! 3209 2018-08-27 16:58:37Z suehring
85! Rename flags indicating outflow boundary conditions
86!
87! 3182 2018-07-27 13:36:03Z suehring
88! Revise output of surface quantities in case of overhanging structures
89!
90! 3045 2018-05-28 07:55:41Z Giersch
91! error messages revised
92!
93! 3014 2018-05-09 08:42:38Z maronga
94! Bugfix: nzb_do and nzt_do were not used for 3d data output
95!
96! 3004 2018-04-27 12:33:25Z Giersch
97! Comment concerning averaged data output added
98!
99! 2932 2018-03-26 09:39:22Z maronga
100! renamed chemistry_par to chemistry_parameters
101!
102! 2894 2018-03-15 09:17:58Z Giersch
103! Calculations of the index range of the subdomain on file which overlaps with
104! the current subdomain are already done in read_restart_data_mod,
105! chem_last_actions was renamed to chem_wrd_local, chem_read_restart_data was
106! renamed to chem_rrd_local, chem_write_var_list was renamed to
107! chem_wrd_global, chem_read_var_list was renamed to chem_rrd_global,
108! chem_skip_var_list has been removed, variable named found has been
109! introduced for checking if restart data was found, reading of restart strings
110! has been moved completely to read_restart_data_mod, chem_rrd_local is already
111! inside the overlap loop programmed in read_restart_data_mod, todo list has
112! bees extended, redundant characters in chem_wrd_local have been removed,
113! the marker *** end chemistry *** is not necessary anymore, strings and their
114! respective lengths are written out and read now in case of restart runs to
115! get rid of prescribed character lengths
116!
117! 2815 2018-02-19 11:29:57Z suehring
118! Bugfix in restart mechanism,
119! rename chem_tendency to chem_prognostic_equations,
120! implement vector-optimized version of chem_prognostic_equations,
121! some clean up (incl. todo list)
122!
123! 2773 2018-01-30 14:12:54Z suehring
124! Declare variables required for nesting as public
125!
126! 2772 2018-01-29 13:10:35Z suehring
127! Bugfix in string handling
128!
129! 2768 2018-01-24 15:38:29Z kanani
130! Shorten lines to maximum length of 132 characters
131!
132! 2766 2018-01-22 17:17:47Z kanani
133! Removed preprocessor directive __chem
134!
135! 2756 2018-01-16 18:11:14Z suehring
136! Fill values in 3D output introduced.
137!
138! 2718 2018-01-02 08:49:38Z maronga
139! Initial revision
140!
141!
142!
143!
144! Authors:
145! --------
146! @author Renate Forkel
147! @author Farah Kanani-Suehring
148! @author Klaus Ketelsen
149! @author Basit Khan
150!
151!
152!------------------------------------------------------------------------------!
153! Description:
154! ------------
155!> Chemistry model for PALM-4U
156!> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not
157!>       allowed to use the chemistry model in a precursor run and additionally
158!>       not using it in a main run
159!> @todo Update/clean-up todo list! (FK)
160!> @todo Set proper fill values (/= 0) for chem output arrays! (FK)
161!> @todo Add routine chem_check_parameters, add checks for inconsistent or
162!>       unallowed parameter settings.
163!>       CALL of chem_check_parameters from check_parameters. (FK)
164!> @todo Make routine chem_header available, CALL from header.f90
165!>       (see e.g. how it is done in routine lsm_header in
166!>        land_surface_model_mod.f90). chem_header should include all setup
167!>        info about chemistry parameter settings. (FK)
168!> @todo Implement turbulent inflow of chem spcs in inflow_turbulence.
169!> @todo Separate boundary conditions for each chem spcs to be implemented
170!> @todo Currently only total concentration are calculated. Resolved, parameterized
171!>       and chemistry fluxes although partially and some completely coded but
172!>       are not operational/activated in this version. bK.
173!> @todo slight differences in passive scalar and chem spcs when chem reactions
174!>       turned off. Need to be fixed. bK
175!> @todo test nesting for chem spcs, was implemented by suehring (kanani)
176!> @todo chemistry error messages
177!> @todo Format this module according to PALM coding standard (see e.g. module
178!>       template under http://palm.muk.uni-hannover.de/mosaik/downloads/8 or
179!>       D3_coding_standard.pdf under https://palm.muk.uni-hannover.de/trac/downloads/16)
180!
181!------------------------------------------------------------------------------!
182
183 MODULE chemistry_model_mod
184
185    USE kinds,              ONLY: wp, iwp
186    USE indices,            ONLY: nz, nzb,nzt,nysg,nyng,nxlg,nxrg,                                 &
187                                 nys,nyn,nx,nxl,nxr,ny,wall_flags_0
188    USE pegrid,             ONLY: myid, threads_per_task
189
190   USE bulk_cloud_model_mod,                                               &
191        ONLY:  bulk_cloud_model
192
193    USE control_parameters, ONLY: bc_lr_cyc, bc_ns_cyc, dt_3d, humidity,                           &
194                                 initializing_actions, message_string,                             &
195                                 omega, tsc, intermediate_timestep_count,                          &
196                                 intermediate_timestep_count_max,                                  &
197                                 timestep_scheme, use_prescribed_profile_data, ws_scheme_sca 
198   USE arrays_3d,          ONLY: exner, hyp, pt, q, ql, rdf_sc, tend, zu
199    USE chem_gasphase_mod,  ONLY: nspec, spc_names, nkppctrl, nmaxfixsteps,                        &
200                                 t_steps, chem_gasphase_integrate, vl_dim,                         &
201                                 nvar, nreact,  atol, rtol, nphot, phot_names
202    USE cpulog,             ONLY: cpu_log, log_point
203
204    USE chem_modules
205 
206    USE statistics
207
208    IMPLICIT   NONE
209    PRIVATE
210    SAVE
211
212    LOGICAL ::  nest_chemistry = .TRUE. !< flag for nesting mode of chemical species, independent on parent or not
213!
214!- Define chemical variables
215    TYPE   species_def
216       CHARACTER(LEN=8)                                   :: name
217       CHARACTER(LEN=16)                                  :: unit
218       REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: conc
219       REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: conc_av
220       REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: conc_p
221       REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: tconc_m
222       REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:)           :: cssws_av           
223       REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:)           :: flux_s_cs
224       REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:)           :: diss_s_cs
225       REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:)         :: flux_l_cs
226       REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:)         :: diss_l_cs
227       REAL(kind=wp),ALLOCATABLE,DIMENSION(:)             :: conc_pr_init
228    END TYPE species_def
229
230
231    TYPE   photols_def                                                           
232       CHARACTER(LEN=8)                                   :: name
233       CHARACTER(LEN=16)                                  :: unit
234       REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: freq
235    END TYPE photols_def
236
237
238    PUBLIC  species_def
239    PUBLIC  photols_def
240
241
242    TYPE(species_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC     :: chem_species
243    TYPE(photols_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC     :: phot_frequen 
244
245    REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: spec_conc_1
246    REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: spec_conc_2
247    REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: spec_conc_3
248    REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: spec_conc_av       
249    REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: freq_1
250
251    INTEGER,DIMENSION(nkppctrl)                           :: icntrl                            ! Fine tuning kpp
252    REAL(kind=wp),DIMENSION(nkppctrl)                     :: rcntrl                            ! Fine tuning kpp
253    LOGICAL :: decycle_chem_lr    = .FALSE.
254    LOGICAL :: decycle_chem_ns    = .FALSE.
255    CHARACTER (LEN=20), DIMENSION(4) ::  decycle_method = &
256                  (/'dirichlet','dirichlet','dirichlet','dirichlet'/)
257                                 !< Decycling method at horizontal boundaries,
258                                 !< 1=left, 2=right, 3=south, 4=north
259                                 !< dirichlet = initial size distribution and
260                                 !< chemical composition set for the ghost and
261                                 !< first three layers
262                                 !< neumann = zero gradient
263
264    REAL(kind=wp), PUBLIC :: cs_time_step
265    CHARACTER(10), PUBLIC :: photolysis_scheme
266                                         ! 'constant',
267                                         ! 'simple' (Simple parameterisation from MCM, Saunders et al., 2003, Atmos. Chem. Phys., 3, 161-180
268                                         ! 'fastj'  (Wild et al., 2000, J. Atmos. Chem., 37, 245-282) STILL NOT IMPLEMENTED
269
270    PUBLIC nest_chemistry
271    PUBLIC nspec
272    PUBLIC nreact
273    PUBLIC nvar     
274    PUBLIC spc_names
275    PUBLIC spec_conc_2 
276
277!-  Interface section
278    INTERFACE chem_3d_data_averaging
279       MODULE PROCEDURE chem_3d_data_averaging 
280    END INTERFACE chem_3d_data_averaging
281
282    INTERFACE chem_boundary_conds
283       MODULE PROCEDURE chem_boundary_conds
284       MODULE PROCEDURE chem_boundary_conds_decycle
285    END INTERFACE chem_boundary_conds
286
287    INTERFACE chem_check_data_output
288       MODULE PROCEDURE chem_check_data_output
289    END INTERFACE chem_check_data_output
290
291    INTERFACE chem_data_output_2d
292       MODULE PROCEDURE chem_data_output_2d
293    END INTERFACE chem_data_output_2d
294
295    INTERFACE chem_data_output_3d
296       MODULE PROCEDURE chem_data_output_3d
297    END INTERFACE chem_data_output_3d
298
299    INTERFACE chem_data_output_mask
300       MODULE PROCEDURE chem_data_output_mask
301    END INTERFACE chem_data_output_mask
302
303    INTERFACE chem_check_data_output_pr
304       MODULE PROCEDURE chem_check_data_output_pr
305    END INTERFACE chem_check_data_output_pr
306
307    INTERFACE chem_check_parameters
308       MODULE PROCEDURE chem_check_parameters
309    END INTERFACE chem_check_parameters
310
311    INTERFACE chem_define_netcdf_grid
312       MODULE PROCEDURE chem_define_netcdf_grid
313    END INTERFACE chem_define_netcdf_grid
314
315    INTERFACE chem_header
316       MODULE PROCEDURE chem_header
317    END INTERFACE chem_header
318
319    INTERFACE chem_init
320       MODULE PROCEDURE chem_init
321    END INTERFACE chem_init
322
323    INTERFACE chem_init_profiles
324       MODULE PROCEDURE chem_init_profiles
325    END INTERFACE chem_init_profiles
326
327    INTERFACE chem_integrate
328       MODULE PROCEDURE chem_integrate_ij
329    END INTERFACE chem_integrate
330
331    INTERFACE chem_parin
332       MODULE PROCEDURE chem_parin
333    END INTERFACE chem_parin
334
335    INTERFACE chem_prognostic_equations
336       MODULE PROCEDURE chem_prognostic_equations
337       MODULE PROCEDURE chem_prognostic_equations_ij
338    END INTERFACE chem_prognostic_equations
339
340    INTERFACE chem_rrd_local
341       MODULE PROCEDURE chem_rrd_local
342    END INTERFACE chem_rrd_local
343
344    INTERFACE chem_statistics
345       MODULE PROCEDURE chem_statistics
346    END INTERFACE chem_statistics
347
348    INTERFACE chem_swap_timelevel
349       MODULE PROCEDURE chem_swap_timelevel
350    END INTERFACE chem_swap_timelevel
351
352    INTERFACE chem_wrd_local
353       MODULE PROCEDURE chem_wrd_local 
354    END INTERFACE chem_wrd_local
355
356
357    PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_check_data_output, &
358           chem_check_data_output_pr, chem_check_parameters,                    &
359           chem_data_output_2d, chem_data_output_3d, chem_data_output_mask,     &
360           chem_define_netcdf_grid, chem_header, chem_init,                     &
361           chem_init_profiles, chem_integrate, chem_parin,                      &
362           chem_prognostic_equations, chem_rrd_local,                           &
363           chem_statistics, chem_swap_timelevel, chem_wrd_local 
364
365
366
367 CONTAINS
368
369!
370!------------------------------------------------------------------------------!
371!
372! Description:
373! ------------
374!> Subroutine for averaging 3D data of chemical species. Due to the fact that
375!> the averaged chem arrays are allocated in chem_init, no if-query concerning
376!> the allocation is required (in any mode). Attention: If you just specify an
377!> averaged output quantity in the _p3dr file during restarts the first output
378!> includes the time between the beginning of the restart run and the first
379!> output time (not necessarily the whole averaging_interval you have
380!> specified in your _p3d/_p3dr file )
381!------------------------------------------------------------------------------!
382
383    SUBROUTINE chem_3d_data_averaging ( mode, variable )
384
385       USE control_parameters
386
387       USE indices
388
389       USE kinds
390
391       USE surface_mod,                                                         &
392          ONLY: surf_def_h, surf_lsm_h, surf_usm_h   
393 
394       IMPLICIT NONE
395 
396       CHARACTER (LEN=*) ::  mode    !<
397       CHARACTER (LEN=*) :: variable !<
398 
399       LOGICAL      ::  match_def !< flag indicating natural-type surface
400       LOGICAL      ::  match_lsm !< flag indicating natural-type surface
401       LOGICAL      ::  match_usm !< flag indicating urban-type surface
402
403       INTEGER(iwp) ::  i                  !< grid index x direction
404       INTEGER(iwp) ::  j                  !< grid index y direction
405       INTEGER(iwp) ::  k                  !< grid index z direction
406       INTEGER(iwp) ::  m                  !< running index surface type
407       INTEGER(iwp) :: lsp                 !< running index for chem spcs
408
409
410       IF ( mode == 'allocate' )  THEN
411          DO lsp = 1, nspec
412             IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
413                chem_species(lsp)%conc_av = 0.0_wp
414               
415             ENDIF
416          ENDDO
417
418       ELSEIF ( mode == 'sum' )  THEN
419
420
421          DO lsp = 1, nspec
422             IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
423                DO  i = nxlg, nxrg
424                   DO  j = nysg, nyng
425                      DO  k = nzb, nzt+1
426                          chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) + &
427                                                             chem_species(lsp)%conc(k,j,i)
428                      ENDDO
429                   ENDDO
430                ENDDO
431             ELSEIF ( TRIM(variable(4:)) == TRIM('cssws*') )        THEN
432                DO  i = nxl, nxr
433                   DO  j = nys, nyn
434                      match_def = surf_def_h(0)%start_index(j,i) <=               &
435                                  surf_def_h(0)%end_index(j,i)
436                      match_lsm = surf_lsm_h%start_index(j,i) <=                  &
437                                  surf_lsm_h%end_index(j,i)
438                      match_usm = surf_usm_h%start_index(j,i) <=                  &
439                                  surf_usm_h%end_index(j,i)
440
441                      IF ( match_def )  THEN
442                         m = surf_def_h(0)%end_index(j,i)
443                         chem_species(lsp)%cssws_av(j,i) =                        &
444                                                chem_species(lsp)%cssws_av(j,i) + &
445                                                surf_def_h(0)%cssws(lsp,m)
446                      ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
447                         m = surf_lsm_h%end_index(j,i)
448                         chem_species(lsp)%cssws_av(j,i) =                        &
449                                                chem_species(lsp)%cssws_av(j,i) + &
450                                                surf_lsm_h%cssws(lsp,m)
451                      ELSEIF ( match_usm )  THEN
452                         m = surf_usm_h%end_index(j,i)
453                         chem_species(lsp)%cssws_av(j,i) =                        &
454                                                chem_species(lsp)%cssws_av(j,i) + &
455                                                surf_usm_h%cssws(lsp,m)
456                      ENDIF
457                   ENDDO
458                ENDDO
459             ENDIF
460          ENDDO
461 
462       ELSEIF ( mode == 'average' )  THEN
463          DO lsp = 1, nspec
464             IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
465                DO  i = nxlg, nxrg
466                   DO  j = nysg, nyng
467                      DO  k = nzb, nzt+1
468                          chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
469                      ENDDO
470                   ENDDO
471                ENDDO
472
473             ELSEIF (TRIM(variable(4:)) == TRIM('cssws*') )            THEN
474                DO i = nxlg, nxrg
475                   DO  j = nysg, nyng
476                        chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp )
477                   ENDDO
478                ENDDO
479                   CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp )                 
480             ENDIF
481          ENDDO
482         
483       ENDIF     
484
485    END SUBROUTINE chem_3d_data_averaging
486
487!   
488!------------------------------------------------------------------------------!
489!
490! Description:
491! ------------
492!> Subroutine to initialize and set all boundary conditions for chemical species
493!------------------------------------------------------------------------------!
494
495    SUBROUTINE chem_boundary_conds( mode )                                           
496                                                                                     
497       USE control_parameters,                                                    & 
498          ONLY:  air_chemistry, bc_radiation_l, bc_radiation_n, bc_radiation_r,  &
499                 bc_radiation_s               
500       USE indices,                                                               & 
501          ONLY:  nxl, nxr,  nxlg, nxrg, nyng, nysg, nzt                             
502                                                                                     
503
504       USE arrays_3d,                                                             &     
505          ONLY:  dzu                                               
506       USE surface_mod,                                                           &
507          ONLY:  bc_h                                                           
508
509       CHARACTER (len=*), INTENT(IN) :: mode
510       INTEGER(iwp) ::  i                                                            !< grid index x direction.
511       INTEGER(iwp) ::  j                                                            !< grid index y direction.
512       INTEGER(iwp) ::  k                                                            !< grid index z direction.
513       INTEGER(iwp) ::  kb                                                           !< variable to set respective boundary value, depends on facing.
514       INTEGER(iwp) ::  l                                                            !< running index boundary type, for up- and downward-facing walls.
515       INTEGER(iwp) ::  m                                                            !< running index surface elements.
516       INTEGER(iwp) ::  lsp                                                          !< running index for chem spcs.
517       INTEGER(iwp) ::  lph                                                          !< running index for photolysis frequencies
518
519
520       SELECT CASE ( TRIM( mode ) )       
521          CASE ( 'init' )
522
523             IF ( bc_cs_b == 'dirichlet' )    THEN
524                ibc_cs_b = 0 
525             ELSEIF ( bc_cs_b == 'neumann' )  THEN
526                ibc_cs_b = 1 
527             ELSE
528                message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"' 
529                CALL message( 'chem_boundary_conds', 'CM0429', 1, 2, 0, 6, 0 )     !<
530             ENDIF                                                                 
531!
532!--          Set Integer flags and check for possible erroneous settings for top
533!--          boundary condition.
534             IF ( bc_cs_t == 'dirichlet' )             THEN
535                ibc_cs_t = 0 
536             ELSEIF ( bc_cs_t == 'neumann' )           THEN
537                ibc_cs_t = 1
538             ELSEIF ( bc_cs_t == 'initial_gradient' )  THEN
539                ibc_cs_t = 2
540             ELSEIF ( bc_cs_t == 'nested' )            THEN         
541                ibc_cs_t = 3
542             ELSE
543                message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"'     
544                CALL message( 'check_parameters', 'CM0430', 1, 2, 0, 6, 0 )
545             ENDIF
546
547         
548          CASE ( 'set_bc_bottomtop' )                   
549!--          Bottom boundary condtions for chemical species     
550             DO lsp = 1, nspec                                                     
551                IF ( ibc_cs_b == 0 ) THEN                   
552                   DO l = 0, 1 
553!--                Set index kb: For upward-facing surfaces (l=0), kb=-1, i.e.
554!--                the chem_species(nspec)%conc_p value at the topography top (k-1)
555!--                is set; for downward-facing surfaces (l=1), kb=1, i.e. the
556!--                value at the topography bottom (k+1) is set.
557
558                      kb = MERGE( -1, 1, l == 0 )
559                      !$OMP PARALLEL DO PRIVATE( i, j, k )
560                      DO m = 1, bc_h(l)%ns
561                         i = bc_h(l)%i(m)           
562                         j = bc_h(l)%j(m)
563                         k = bc_h(l)%k(m)
564                         chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc(k+kb,j,i) 
565                      ENDDO                                       
566                   ENDDO                                       
567             
568                ELSEIF ( ibc_cs_b == 1 ) THEN
569!--             in boundary_conds there is som extra loop over m here for passive tracer
570                   DO l = 0, 1
571                      kb = MERGE( -1, 1, l == 0 )
572                      !$OMP PARALLEL DO PRIVATE( i, j, k )                                           
573                      DO m = 1, bc_h(l)%ns
574                         i = bc_h(l)%i(m)           
575                         j = bc_h(l)%j(m)
576                         k = bc_h(l)%k(m)
577                         chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc_p(k,j,i)
578
579                      ENDDO
580                   ENDDO
581                ENDIF
582          ENDDO    ! end lsp loop 
583
584!--    Top boundary conditions for chemical species - Should this not be done for all species?
585             IF ( ibc_cs_t == 0 )  THEN                   
586                DO lsp = 1, nspec
587                   chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)       
588                ENDDO
589             ELSEIF ( ibc_cs_t == 1 )  THEN
590                DO lsp = 1, nspec
591                   chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
592                ENDDO
593             ELSEIF ( ibc_cs_t == 2 )  THEN
594                DO lsp = 1, nspec
595                   chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1)
596                ENDDO
597             ENDIF
598!
599          CASE ( 'set_bc_lateral' )                       
600!--             Lateral boundary conditions for chem species at inflow boundary
601!--             are automatically set when chem_species concentration is
602!--             initialized. The initially set value at the inflow boundary is not
603!--             touched during time integration, hence, this boundary value remains
604!--             at a constant value, which is the concentration that flows into the
605!--             domain.                                                           
606!--             Lateral boundary conditions for chem species at outflow boundary
607
608             IF ( bc_radiation_s )  THEN
609                DO lsp = 1, nspec
610                   chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
611                ENDDO
612             ELSEIF ( bc_radiation_n )  THEN
613                DO lsp = 1, nspec
614                   chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
615                ENDDO
616             ELSEIF ( bc_radiation_l )  THEN
617                DO lsp = 1, nspec
618                   chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
619                ENDDO
620             ELSEIF ( bc_radiation_r )  THEN
621                DO lsp = 1, nspec
622                   chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
623                ENDDO
624             ENDIF
625           
626       END SELECT
627
628    END SUBROUTINE chem_boundary_conds
629
630!
631!------------------------------------------------------------------------------!
632! Description:
633! ------------
634!> Boundary conditions for prognostic variables in chemistry: decycling in the
635!> x-direction
636!-----------------------------------------------------------------------------
637    SUBROUTINE chem_boundary_conds_decycle (cs_3d, cs_pr_init )
638       USE pegrid,                                                             &
639                 ONLY: myid
640
641       IMPLICIT NONE
642       INTEGER(iwp) ::  boundary !<
643       INTEGER(iwp) ::  ee !<
644       INTEGER(iwp) ::  copied !<
645       INTEGER(iwp) ::  i !<
646       INTEGER(iwp) ::  j !<
647       INTEGER(iwp) ::  k !<
648       INTEGER(iwp) ::  ss !<
649       REAL(wp), DIMENSION(nzb:nzt+1) ::  cs_pr_init
650       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  cs_3d
651       REAL(wp) ::  flag !< flag to mask topography grid points
652
653       flag = 0.0_wp
654
655       
656!--    Left and right boundaries
657       IF ( decycle_chem_lr  .AND.  bc_lr_cyc )  THEN
658
659          DO  boundary = 1, 2
660
661             IF ( decycle_method(boundary) == 'dirichlet' )  THEN
662!   
663!--             Initial profile is copied to ghost and first three layers         
664                ss = 1
665                ee = 0
666                IF ( boundary == 1  .AND.  nxl == 0 )  THEN
667                   ss = nxlg
668                   ee = nxl+2
669                ELSEIF ( boundary == 2  .AND.  nxr == nx )  THEN
670                   ss = nxr-2
671                   ee = nxrg
672                ENDIF
673
674                DO  i = ss, ee
675                   DO  j = nysg, nyng
676                      DO  k = nzb+1, nzt
677                         flag = MERGE( 1.0_wp, 0.0_wp,                            &
678                                       BTEST( wall_flags_0(k,j,i), 0 ) )
679                         cs_3d(k,j,i) = cs_pr_init(k) * flag
680                      ENDDO
681                   ENDDO
682                ENDDO
683
684           ELSEIF ( decycle_method(boundary) == 'neumann' )  THEN
685!
686!--             The value at the boundary is copied to the ghost layers to simulate
687!--             an outlet with zero gradient
688                ss = 1
689                ee = 0
690                IF ( boundary == 1  .AND.  nxl == 0 )  THEN
691                   ss = nxlg
692                   ee = nxl-1
693                   copied = nxl
694                ELSEIF ( boundary == 2  .AND.  nxr == nx )  THEN
695                   ss = nxr+1
696                   ee = nxrg
697                   copied = nxr
698                ENDIF
699
700                 DO  i = ss, ee
701                   DO  j = nysg, nyng
702                      DO  k = nzb+1, nzt
703                         flag = MERGE( 1.0_wp, 0.0_wp,                            &
704                                       BTEST( wall_flags_0(k,j,i), 0 ) )
705                        cs_3d(k,j,i) = cs_3d(k,j,copied) * flag
706                      ENDDO
707                   ENDDO
708                ENDDO
709
710             ELSE
711                WRITE(message_string,*)                                           &
712                                    'unknown decycling method: decycle_method (', &
713                        boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
714                CALL message( 'chem_boundary_conds_decycle', 'CM0431',           &
715                              1, 2, 0, 6, 0 )
716             ENDIF
717          ENDDO
718       ENDIF
719
720       
721!--    South and north boundaries
722       IF ( decycle_chem_ns  .AND.  bc_ns_cyc )  THEN
723
724          DO  boundary = 3, 4
725
726             IF ( decycle_method(boundary) == 'dirichlet' )  THEN
727!   
728!--             Initial profile is copied to ghost and first three layers         
729                ss = 1
730                ee = 0
731                IF ( boundary == 3  .AND.  nys == 0 )  THEN
732                   ss = nysg
733                   ee = nys+2
734                ELSEIF ( boundary == 4  .AND.  nyn == ny )  THEN
735                   ss = nyn-2
736                   ee = nyng
737                ENDIF
738
739                DO  i = nxlg, nxrg
740                   DO  j = ss, ee
741                      DO  k = nzb+1, nzt
742                         flag = MERGE( 1.0_wp, 0.0_wp,                            &
743                                       BTEST( wall_flags_0(k,j,i), 0 ) )
744                         cs_3d(k,j,i) = cs_pr_init(k) * flag
745                      ENDDO
746                   ENDDO
747                ENDDO
748       
749       
750           ELSEIF ( decycle_method(boundary) == 'neumann' )  THEN
751!
752!--             The value at the boundary is copied to the ghost layers to simulate
753!--             an outlet with zero gradient
754                ss = 1
755                ee = 0
756                IF ( boundary == 3  .AND.  nys == 0 )  THEN
757                   ss = nysg
758                   ee = nys-1
759                   copied = nys
760                ELSEIF ( boundary == 4  .AND.  nyn == ny )  THEN
761                   ss = nyn+1
762                   ee = nyng
763                   copied = nyn
764                ENDIF
765
766                 DO  i = nxlg, nxrg
767                   DO  j = ss, ee
768                      DO  k = nzb+1, nzt
769                         flag = MERGE( 1.0_wp, 0.0_wp,                            &
770                                       BTEST( wall_flags_0(k,j,i), 0 ) )
771                         cs_3d(k,j,i) = cs_3d(k,copied,i) * flag
772                      ENDDO
773                   ENDDO
774                ENDDO
775
776             ELSE
777                WRITE(message_string,*)                                           &
778                                    'unknown decycling method: decycle_method (', &
779                        boundary, ') ="' // TRIM( decycle_method(boundary) ) // '"'
780                CALL message( 'chem_boundary_conds_decycle', 'CM0432',           &
781                              1, 2, 0, 6, 0 )
782             ENDIF
783          ENDDO
784       ENDIF
785    END SUBROUTINE chem_boundary_conds_decycle
786   !
787!------------------------------------------------------------------------------!
788!
789! Description:
790! ------------
791!> Subroutine for checking data output for chemical species
792!------------------------------------------------------------------------------!
793
794    SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
795
796
797       USE control_parameters,                                                 &
798          ONLY: data_output, message_string
799
800       IMPLICIT NONE
801
802       CHARACTER (LEN=*) ::  unit     !<
803       CHARACTER (LEN=*) ::  var      !<
804
805       INTEGER(iwp) :: i, lsp
806       INTEGER(iwp) :: ilen
807       INTEGER(iwp) :: k
808
809       CHARACTER(len=16)    ::  spec_name
810
811       unit = 'illegal'
812
813       spec_name = TRIM(var(4:))             !< var 1:3 is 'kc_'.
814
815          DO lsp=1,nspec
816             IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
817                unit = 'ppm'
818             ENDIF
819!            It is possible  to plant PM10 and PM25 into the gasphase chemistry code
820!            as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
821!            set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
822             IF (spec_name(1:2) == 'PM')   THEN 
823            unit = 'kg m-3'
824             ENDIF
825          ENDDO
826
827          DO lsp=1,nphot                                                     
828             IF (TRIM(spec_name) == TRIM(phot_frequen(lsp)%name))   THEN
829                unit = 'sec-1'
830             ENDIF
831          ENDDO
832
833
834       RETURN
835    END SUBROUTINE chem_check_data_output
836!
837!------------------------------------------------------------------------------!
838!
839! Description:
840! ------------
841!> Subroutine for checking data output of profiles for chemistry model
842!------------------------------------------------------------------------------!
843
844    SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
845
846       USE arrays_3d
847       USE control_parameters,                                                    &
848           ONLY: data_output_pr, message_string, air_chemistry
849       USE indices
850       USE profil_parameter
851       USE statistics
852
853
854       IMPLICIT NONE
855
856       CHARACTER (LEN=*) ::  unit     !<
857       CHARACTER (LEN=*) ::  variable !<
858       CHARACTER (LEN=*) ::  dopr_unit
859       CHARACTER(len=16) ::  spec_name
860   
861       INTEGER(iwp) ::  var_count, lsp  !<
862       
863
864       spec_name = TRIM(variable(4:))   
865
866          IF (  .NOT.  air_chemistry )  THEN
867             message_string = 'data_output_pr = ' //                        &
868             TRIM( data_output_pr(var_count) ) // ' is not imp' // &
869                         'lemented for air_chemistry = .FALSE.'
870             CALL message( 'chem_check_parameters', 'CM0433', 1, 2, 0, 6, 0 )             
871
872          ELSE
873             DO lsp = 1, nspec
874                IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN
875                   cs_pr_count = cs_pr_count+1
876                   cs_pr_index(cs_pr_count) = lsp
877                   dopr_index(var_count) = pr_palm+cs_pr_count 
878                   dopr_unit  = 'ppm'
879                   IF (spec_name(1:2) == 'PM')   THEN
880                        dopr_unit = 'kg m-3'
881                   ENDIF
882                      hom(:,2, dopr_index(var_count),:) = SPREAD( zu, 2, statistic_regions+1 )
883                      unit = dopr_unit 
884                ENDIF
885             ENDDO
886          ENDIF
887
888    END SUBROUTINE chem_check_data_output_pr
889
890!
891!------------------------------------------------------------------------------!
892! Description:
893! ------------
894!> Check parameters routine for chemistry_model_mod
895!------------------------------------------------------------------------------!
896    SUBROUTINE chem_check_parameters
897
898       IMPLICIT NONE
899
900       LOGICAL                        ::  found
901       INTEGER (iwp)                  ::  lsp_usr      !< running index for user defined chem spcs
902       INTEGER (iwp)                  ::  lsp          !< running index for chem spcs.
903
904
905!!--   check for chemical reactions status
906       IF ( chem_gasphase_on )  THEN
907          message_string = 'Chemical reactions: ON'
908          CALL message( 'chem_check_parameters', 'CM0421', 0, 0, 0, 6, 0 )
909       ELSEIF ( .not. (chem_gasphase_on) ) THEN
910          message_string = 'Chemical reactions: OFF'
911          CALL message( 'chem_check_parameters', 'CM0422', 0, 0, 0, 6, 0 )
912       ENDIF
913
914!--    check for chemistry time-step
915       IF ( call_chem_at_all_substeps )  THEN
916          message_string = 'Chemistry is calculated at all meteorology time-step'
917          CALL message( 'chem_check_parameters', 'CM0423', 0, 0, 0, 6, 0 )
918       ELSEIF ( .not. (call_chem_at_all_substeps) ) THEN
919          message_string = 'Sub-time-steps are skipped for chemistry time-steps'
920          CALL message( 'chem_check_parameters', 'CM0424', 0, 0, 0, 6, 0 )
921       ENDIF
922
923!--    check for photolysis scheme
924       IF ( (photolysis_scheme /= 'simple') .AND. (photolysis_scheme /= 'constant')  )  THEN
925          message_string = 'Incorrect photolysis scheme selected, please check spelling'
926          CALL message( 'chem_check_parameters', 'CM0425', 1, 2, 0, 6, 0 )
927       ENDIF
928
929!--    check for  decycling of chem species
930       IF ( (.not. any(decycle_method == 'neumann') ) .AND. (.not. any(decycle_method == 'dirichlet') ) )   THEN
931          message_string = 'Incorrect boundary conditions. Only neumann or '   &
932                   // 'dirichlet &available for decycling chemical species '
933          CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 )
934       ENDIF
935
936!--    check for initial chem species input
937       lsp_usr = 1
938       lsp     = 1
939       DO WHILE ( cs_name (lsp_usr) /= 'novalue')
940          found = .FALSE.
941          DO lsp = 1, nvar
942             IF ( TRIM(cs_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN
943                found = .TRUE.
944                EXIT
945             ENDIF
946          ENDDO
947          IF ( .not. found ) THEN
948             message_string = 'Incorrect input for initial surface vaue: ' // TRIM(cs_name(lsp_usr))
949             CALL message( 'chem_check_parameters', 'CM0427', 0, 1, 0, 6, 0 )
950          ENDIF
951          lsp_usr = lsp_usr + 1
952       ENDDO
953
954!--    check for surface  emission flux chem species
955
956       lsp_usr = 1
957       lsp     = 1
958       DO WHILE ( surface_csflux_name (lsp_usr) /= 'novalue')
959          found = .FALSE.
960          DO lsp = 1, nvar
961             IF ( TRIM(surface_csflux_name (lsp_usr)) == TRIM(chem_species(lsp)%name) ) THEN
962                found = .TRUE.
963                EXIT
964             ENDIF
965          ENDDO
966          IF ( .not. found ) THEN
967             message_string = 'Incorrect input of chemical species for surface emission fluxes: '  &
968                               // TRIM(surface_csflux_name(lsp_usr))
969             CALL message( 'chem_check_parameters', 'CM0428', 0, 1, 0, 6, 0 )
970          ENDIF
971          lsp_usr = lsp_usr + 1
972       ENDDO
973
974    END SUBROUTINE chem_check_parameters
975
976!
977!------------------------------------------------------------------------------!
978!
979! Description:
980! ------------
981!> Subroutine defining 2D output variables for chemical species
982!> @todo: Remove "mode" from argument list, not used.
983!------------------------------------------------------------------------------!
984   
985    SUBROUTINE chem_data_output_2d( av, variable, found, grid, mode, local_pf,   &
986                                     two_d, nzb_do, nzt_do, fill_value )
987   
988       USE indices
989
990       USE kinds
991
992       USE pegrid,             ONLY: myid, threads_per_task
993
994       IMPLICIT NONE
995
996       CHARACTER (LEN=*) ::  grid       !<
997       CHARACTER (LEN=*) ::  mode       !<
998       CHARACTER (LEN=*) ::  variable   !<
999       INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or time-averaged data
1000       INTEGER(iwp) ::  nzb_do          !< lower limit of the domain (usually nzb)
1001       INTEGER(iwp) ::  nzt_do          !< upper limit of the domain (usually nzt+1)
1002       LOGICAL      ::  found           !<
1003       LOGICAL      ::  two_d           !< flag parameter that indicates 2D variables (horizontal cross sections)
1004       REAL(wp)     ::  fill_value 
1005       REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf !<
1006
1007       !-- local variables.
1008       CHARACTER(len=16)    ::  spec_name
1009       INTEGER(iwp) ::  lsp
1010       INTEGER(iwp) ::  i               !< grid index along x-direction
1011       INTEGER(iwp) ::  j               !< grid index along y-direction
1012       INTEGER(iwp) ::  k               !< grid index along z-direction
1013       INTEGER(iwp) ::  m               !< running index surface elements
1014       INTEGER(iwp) ::  char_len        !< length of a character string
1015       found = .TRUE.
1016       char_len  = LEN_TRIM(variable)
1017
1018       spec_name = TRIM( variable(4:char_len-3) )
1019
1020          DO lsp=1,nspec
1021             IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name)    .AND.                             &
1022                   ( (variable(char_len-2:) == '_xy')               .OR.                              &
1023                     (variable(char_len-2:) == '_xz')               .OR.                              &
1024                     (variable(char_len-2:) == '_yz') ) )               THEN             
1025!
1026!--   todo: remove or replace by "CALL message" mechanism (kanani)
1027!                    IF(myid == 0)  WRITE(6,*) 'Output of species ' // TRIM(variable)  //               &
1028!                                                             TRIM(chem_species(lsp)%name)       
1029                IF (av == 0) THEN
1030                   DO  i = nxl, nxr
1031                      DO  j = nys, nyn
1032                         DO  k = nzb_do, nzt_do
1033                              local_pf(i,j,k) = MERGE(                                                &
1034                                                 chem_species(lsp)%conc(k,j,i),                       &
1035                                                 REAL( fill_value, KIND = wp ),                       &
1036                                                 BTEST( wall_flags_0(k,j,i), 0 ) )
1037
1038
1039                         ENDDO
1040                      ENDDO
1041                   ENDDO
1042             
1043                ELSE
1044                   DO  i = nxl, nxr
1045                      DO  j = nys, nyn
1046                         DO  k = nzb_do, nzt_do
1047                              local_pf(i,j,k) = MERGE(                                                &
1048                                                 chem_species(lsp)%conc(k,j,i),                       &
1049                                                 REAL( fill_value, KIND = wp ),                       &
1050                                                 BTEST( wall_flags_0(k,j,i), 0 ) )
1051                         ENDDO
1052                      ENDDO
1053                   ENDDO
1054                ENDIF
1055                 grid = 'zu'           
1056             ENDIF
1057          ENDDO
1058
1059          RETURN
1060     
1061    END SUBROUTINE chem_data_output_2d     
1062
1063!
1064!------------------------------------------------------------------------------!
1065!
1066! Description:
1067! ------------
1068!> Subroutine defining 3D output variables for chemical species
1069!------------------------------------------------------------------------------!
1070
1071    SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
1072
1073
1074       USE indices
1075
1076       USE kinds
1077
1078
1079       IMPLICIT NONE
1080
1081       CHARACTER (LEN=*)    ::  variable     !<
1082       INTEGER(iwp)         ::  av           !<
1083       INTEGER(iwp) ::  nzb_do               !< lower limit of the data output (usually 0)
1084       INTEGER(iwp) ::  nzt_do               !< vertical upper limit of the data output (usually nz_do3d)
1085
1086       LOGICAL      ::  found                !<
1087
1088       REAL(wp)             ::  fill_value !<
1089       REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf 
1090
1091
1092       !-- local variables
1093
1094       INTEGER              ::  i, j, k, lsp
1095       CHARACTER(len=16)    ::  spec_name
1096
1097
1098       found = .FALSE.
1099
1100       spec_name = TRIM(variable(4:))
1101
1102       DO lsp=1,nspec
1103          IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
1104!
1105!--   todo: remove or replace by "CALL message" mechanism (kanani)
1106!              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM(variable)  //  &
1107!                                                           TRIM(chem_species(lsp)%name)       
1108             
1109             IF (av == 0) THEN
1110                DO  i = nxl, nxr
1111                   DO  j = nys, nyn
1112                      DO  k = nzb_do, nzt_do
1113                          local_pf(i,j,k) = MERGE(                             &
1114                                              chem_species(lsp)%conc(k,j,i),   &
1115                                              REAL( fill_value, KIND = wp ),   &
1116                                              BTEST( wall_flags_0(k,j,i), 0 ) )
1117                      ENDDO
1118                   ENDDO
1119                ENDDO
1120           
1121             ELSE
1122                DO  i = nxl, nxr
1123                   DO  j = nys, nyn
1124                      DO  k = nzb_do, nzt_do
1125                          local_pf(i,j,k) = MERGE(                             &
1126                                              chem_species(lsp)%conc_av(k,j,i),&
1127                                              REAL( fill_value, KIND = wp ),   &
1128                                              BTEST( wall_flags_0(k,j,i), 0 ) )
1129                      ENDDO
1130                   ENDDO
1131                ENDDO
1132             ENDIF
1133
1134             found = .TRUE.
1135          ENDIF
1136       ENDDO
1137
1138       RETURN
1139    END SUBROUTINE chem_data_output_3d
1140!
1141!------------------------------------------------------------------------------!
1142!
1143! Description:
1144! ------------
1145!> Subroutine defining mask output variables for chemical species
1146!------------------------------------------------------------------------------!
1147   
1148    SUBROUTINE chem_data_output_mask( av, variable, found, local_pf )
1149   
1150       USE control_parameters
1151       USE indices
1152       USE kinds
1153       USE pegrid,             ONLY: myid, threads_per_task
1154       USE surface_mod,        ONLY: get_topography_top_index_ji
1155
1156       IMPLICIT NONE
1157
1158       CHARACTER(LEN=5) ::  grid        !< flag to distinquish between staggered grids
1159
1160       CHARACTER (LEN=*)::  variable    !<
1161       INTEGER(iwp)     ::  av          !< flag to control data output of instantaneous or time-averaged data
1162       LOGICAL          ::  found       !<
1163       REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
1164                 local_pf   !<
1165
1166
1167       !-- local variables.
1168       CHARACTER(len=16)    ::  spec_name
1169       INTEGER(iwp) ::  lsp
1170       INTEGER(iwp) ::  i               !< grid index along x-direction
1171       INTEGER(iwp) ::  j               !< grid index along y-direction
1172       INTEGER(iwp) ::  k               !< grid index along z-direction
1173       INTEGER(iwp) ::  topo_top_ind    !< k index of highest horizontal surface
1174
1175       found = .TRUE.
1176       grid  = 's'
1177
1178       spec_name = TRIM( variable(4:) )
1179   !av == 0
1180
1181       DO lsp=1,nspec
1182          IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name) )               THEN             
1183!
1184!--   todo: remove or replace by "CALL message" mechanism (kanani)
1185!                 IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM(variable)  // &
1186!                                                           TRIM(chem_species(lsp)%name)       
1187             IF (av == 0) THEN
1188                IF ( .NOT. mask_surface(mid) )  THEN
1189
1190                   DO  i = 1, mask_size_l(mid,1)
1191                      DO  j = 1, mask_size_l(mid,2) 
1192                         DO  k = 1, mask_size(mid,3) 
1193                             local_pf(i,j,k) = chem_species(lsp)%conc(  &
1194                                                  mask_k(mid,k),        &
1195                                                  mask_j(mid,j),        &
1196                                                  mask_i(mid,i)      )
1197                         ENDDO
1198                      ENDDO
1199                   ENDDO
1200
1201                ELSE
1202!             
1203!--                Terrain-following masked output
1204                   DO  i = 1, mask_size_l(mid,1)
1205                      DO  j = 1, mask_size_l(mid,2)
1206!             
1207!--                      Get k index of highest horizontal surface
1208                         topo_top_ind = get_topography_top_index_ji( &
1209                                           mask_j(mid,j),  &
1210                                           mask_i(mid,i),  &
1211                                           grid                    )
1212!             
1213!--                      Save output array
1214                         DO  k = 1, mask_size_l(mid,3)
1215                            local_pf(i,j,k) = chem_species(lsp)%conc( &
1216                                                 MIN( topo_top_ind+mask_k(mid,k), &
1217                                                      nzt+1 ),        &
1218                                                 mask_j(mid,j),       &
1219                                                 mask_i(mid,i)      )
1220                         ENDDO
1221                      ENDDO
1222                   ENDDO
1223
1224                ENDIF
1225             ELSE
1226                IF ( .NOT. mask_surface(mid) )  THEN
1227
1228                   DO  i = 1, mask_size_l(mid,1)
1229                      DO  j = 1, mask_size_l(mid,2)
1230                         DO  k =  1, mask_size_l(mid,3)
1231                             local_pf(i,j,k) = chem_species(lsp)%conc_av(  &
1232                                                  mask_k(mid,k),           &
1233                                                  mask_j(mid,j),           &
1234                                                  mask_i(mid,i)         )
1235                         ENDDO
1236                      ENDDO
1237                   ENDDO
1238
1239                ELSE
1240!             
1241!--                Terrain-following masked output
1242                   DO  i = 1, mask_size_l(mid,1)
1243                      DO  j = 1, mask_size_l(mid,2)
1244!             
1245!--                      Get k index of highest horizontal surface
1246                         topo_top_ind = get_topography_top_index_ji( &
1247                                           mask_j(mid,j),  &
1248                                           mask_i(mid,i),  &
1249                                           grid                    )
1250!             
1251!--                      Save output array
1252                         DO  k = 1, mask_size_l(mid,3)
1253                            local_pf(i,j,k) = chem_species(lsp)%conc_av(  &
1254                                                 MIN( topo_top_ind+mask_k(mid,k), &
1255                                                      nzt+1 ),            &
1256                                                 mask_j(mid,j),           &
1257                                                 mask_i(mid,i)         )
1258                         ENDDO
1259                      ENDDO
1260                   ENDDO
1261
1262                ENDIF
1263
1264
1265             ENDIF
1266             found = .FALSE.
1267          ENDIF
1268       ENDDO
1269
1270       RETURN
1271     
1272    END SUBROUTINE chem_data_output_mask     
1273
1274!
1275!------------------------------------------------------------------------------!
1276!
1277! Description:
1278! ------------
1279!> Subroutine defining appropriate grid for netcdf variables.
1280!> It is called out from subroutine netcdf.
1281!------------------------------------------------------------------------------!
1282    SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
1283
1284       IMPLICIT NONE
1285
1286       CHARACTER (LEN=*), INTENT(IN)  ::  var          !<
1287       LOGICAL, INTENT(OUT)           ::  found        !<
1288       CHARACTER (LEN=*), INTENT(OUT) ::  grid_x       !<
1289       CHARACTER (LEN=*), INTENT(OUT) ::  grid_y       !<
1290       CHARACTER (LEN=*), INTENT(OUT) ::  grid_z       !<
1291
1292       found  = .TRUE.
1293
1294       IF (var(1:3) == 'kc_')   THEN                   !< always the same grid for chemistry variables
1295          grid_x = 'x'
1296          grid_y = 'y'
1297          grid_z = 'zu'                             
1298       ELSE
1299          found  = .FALSE.
1300          grid_x = 'none'
1301          grid_y = 'none'
1302          grid_z = 'none'
1303       ENDIF
1304
1305
1306    END SUBROUTINE chem_define_netcdf_grid
1307!
1308!------------------------------------------------------------------------------!
1309!
1310! Description:
1311! ------------
1312!> Subroutine defining header output for chemistry model
1313!------------------------------------------------------------------------------!
1314    SUBROUTINE chem_header ( io )
1315       
1316       IMPLICIT NONE
1317 
1318       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
1319       INTEGER(iwp)             :: lsp            !< running index for chem spcs
1320       INTEGER(iwp)             :: cs_fixed 
1321       CHARACTER (LEN=80)       :: docsflux_chr
1322       CHARACTER (LEN=80)       :: docsinit_chr
1323
1324!
1325!--    Write chemistry model  header
1326       WRITE( io, 1 )
1327
1328!--    Gasphase reaction status
1329       IF ( chem_gasphase_on )  THEN
1330          WRITE( io, 2 )
1331       ELSE
1332          WRITE( io, 3 )
1333       ENDIF
1334
1335!      Chemistry time-step
1336       WRITE ( io, 4 ) cs_time_step
1337
1338!--    Emission mode info
1339       IF ( mode_emis == "DEFAULT" )  THEN
1340          WRITE( io, 5 ) 
1341       ELSEIF ( mode_emis == "PARAMETERIZED" )  THEN
1342          WRITE( io, 6 )
1343       ELSEIF ( mode_emis == "PRE-PROCESSED" )  THEN
1344          WRITE( io, 7 )
1345       ENDIF 
1346
1347!--    Photolysis scheme info
1348       IF ( photolysis_scheme == "simple" )      THEN
1349          WRITE( io, 8 ) 
1350       ELSEIF (photolysis_scheme == "conastant" ) THEN
1351          WRITE( io, 9 )
1352       ENDIF
1353 
1354 !--   Emission flux info
1355       lsp = 1
1356       docsflux_chr ='Chemical species for surface emission flux: ' 
1357       DO WHILE ( surface_csflux_name(lsp) /= 'novalue' )
1358          docsflux_chr = TRIM( docsflux_chr ) // ' ' // TRIM( surface_csflux_name(lsp) ) // ',' 
1359          IF ( LEN_TRIM( docsflux_chr ) >= 75 )  THEN
1360           WRITE ( io, 10 )  docsflux_chr
1361           docsflux_chr = '       '
1362          ENDIF
1363          lsp = lsp + 1
1364       ENDDO
1365 
1366       IF ( docsflux_chr /= '' )  THEN
1367          WRITE ( io, 10 )  docsflux_chr
1368       ENDIF
1369
1370
1371!--    initializatoin of Surface and profile chemical species
1372
1373       lsp = 1
1374       docsinit_chr ='Chemical species for initial surface and profile emissions: ' 
1375       DO WHILE ( cs_name(lsp) /= 'novalue' )
1376          docsinit_chr = TRIM( docsinit_chr ) // ' ' // TRIM( cs_name(lsp) ) // ',' 
1377          IF ( LEN_TRIM( docsinit_chr ) >= 75 )  THEN
1378           WRITE ( io, 11 )  docsinit_chr
1379           docsinit_chr = '       '
1380          ENDIF
1381          lsp = lsp + 1
1382       ENDDO
1383 
1384       IF ( docsinit_chr /= '' )  THEN
1385          WRITE ( io, 11 )  docsinit_chr
1386       ENDIF
1387
1388!--    number of variable and fix chemical species and number of reactions
1389       cs_fixed = nspec - nvar
1390       WRITE ( io, * ) '   --> Chemical species, variable: ', nvar
1391       WRITE ( io, * ) '   --> Chemical species, fixed   : ', cs_fixed
1392       WRITE ( io, * ) '   --> Total number of reactions : ', nreact
1393
1394
13951   FORMAT (//' Chemistry model information:'/                                  &
1396              ' ----------------------------'/)
13972   FORMAT ('    --> Chemical reactions are turned on')
13983   FORMAT ('    --> Chemical reactions are turned off')
13994   FORMAT ('    --> Time-step for chemical species: ',F6.2, ' s')
14005   FORMAT ('    --> Emission mode = DEFAULT ')
14016   FORMAT ('    --> Emission mode = PARAMETERIZED ')
14027   FORMAT ('    --> Emission mode = PRE-PROCESSED ')
14038   FORMAT ('    --> Photolysis scheme used =  simple ')
14049   FORMAT ('    --> Photolysis scheme used =  constant ')
140510  FORMAT (/'    ',A) 
140611  FORMAT (/'    ',A) 
1407!
1408!
1409    END SUBROUTINE chem_header
1410
1411!
1412!------------------------------------------------------------------------------!
1413!
1414! Description:
1415! ------------
1416!> Subroutine initializating chemistry_model_mod
1417!------------------------------------------------------------------------------!
1418    SUBROUTINE chem_init
1419
1420       USE control_parameters,                                                  &
1421          ONLY: message_string, io_blocks, io_group, turbulent_inflow
1422       USE arrays_3d,                                                           &
1423           ONLY: mean_inflow_profiles
1424       USE pegrid
1425
1426       IMPLICIT   none
1427   !-- local variables
1428       INTEGER                 :: i,j               !< running index for for horiz numerical grid points
1429       INTEGER                 :: lsp               !< running index for chem spcs
1430       INTEGER                 :: lpr_lev           !< running index for chem spcs profile level
1431!
1432!-- NOPOINTER version not implemented yet
1433! #if defined( __nopointer )
1434!     message_string = 'The chemistry module only runs with POINTER version'
1435!     CALL message( 'chemistry_model_mod', 'CM0001', 1, 2, 0, 6, 0 )     
1436! #endif
1437!
1438!-- Allocate memory for chemical species
1439       ALLOCATE( chem_species(nspec) )
1440       ALLOCATE( spec_conc_1 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1441       ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1442       ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
1443       ALLOCATE( spec_conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) ) 
1444       ALLOCATE( phot_frequen(nphot) ) 
1445       ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) )
1446       ALLOCATE( bc_cs_t_val(nspec) )
1447!
1448!-- Initialize arrays
1449       spec_conc_1 (:,:,:,:) = 0.0_wp
1450       spec_conc_2 (:,:,:,:) = 0.0_wp
1451       spec_conc_3 (:,:,:,:) = 0.0_wp
1452       spec_conc_av(:,:,:,:) = 0.0_wp
1453
1454
1455       DO lsp = 1, nspec
1456          chem_species(lsp)%name    = spc_names(lsp)
1457
1458          chem_species(lsp)%conc   (nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_1 (:,:,:,lsp)
1459          chem_species(lsp)%conc_p (nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_2 (:,:,:,lsp)
1460          chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_3 (:,:,:,lsp)
1461          chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_av(:,:,:,lsp)     
1462
1463          ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))                   
1464          chem_species(lsp)%cssws_av    = 0.0_wp
1465!
1466!--       The following block can be useful when emission module is not applied. &
1467!--       if emission module is applied the following block will be overwritten.
1468          ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))   
1469          ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))   
1470          ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1)) 
1471          ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))   
1472          chem_species(lsp)%flux_s_cs = 0.0_wp                                     
1473          chem_species(lsp)%flux_l_cs = 0.0_wp                                     
1474          chem_species(lsp)%diss_s_cs = 0.0_wp                                     
1475          chem_species(lsp)%diss_l_cs = 0.0_wp                                     
1476!
1477!--   Allocate memory for initial concentration profiles
1478!--   (concentration values come from namelist)
1479!--   (@todo (FK): Because of this, chem_init is called in palm before
1480!--               check_parameters, since conc_pr_init is used there.
1481!--               We have to find another solution since chem_init should
1482!--               eventually be called from init_3d_model!!)
1483          ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) )
1484          chem_species(lsp)%conc_pr_init(:) = 0.0_wp
1485
1486       ENDDO
1487
1488!
1489!--    Initial concentration of profiles is prescribed by parameters cs_profile
1490!--    and cs_heights in the namelist &chemistry_parameters
1491       CALL chem_init_profiles     
1492           
1493           
1494!
1495!-- Initialize model variables
1496
1497
1498       IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
1499            TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1500
1501
1502!--    First model run of a possible job queue.
1503!--    Initial profiles of the variables must be computed.
1504          IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
1505            CALL location_message( 'initializing with 1D chemistry model profiles', .FALSE. )
1506!
1507!--        Transfer initial profiles to the arrays of the 3D model
1508             DO lsp = 1, nspec
1509                DO  i = nxlg, nxrg
1510                   DO  j = nysg, nyng
1511                      DO lpr_lev = 1, nz + 1
1512                         chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)
1513                      ENDDO
1514                   ENDDO
1515                ENDDO   
1516             ENDDO
1517         
1518          ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
1519          THEN
1520             CALL location_message( 'initializing with constant chemistry profiles', .FALSE. )
1521
1522             DO lsp = 1, nspec 
1523                DO  i = nxlg, nxrg
1524                   DO  j = nysg, nyng
1525                      chem_species(lsp)%conc(:,j,i) = chem_species(lsp)%conc_pr_init   
1526                   ENDDO
1527                ENDDO
1528             ENDDO
1529
1530          ENDIF
1531
1532!
1533!--       If required, change the surface chem spcs at the start of the 3D run
1534          IF ( cs_surface_initial_change(1) /= 0.0_wp ) THEN           
1535             DO lsp = 1, nspec 
1536                chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) +  &
1537                                                  cs_surface_initial_change(lsp)
1538             ENDDO
1539          ENDIF 
1540!
1541!--      Initiale old and new time levels.
1542          DO lsp = 1, nvar
1543             chem_species(lsp)%tconc_m = 0.0_wp                     
1544             chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
1545          ENDDO
1546
1547       ENDIF
1548
1549
1550
1551!---   new code add above this line
1552       DO lsp = 1, nphot
1553          phot_frequen(lsp)%name = phot_names(lsp)
1554!
1555!--   todo: remove or replace by "CALL message" mechanism (kanani)
1556!         IF( myid == 0 )  THEN
1557!            WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,trim(phot_names(lsp))
1558!         ENDIF
1559             phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => freq_1(:,:,:,lsp)
1560       ENDDO
1561
1562       RETURN
1563
1564 END SUBROUTINE chem_init
1565
1566!
1567!------------------------------------------------------------------------------!
1568!
1569! Description:
1570! ------------
1571!> Subroutine defining initial vertical profiles of chemical species (given by
1572!> namelist parameters chem_profiles and chem_heights)  --> which should work
1573!> analogue to parameters u_profile, v_profile and uv_heights)
1574!------------------------------------------------------------------------------!
1575    SUBROUTINE chem_init_profiles              !< SUBROUTINE is called from chem_init in case of
1576                                               !< TRIM( initializing_actions ) /= 'read_restart_data'
1577                                               !< We still need to see what has to be done in case of restart run
1578       USE chem_modules
1579
1580       IMPLICIT NONE
1581
1582   !-- Local variables
1583       INTEGER ::  lsp        !< running index for number of species in derived data type species_def
1584       INTEGER ::  lsp_usr     !< running index for number of species (user defined)  in cs_names, cs_profiles etc
1585       INTEGER ::  lpr_lev    !< running index for profile level for each chem spcs.
1586       INTEGER ::  npr_lev    !< the next available profile lev
1587
1588!-----------------
1589!-- Parameter "cs_profile" and "cs_heights" are used to prescribe user defined initial profiles
1590!-- and heights. If parameter "cs_profile" is not prescribed then initial surface values
1591!-- "cs_surface" are used as constant initial profiles for each species. If "cs_profile" and
1592!-- "cs_heights" are prescribed, their values will!override the constant profile given by
1593!-- "cs_surface".
1594!
1595       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1596          lsp_usr = 1
1597          DO  WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )   !'novalue' is the default
1598             DO  lsp = 1, nspec                                !
1599!--             create initial profile (conc_pr_init) for each chemical species
1600                IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN   !
1601                   IF ( cs_profile(lsp_usr,1) == 9999999.9_wp ) THEN
1602!--                set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
1603                      DO lpr_lev = 0, nzt+1
1604                         chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
1605                      ENDDO
1606                   ELSE
1607                      IF ( cs_heights(1,1) /= 0.0_wp )  THEN
1608                         message_string = 'The surface value of cs_heights must be 0.0'
1609                         CALL message( 'chem_check_parameters', 'CM0434', 1, 2, 0, 6, 0 )
1610                      ENDIF
1611     
1612                      use_prescribed_profile_data = .TRUE.
1613   
1614                      npr_lev = 1
1615!                     chem_species(lsp)%conc_pr_init(0) = 0.0_wp
1616                      DO  lpr_lev = 1, nz+1
1617                         IF ( npr_lev < 100 )  THEN
1618                            DO  WHILE ( cs_heights(lsp_usr, npr_lev+1) <= zu(lpr_lev) )
1619                               npr_lev = npr_lev + 1
1620                               IF ( npr_lev == 100 )  THEN
1621                                   message_string = 'number of chem spcs exceeding the limit'
1622                                   CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )               
1623                                   EXIT
1624                               ENDIF
1625                            ENDDO
1626                         ENDIF
1627                         IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr, npr_lev + 1) /= 9999999.9_wp )  THEN
1628                            chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +                          &
1629                                                    ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                          &
1630                                                    ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) *  &
1631                                                    ( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
1632                         ELSE
1633                               chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev)
1634                         ENDIF
1635                      ENDDO
1636                   ENDIF
1637!-- If a profile is prescribed explicity using cs_profiles and cs_heights, then 
1638!-- chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
1639!-- on the cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:).
1640                ENDIF
1641             ENDDO
1642             lsp_usr = lsp_usr + 1
1643          ENDDO
1644       ENDIF
1645
1646    END SUBROUTINE chem_init_profiles
1647!
1648!------------------------------------------------------------------------------!
1649!
1650! Description:
1651! ------------
1652!> Subroutine to integrate chemical species in the given chemical mechanism
1653!------------------------------------------------------------------------------!
1654
1655    SUBROUTINE chem_integrate_ij (i, j)
1656
1657       USE cpulog,                                                              &
1658          ONLY: cpu_log, log_point
1659       USE statistics,                                                          &
1660          ONLY:  weight_pres
1661        USE control_parameters,                                                 &
1662          ONLY:  dt_3d, intermediate_timestep_count,simulated_time
1663
1664       IMPLICIT   none
1665       INTEGER,INTENT(IN)       :: i,j
1666
1667 !--   local variables
1668       INTEGER(iwp) ::  lsp                                                     !< running index for chem spcs.
1669       INTEGER(iwp) ::  lph                                                     !< running index for photolysis frequencies
1670       INTEGER                  :: k,m,istatf
1671       INTEGER,DIMENSION(20)    :: istatus
1672       REAL(kind=wp),DIMENSION(nzb+1:nzt,nspec)                :: tmp_conc
1673       REAL(kind=wp),DIMENSION(nzb+1:nzt)                      :: tmp_temp
1674       REAL(kind=wp),DIMENSION(nzb+1:nzt)                      :: tmp_qvap
1675       REAL(kind=wp),DIMENSION(nzb+1:nzt,nphot)                :: tmp_phot
1676       REAL(kind=wp),DIMENSION(nzb+1:nzt)                      :: tmp_fact
1677       REAL(kind=wp),DIMENSION(nzb+1:nzt)                      :: tmp_fact_i    !< conversion factor between
1678                                                                                !<    molecules cm^{-3} and ppm
1679
1680       INTEGER,DIMENSION(nzb+1:nzt)                            :: nacc          !< Number of accepted steps
1681       INTEGER,DIMENSION(nzb+1:nzt)                            :: nrej          !< Number of rejected steps
1682
1683       REAL(wp)                         ::  conv                                !< conversion factor
1684       REAL(wp), PARAMETER              ::  ppm2fr  = 1.0e-6_wp                 !< Conversion factor ppm to fraction
1685       REAL(wp), PARAMETER              ::  xm_air  = 28.96_wp                  !< Mole mass of dry air
1686       REAL(wp), PARAMETER              ::  xm_h2o  = 18.01528_wp               !< Mole mass of water vapor
1687       REAL(wp), PARAMETER              ::  pref_i  = 1._wp / 100000.0_wp       !< inverse reference pressure (1/Pa)
1688       REAL(wp), PARAMETER              ::  t_std   = 273.15_wp                 !< standard pressure (Pa)
1689       REAL(wp), PARAMETER              ::  p_std   = 101325.0_wp               !< standard pressure (Pa)
1690       REAL(wp), PARAMETER              ::  vmolcm  = 22.414e3_wp               !< Mole volume (22.414 l) in cm^{-3}
1691       REAL(wp), PARAMETER              ::  xna     = 6.022e23_wp               !< Avogadro number (molecules/mol)
1692
1693       REAL(wp),DIMENSION(size(rcntrl)) :: rcntrl_local
1694
1695
1696       REAL(kind=wp)  :: dt_chem                                             
1697
1698       CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'start' )
1699!<     set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry
1700       IF (chem_gasphase_on) THEN
1701          nacc = 0
1702          nrej = 0
1703
1704       tmp_temp(:) = pt(nzb+1:nzt,j,i) * exner(nzb+1:nzt)
1705!         ppm to molecules/cm**3
1706!         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 ) 
1707          conv = ppm2fr * xna / vmolcm
1708          tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std)
1709          tmp_fact_i = 1.0_wp/tmp_fact
1710
1711          IF ( humidity ) THEN
1712             IF ( bulk_cloud_model )  THEN
1713                tmp_qvap(:) = ( q(nzb+1:nzt,j,i) - ql(nzb+1:nzt,j,i) ) * xm_air/xm_h2o * tmp_fact(:)
1714             ELSE
1715                tmp_qvap(:) = q(nzb+1:nzt,j,i) * xm_air/xm_h2o * tmp_fact(:)
1716             ENDIF
1717          ELSE
1718             tmp_qvap(:) = 0.01 * tmp_fact(:)                          !< Constant value for q if water vapor is not computed
1719          ENDIF
1720
1721          DO lsp = 1,nspec
1722             tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:) 
1723          ENDDO
1724
1725          DO lph = 1,nphot
1726             tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)               
1727          ENDDO
1728!
1729!--   todo: remove (kanani)
1730!           IF(myid == 0 .AND. chem_debug0 ) THEN
1731!              IF (i == 10 .and. j == 10) WRITE(0,*) 'begin chemics step ',dt_3d
1732!           ENDIF
1733
1734!--       Compute length of time step
1735          IF ( call_chem_at_all_substeps )  THEN
1736             dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
1737          ELSE
1738             dt_chem = dt_3d
1739          ENDIF
1740
1741          cs_time_step = dt_chem
1742
1743          CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'start' )
1744
1745          IF(maxval(rcntrl) > 0.0)   THEN    ! Only if rcntrl is set
1746             IF( simulated_time <= 2*dt_3d)  THEN
1747                 rcntrl_local = 0
1748!
1749!--   todo: remove (kanani)
1750!                  WRITE(9,'(a,2f10.3)') 'USE Default rcntrl in the first steps ',simulated_time,dt_3d
1751             ELSE
1752                 rcntrl_local = rcntrl
1753             ENDIF
1754          ELSE
1755             rcntrl_local = 0
1756          END IF
1757
1758          CALL chem_gasphase_integrate (dt_chem, tmp_conc, tmp_temp, tmp_qvap, tmp_fact, tmp_phot, &
1759                             icntrl_i = icntrl, rcntrl_i = rcntrl_local, xnacc = nacc, xnrej = nrej, istatus=istatus)
1760
1761          CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'stop' )
1762
1763          DO lsp = 1,nspec
1764             chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:)
1765          ENDDO
1766
1767
1768       ENDIF
1769          CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'stop' )
1770
1771       RETURN
1772    END SUBROUTINE chem_integrate_ij
1773!
1774!------------------------------------------------------------------------------!
1775!
1776! Description:
1777! ------------
1778!> Subroutine defining parin for &chemistry_parameters for chemistry model
1779!------------------------------------------------------------------------------!
1780    SUBROUTINE chem_parin
1781   
1782       USE chem_modules
1783       USE control_parameters
1784     
1785       USE kinds
1786       USE pegrid
1787       USE statistics
1788           
1789       IMPLICIT   NONE
1790
1791       CHARACTER (LEN=80) ::  line                        !< dummy string that contains the current line of the parameter file
1792       CHARACTER (LEN=3)  ::  cs_prefix 
1793
1794       REAL(wp), DIMENSION(nmaxfixsteps) ::   my_steps    !< List of fixed timesteps   my_step(1) = 0.0 automatic stepping
1795       INTEGER(iwp) ::  i                                 !<
1796       INTEGER(iwp) ::  j                                 !<
1797       INTEGER(iwp) ::  max_pr_cs_tmp                     !<
1798
1799
1800       NAMELIST /chemistry_parameters/  bc_cs_b,                          &
1801                                        bc_cs_t,                          &
1802                                        call_chem_at_all_substeps,        &
1803                                        chem_debug0,                      &
1804                                        chem_debug1,                      &
1805                                        chem_debug2,                      &
1806                                        chem_gasphase_on,                 &
1807                                        cs_heights,                       &
1808                                        cs_name,                          &
1809                                        cs_profile,                       &
1810                                        cs_profile_name,                  & 
1811                                        cs_surface,                       &
1812                                        decycle_chem_lr,                  &
1813                                        decycle_chem_ns,                  &           
1814                                        decycle_method,                   & 
1815                                        emiss_factor_main,                &
1816                                        emiss_factor_side,                &                     
1817                                        icntrl,                           &
1818                                        main_street_id,                   &
1819                                        max_street_id,                    &
1820                                        my_steps,                         &
1821                                        nest_chemistry,                   &
1822                                        rcntrl,                           &
1823                                        side_street_id,                   &
1824                                        photolysis_scheme,                &
1825                                        wall_csflux,                      &
1826                                        cs_vertical_gradient,             &
1827                                        top_csflux,                       & 
1828                                        surface_csflux,                   &
1829                                        surface_csflux_name,              &
1830                                        cs_surface_initial_change,        &
1831                                        cs_vertical_gradient_level,       &
1832!                                       namelist parameters for emissions
1833                                        mode_emis,                        &
1834                                        time_fac_type,                    &
1835                                        daytype_mdh,                      &
1836                                        do_emis                             
1837                             
1838!-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
1839!-- so this way we could prescribe a specific flux value for each species
1840!>  chemistry_parameters for initial profiles
1841!>  cs_names = 'O3', 'NO2', 'NO', ...   to set initial profiles)
1842!>  cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3)
1843!>  cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.)
1844!>  cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, .....  (chem spcs conc at height lvls chem_heights(1,:)) etc.
1845!>  If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)"
1846!>  then write these cs_surface values to chem_species(lsp)%conc_pr_init(:)
1847
1848!
1849!--   Read chem namelist   
1850
1851       INTEGER             :: ier
1852       CHARACTER(LEN=64)   :: text
1853       CHARACTER(LEN=8)    :: solver_type
1854
1855       icntrl    = 0
1856       rcntrl    = 0.0_wp
1857       my_steps  = 0.0_wp
1858       photolysis_scheme = 'simple'
1859       atol = 1.0_wp
1860       rtol = 0.01_wp
1861 !
1862 !--   Try to find chemistry package
1863       REWIND ( 11 )
1864       line = ' '
1865       DO   WHILE ( INDEX( line, '&chemistry_parameters' ) == 0 )
1866          READ ( 11, '(A)', END=20 )  line
1867       ENDDO
1868       BACKSPACE ( 11 )
1869 !
1870 !--   Read chemistry namelist
1871       READ ( 11, chemistry_parameters, ERR = 10, END = 20 )     
1872 !
1873 !--   Enable chemistry model
1874       air_chemistry = .TRUE.                   
1875      GOTO 20
1876
1877 10   BACKSPACE( 11 )
1878      READ( 11 , '(A)') line
1879      CALL parin_fail_message( 'chemistry_parameters', line )
1880
1881 20   CONTINUE
1882
1883!
1884!--    check for  emission mode for chem species
1885       IF ( (mode_emis /= 'PARAMETERIZED')  .AND. ( mode_emis /= 'DEFAULT') .AND. (mode_emis /= 'PRE-PROCESSED'  ) )   THEN
1886            message_string = 'Incorrect mode_emiss  option select. Please check spelling'
1887            CALL message( 'chem_check_parameters', 'CM0436', 1, 2, 0, 6, 0 )
1888       ENDIF
1889
1890       t_steps = my_steps         
1891                                   
1892!--    Determine the number of user-defined profiles and append them to the
1893!--    standard data output (data_output_pr)
1894       max_pr_cs_tmp = 0 
1895       i = 1
1896
1897       DO  WHILE ( data_output_pr(i)  /= ' '  .AND.  i <= 100 )
1898          IF ( TRIM( data_output_pr(i)(1:3)) == 'kc_' ) THEN
1899             max_pr_cs_tmp = max_pr_cs_tmp+1
1900          ENDIF
1901          i = i +1
1902       ENDDO
1903
1904       IF ( max_pr_cs_tmp > 0 ) THEN
1905          cs_pr_namelist_found = .TRUE.
1906          max_pr_cs = max_pr_cs_tmp
1907       ENDIF
1908
1909!      Set Solver Type
1910       IF(icntrl(3) == 0)   THEN
1911          solver_type = 'rodas3'           !Default
1912       ELSE IF(icntrl(3) == 1)   THEN
1913          solver_type = 'ros2'
1914       ELSE IF(icntrl(3) == 2)   THEN
1915          solver_type = 'ros3'
1916       ELSE IF(icntrl(3) == 3)   THEN
1917          solver_type = 'ro4'
1918       ELSE IF(icntrl(3) == 4)   THEN
1919          solver_type = 'rodas3'
1920       ELSE IF(icntrl(3) == 5)   THEN
1921          solver_type = 'rodas4'
1922       ELSE IF(icntrl(3) == 6)   THEN
1923          solver_type = 'Rang3'
1924       ELSE
1925           message_string = 'illegal solver type'
1926           CALL message( 'chem_parin', 'PA0506', 1, 2, 0, 6, 0 )
1927       END IF
1928
1929!
1930!--   todo: remove or replace by "CALL message" mechanism (kanani)
1931!       write(text,*) 'gas_phase chemistry: solver_type = ',trim(solver_type)
1932!kk    Has to be changed to right calling sequence
1933!kk       CALL location_message( trim(text), .FALSE. )
1934!        IF(myid == 0)   THEN
1935!           write(9,*) ' '
1936!           write(9,*) 'kpp setup '
1937!           write(9,*) ' '
1938!           write(9,*) '    gas_phase chemistry: solver_type = ',trim(solver_type)
1939!           write(9,*) ' '
1940!           write(9,*) '    Hstart  = ',rcntrl(3)
1941!           write(9,*) '    FacMin  = ',rcntrl(4)
1942!           write(9,*) '    FacMax  = ',rcntrl(5)
1943!           write(9,*) ' '
1944!           IF(vl_dim > 1)    THEN
1945!              write(9,*) '    Vector mode                   vektor length = ',vl_dim
1946!           ELSE
1947!              write(9,*) '    Scalar mode'
1948!           ENDIF
1949!           write(9,*) ' '
1950!        END IF
1951
1952       RETURN
1953
1954    END SUBROUTINE chem_parin
1955
1956!
1957!------------------------------------------------------------------------------!
1958!
1959! Description:
1960! ------------
1961!> Subroutine calculating prognostic equations for chemical species
1962!> (vector-optimized).
1963!> Routine is called separately for each chemical species over a loop from
1964!> prognostic_equations.
1965!------------------------------------------------------------------------------!
1966    SUBROUTINE chem_prognostic_equations ( cs_scalar_p, cs_scalar, tcs_scalar_m,  &
1967                                           pr_init_cs, ilsp )
1968
1969       USE advec_s_pw_mod,                                                        &
1970           ONLY:  advec_s_pw
1971       USE advec_s_up_mod,                                                        &
1972           ONLY:  advec_s_up
1973       USE advec_ws,                                                              &
1974           ONLY:  advec_s_ws 
1975       USE diffusion_s_mod,                                                       &
1976           ONLY:  diffusion_s
1977       USE indices,                                                               &
1978           ONLY:  nxl, nxr, nyn, nys, wall_flags_0
1979       USE pegrid
1980       USE surface_mod,                                                           &
1981           ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
1982                  surf_usm_v
1983
1984       IMPLICIT NONE
1985
1986       INTEGER ::  i   !< running index
1987       INTEGER ::  j   !< running index
1988       INTEGER ::  k   !< running index
1989
1990       INTEGER(iwp),INTENT(IN) ::  ilsp          !<
1991
1992       REAL(wp), DIMENSION(0:nz+1) ::  pr_init_cs   !<
1993
1994       REAL(wp), DIMENSION(:,:,:), POINTER ::  cs_scalar      !<
1995       REAL(wp), DIMENSION(:,:,:), POINTER ::  cs_scalar_p    !<
1996       REAL(wp), DIMENSION(:,:,:), POINTER ::  tcs_scalar_m   !<
1997
1998
1999   !
2000   !-- Tendency terms for chemical species
2001       tend = 0.0_wp
2002   !   
2003   !-- Advection terms
2004       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2005          IF ( ws_scheme_sca )  THEN
2006             CALL advec_s_ws( cs_scalar, 'kc' )
2007          ELSE
2008             CALL advec_s_pw( cs_scalar )
2009          ENDIF
2010       ELSE
2011            CALL advec_s_up( cs_scalar )
2012       ENDIF
2013   !
2014   !-- Diffusion terms  (the last three arguments are zero)
2015       CALL diffusion_s( cs_scalar,                                               &
2016                         surf_def_h(0)%cssws(ilsp,:),                             &
2017                         surf_def_h(1)%cssws(ilsp,:),                             &
2018                         surf_def_h(2)%cssws(ilsp,:),                             &
2019                         surf_lsm_h%cssws(ilsp,:),                                &
2020                         surf_usm_h%cssws(ilsp,:),                                &
2021                         surf_def_v(0)%cssws(ilsp,:),                             &
2022                         surf_def_v(1)%cssws(ilsp,:),                             &
2023                         surf_def_v(2)%cssws(ilsp,:),                             &
2024                         surf_def_v(3)%cssws(ilsp,:),                             &
2025                         surf_lsm_v(0)%cssws(ilsp,:),                             &
2026                         surf_lsm_v(1)%cssws(ilsp,:),                             &
2027                         surf_lsm_v(2)%cssws(ilsp,:),                             &
2028                         surf_lsm_v(3)%cssws(ilsp,:),                             &
2029                         surf_usm_v(0)%cssws(ilsp,:),                             &
2030                         surf_usm_v(1)%cssws(ilsp,:),                             &
2031                         surf_usm_v(2)%cssws(ilsp,:),                             &
2032                         surf_usm_v(3)%cssws(ilsp,:) )
2033   !   
2034   !-- Prognostic equation for chemical species
2035       DO  i = nxl, nxr
2036          DO  j = nys, nyn   
2037             DO  k = nzb+1, nzt
2038                cs_scalar_p(k,j,i) =   cs_scalar(k,j,i)                           &
2039                                     + ( dt_3d  *                                 &
2040                                         (   tsc(2) * tend(k,j,i)                 &
2041                                           + tsc(3) * tcs_scalar_m(k,j,i)         &
2042                                         )                                        & 
2043                                       - tsc(5) * rdf_sc(k)                       &
2044                                                * ( cs_scalar(k,j,i) - pr_init_cs(k) )    &   
2045                                       )                                          &
2046                                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )       
2047
2048                IF ( cs_scalar_p(k,j,i) < 0.0_wp )  cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i)
2049             ENDDO
2050          ENDDO
2051       ENDDO
2052   !
2053   !-- Calculate tendencies for the next Runge-Kutta step
2054       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2055          IF ( intermediate_timestep_count == 1 )  THEN
2056             DO  i = nxl, nxr
2057                DO  j = nys, nyn   
2058                   DO  k = nzb+1, nzt
2059                      tcs_scalar_m(k,j,i) = tend(k,j,i)
2060                   ENDDO
2061                ENDDO
2062             ENDDO
2063          ELSEIF ( intermediate_timestep_count < &
2064             intermediate_timestep_count_max )  THEN
2065             DO  i = nxl, nxr
2066                DO  j = nys, nyn
2067                   DO  k = nzb+1, nzt
2068                      tcs_scalar_m(k,j,i) = - 9.5625_wp * tend(k,j,i)             &
2069                                            + 5.3125_wp * tcs_scalar_m(k,j,i)
2070                   ENDDO
2071                ENDDO
2072             ENDDO
2073          ENDIF
2074       ENDIF
2075
2076    END SUBROUTINE chem_prognostic_equations
2077
2078!------------------------------------------------------------------------------!
2079!
2080! Description:
2081! ------------
2082!> Subroutine calculating prognostic equations for chemical species
2083!> (cache-optimized).
2084!> Routine is called separately for each chemical species over a loop from
2085!> prognostic_equations.
2086!------------------------------------------------------------------------------!
2087    SUBROUTINE chem_prognostic_equations_ij ( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs,    &
2088                               i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs,                  &
2089                               flux_l_cs, diss_l_cs )
2090       USE pegrid         
2091       USE advec_ws,        ONLY:  advec_s_ws 
2092       USE advec_s_pw_mod,  ONLY:  advec_s_pw
2093       USE advec_s_up_mod,  ONLY:  advec_s_up
2094       USE diffusion_s_mod, ONLY:  diffusion_s
2095       USE indices,         ONLY:  wall_flags_0
2096       USE surface_mod,     ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
2097                                   surf_usm_v
2098
2099
2100       IMPLICIT NONE
2101
2102       REAL(wp), DIMENSION(:,:,:), POINTER   :: cs_scalar_p, cs_scalar, tcs_scalar_m
2103
2104       INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn, ilsp
2105       REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1)         :: flux_s_cs   !<
2106       REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1)         :: diss_s_cs   !<
2107       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: flux_l_cs   !<
2108       REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: diss_l_cs   !<
2109       REAL(wp), DIMENSION(0:nz+1)                                 :: pr_init_cs  !<
2110
2111!-- local variables
2112
2113       INTEGER :: k
2114!
2115!--    Tendency-terms for chem spcs.
2116       tend(:,j,i) = 0.0_wp
2117!   
2118!-- Advection terms
2119       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2120          IF ( ws_scheme_sca )  THEN
2121             CALL advec_s_ws( i, j, cs_scalar, 'kc', flux_s_cs, diss_s_cs,                &
2122                flux_l_cs, diss_l_cs, i_omp_start, tn )
2123          ELSE
2124             CALL advec_s_pw( i, j, cs_scalar )
2125          ENDIF
2126       ELSE
2127            CALL advec_s_up( i, j, cs_scalar )
2128       ENDIF
2129
2130!
2131
2132!-- Diffusion terms (the last three arguments are zero)
2133
2134         CALL diffusion_s( i, j, cs_scalar,                                                 &
2135                           surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:),        &
2136                           surf_def_h(2)%cssws(ilsp,:),                                     &
2137                           surf_lsm_h%cssws(ilsp,:), surf_usm_h%cssws(ilsp,:),              &
2138                           surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:),        &
2139                           surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:),        &
2140                           surf_lsm_v(0)%cssws(ilsp,:), surf_lsm_v(1)%cssws(ilsp,:),        &
2141                           surf_lsm_v(2)%cssws(ilsp,:), surf_lsm_v(3)%cssws(ilsp,:),        &
2142                           surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:),        &
2143                           surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) )
2144   
2145!   
2146!-- Prognostic equation for chem spcs
2147       DO k = nzb+1, nzt
2148          cs_scalar_p(k,j,i) = cs_scalar(k,j,i) + ( dt_3d  *                       &
2149                                                  ( tsc(2) * tend(k,j,i) +         &
2150                                                    tsc(3) * tcs_scalar_m(k,j,i) ) & 
2151                                                  - tsc(5) * rdf_sc(k)             &
2152                                                           * ( cs_scalar(k,j,i) - pr_init_cs(k) )    &   
2153                                                  )                                                  &
2154                                                           * MERGE( 1.0_wp, 0.0_wp,                  &     
2155                                                                   BTEST( wall_flags_0(k,j,i), 0 )   &             
2156                                                                    )       
2157
2158          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
2159       ENDDO
2160
2161!
2162!-- Calculate tendencies for the next Runge-Kutta step
2163       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2164          IF ( intermediate_timestep_count == 1 )  THEN
2165             DO  k = nzb+1, nzt
2166                tcs_scalar_m(k,j,i) = tend(k,j,i)
2167             ENDDO
2168          ELSEIF ( intermediate_timestep_count < &
2169             intermediate_timestep_count_max )  THEN
2170             DO  k = nzb+1, nzt
2171                tcs_scalar_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
2172                   5.3125_wp * tcs_scalar_m(k,j,i)
2173             ENDDO
2174          ENDIF
2175       ENDIF
2176
2177    END SUBROUTINE chem_prognostic_equations_ij
2178
2179!
2180!------------------------------------------------------------------------------!
2181!
2182! Description:
2183! ------------
2184!> Subroutine to read restart data of chemical species
2185!------------------------------------------------------------------------------!
2186
2187    SUBROUTINE chem_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,         &
2188                               nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc,  &
2189                               nys_on_file, tmp_3d, found )   
2190                                       
2191       USE control_parameters
2192               
2193       USE indices
2194       
2195       USE pegrid
2196
2197       IMPLICIT NONE
2198
2199       CHARACTER (LEN=20) :: spc_name_av !<   
2200         
2201       INTEGER(iwp) ::  i, lsp          !<
2202       INTEGER(iwp) ::  k               !<
2203       INTEGER(iwp) ::  nxlc            !<
2204       INTEGER(iwp) ::  nxlf            !<
2205       INTEGER(iwp) ::  nxl_on_file     !<   
2206       INTEGER(iwp) ::  nxrc            !<
2207       INTEGER(iwp) ::  nxrf            !<
2208       INTEGER(iwp) ::  nxr_on_file     !<   
2209       INTEGER(iwp) ::  nync            !<
2210       INTEGER(iwp) ::  nynf            !<
2211       INTEGER(iwp) ::  nyn_on_file     !<   
2212       INTEGER(iwp) ::  nysc            !<
2213       INTEGER(iwp) ::  nysf            !<
2214       INTEGER(iwp) ::  nys_on_file     !<   
2215       
2216       LOGICAL, INTENT(OUT)  :: found 
2217
2218       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
2219
2220
2221       found = .FALSE. 
2222
2223
2224          IF ( ALLOCATED(chem_species) )  THEN
2225
2226             DO  lsp = 1, nspec
2227
2228                 !< for time-averaged chemical conc.
2229                spc_name_av  =  TRIM(chem_species(lsp)%name)//'_av'
2230
2231                IF (restart_string(1:length) == TRIM(chem_species(lsp)%name) )    &
2232                THEN
2233                   !< read data into tmp_3d
2234                   IF ( k == 1 )  READ ( 13 )  tmp_3d 
2235                   !< fill ..%conc in the restart run   
2236                   chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp,                  &
2237                                          nxlc-nbgp:nxrc+nbgp) =                  & 
2238                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2239                   found = .TRUE.
2240                ELSEIF (restart_string(1:length) == spc_name_av )  THEN
2241                   IF ( k == 1 )  READ ( 13 )  tmp_3d
2242                   chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp,               &
2243                                             nxlc-nbgp:nxrc+nbgp) =               &
2244                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2245                   found = .TRUE.
2246                ENDIF
2247
2248             ENDDO
2249
2250          ENDIF
2251
2252
2253    END SUBROUTINE chem_rrd_local
2254!
2255!-------------------------------------------------------------------------------!
2256!> Description:
2257!> Calculation of horizontally averaged profiles
2258!> This routine is called for every statistic region (sr) defined by the user,
2259!> but at least for the region "total domain" (sr=0).
2260!> quantities.
2261!------------------------------------------------------------------------------!
2262    SUBROUTINE chem_statistics( mode, sr, tn )
2263
2264
2265       USE arrays_3d
2266       USE indices
2267       USE kinds
2268       USE pegrid
2269       USE statistics
2270
2271    USE user
2272
2273    IMPLICIT NONE
2274
2275    CHARACTER (LEN=*) ::  mode   !<
2276
2277    INTEGER(iwp) :: i    !< running index on x-axis
2278    INTEGER(iwp) :: j    !< running index on y-axis
2279    INTEGER(iwp) :: k    !< vertical index counter
2280    INTEGER(iwp) :: sr   !< statistical region
2281    INTEGER(iwp) :: tn   !< thread number
2282    INTEGER(iwp) :: n    !<
2283    INTEGER(iwp) :: m    !<
2284    INTEGER(iwp) :: lpr  !< running index chem spcs
2285!    REAL(wp),                                                                                      &
2286!    DIMENSION(dots_num_palm+1:dots_max) ::                                                         &
2287!          ts_value_l   !<
2288
2289    IF ( mode == 'profiles' )  THEN
2290!
2291!--    Sample on how to calculate horizontally averaged profiles of user-
2292!--    defined quantities. Each quantity is identified by the index
2293!--    "pr_palm+#" where "#" is an integer starting from 1. These
2294!--    user-profile-numbers must also be assigned to the respective strings
2295!--    given by data_output_pr_cs in routine user_check_data_output_pr.
2296!--    hom(:,:,:,:) =  dim-1 = vertical level, dim-2= 1: met-species,2:zu/zw, dim-3 = quantity( e.g.
2297!                       w*pt*), dim-4 = statistical region.
2298
2299       !$OMP DO
2300       DO  i = nxl, nxr
2301          DO  j = nys, nyn
2302              DO  k = nzb, nzt+1
2303                DO lpr = 1, cs_pr_count
2304
2305                   sums_l(k,pr_palm+lpr,tn) = sums_l(k,pr_palm+lpr,tn) +    &
2306                                                 chem_species(cs_pr_index(lpr))%conc(k,j,i) *       &
2307                                                 rmask(j,i,sr)  *                                   &
2308                                                 MERGE( 1.0_wp, 0.0_wp,                             &
2309                                                 BTEST( wall_flags_0(k,j,i), 22 ) )
2310                ENDDO                                                                         
2311             ENDDO
2312          ENDDO
2313       ENDDO
2314    ELSEIF ( mode == 'time_series' )  THEN
2315       CALL location_message( 'Time series not calculated for chemistry', .TRUE. )
2316    ENDIF
2317
2318    END SUBROUTINE chem_statistics
2319!
2320!------------------------------------------------------------------------------!
2321!
2322! Description:
2323! ------------
2324!> Subroutine for swapping of timelevels for chemical species
2325!> called out from subroutine swap_timelevel
2326!------------------------------------------------------------------------------!
2327
2328    SUBROUTINE chem_swap_timelevel (level)
2329
2330          IMPLICIT   none
2331
2332          INTEGER,INTENT(IN)        :: level
2333!--       local variables
2334          INTEGER                   :: lsp
2335
2336
2337          IF ( level == 0 ) THEN
2338             DO lsp=1, nvar                                       
2339                chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_1(:,:,:,lsp)
2340                chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_2(:,:,:,lsp)
2341             ENDDO
2342          ELSE
2343             DO lsp=1, nvar                                       
2344                chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_2(:,:,:,lsp)
2345                chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_1(:,:,:,lsp)
2346             ENDDO
2347          ENDIF
2348
2349       RETURN
2350    END SUBROUTINE chem_swap_timelevel
2351!
2352!------------------------------------------------------------------------------!
2353!
2354! Description:
2355! ------------
2356!> Subroutine to write restart data for chemistry model
2357!------------------------------------------------------------------------------!
2358    SUBROUTINE chem_wrd_local
2359
2360       IMPLICIT NONE
2361
2362       INTEGER(iwp) ::  lsp !<
2363
2364       DO  lsp = 1, nspec
2365          CALL wrd_write_string( TRIM( chem_species(lsp)%name ))
2366          WRITE ( 14 )  chem_species(lsp)%conc
2367          CALL wrd_write_string( TRIM( chem_species(lsp)%name )//'_av' )
2368          WRITE ( 14 )  chem_species(lsp)%conc_av
2369       ENDDO
2370
2371    END SUBROUTINE chem_wrd_local
2372 END MODULE chemistry_model_mod
2373
Note: See TracBrowser for help on using the repository browser.