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

Last change on this file since 3449 was 3449, checked in by suehring, 4 years ago

Branch resler -r 3439 re-integrated into current trunk: RTM 3.0, transpiration of plant canopy, output fixes in USM

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