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

Last change on this file since 3298 was 3298, checked in by kanani, 3 years ago

Merge chemistry branch at r3297 to trunk

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