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

Last change on this file since 3348 was 3318, checked in by sward, 6 years ago

Fixed faulty syntax in message string

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