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

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

Fix cpp directives and error messages

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