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

Last change on this file since 3254 was 3248, checked in by sward, 6 years ago

Minor format changes

  • Property svn:keywords set to Id
File size: 75.6 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 1997-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 3248 2018-09-14 09:42:06Z suehring $
29! Minor formating changes
30!
31! 3246 2018-09-13 15:14:50Z sward
32! Added error handling for input namelist via parin_fail_message
33!
34! 3241 2018-09-12 15:02:00Z raasch
35! +nest_chemistry
36!
37! 3209 2018-08-27 16:58:37Z suehring
38! Rename flags indicating outflow boundary conditions
39!
40! 3182 2018-07-27 13:36:03Z suehring
41! Revise output of surface quantities in case of overhanging structures
42!
43! 3045 2018-05-28 07:55:41Z Giersch
44! error messages revised
45!
46! 3014 2018-05-09 08:42:38Z maronga
47! Bugfix: nzb_do and nzt_do were not used for 3d data output
48!
49! 3004 2018-04-27 12:33:25Z Giersch
50! Comment concerning averaged data output added
51!
52! 2932 2018-03-26 09:39:22Z maronga
53! renamed chemistry_par to chemistry_parameters
54!
55! 2894 2018-03-15 09:17:58Z Giersch
56! Calculations of the index range of the subdomain on file which overlaps with
57! the current subdomain are already done in read_restart_data_mod,
58! chem_last_actions was renamed to chem_wrd_local, chem_read_restart_data was
59! renamed to chem_rrd_local, chem_write_var_list was renamed to
60! chem_wrd_global, chem_read_var_list was renamed to chem_rrd_global,
61! chem_skip_var_list has been removed, variable named found has been
62! introduced for checking if restart data was found, reading of restart strings
63! has been moved completely to read_restart_data_mod, chem_rrd_local is already
64! inside the overlap loop programmed in read_restart_data_mod, todo list has
65! bees extended, redundant characters in chem_wrd_local have been removed,
66! the marker *** end chemistry *** is not necessary anymore, strings and their
67! respective lengths are written out and read now in case of restart runs to
68! get rid of prescribed character lengths
69!
70! 2815 2018-02-19 11:29:57Z suehring
71! Bugfix in restart mechanism,
72! rename chem_tendency to chem_prognostic_equations,
73! implement vector-optimized version of chem_prognostic_equations,
74! some clean up (incl. todo list)
75!
76! 2773 2018-01-30 14:12:54Z suehring
77! Declare variables required for nesting as public
78!
79! 2772 2018-01-29 13:10:35Z suehring
80! Bugfix in string handling
81!
82! 2768 2018-01-24 15:38:29Z kanani
83! Shorten lines to maximum length of 132 characters
84!
85! 2766 2018-01-22 17:17:47Z kanani
86! Removed preprocessor directive __chem
87!
88! 2756 2018-01-16 18:11:14Z suehring
89! Fill values in 3D output introduced.
90!
91! 2718 2018-01-02 08:49:38Z maronga
92! Initial revision
93!
94!
95!
96!
97! Authors:
98! --------
99! @author Renate Forkel
100! @author Farah Kanani-Suehring
101! @author Klaus Ketelsen
102! @author Basit Khan
103!
104!
105!------------------------------------------------------------------------------!
106! Description:
107! ------------
108!> Chemistry model for PALM-4U
109!> @todo Adjust chem_rrd_local to CASE structure of others modules. It is not
110!>       allowed to use the chemistry model in a precursor run and additionally
111!>       not using it in a main run
112!> @todo Update/clean-up todo list! (FK)
113!> @todo Set proper fill values (/= 0) for chem output arrays! (FK)
114!> @todo Add routine chem_check_parameters, add checks for inconsistent or
115!>       unallowed parameter settings.
116!>       CALL of chem_check_parameters from check_parameters. (FK)
117!> @todo Make routine chem_header available, CALL from header.f90
118!>       (see e.g. how it is done in routine lsm_header in
119!>        land_surface_model_mod.f90). chem_header should include all setup
120!>        info about chemistry parameter settings. (FK)
121!> @todo Implement turbulent inflow of chem spcs in inflow_turbulence.
122!> @todo Separate boundary conditions for each chem spcs to be implemented
123!> @todo Currently only total concentration are calculated. Resolved, parameterized
124!>       and chemistry fluxes although partially and some completely coded but
125!>       are not operational/activated in this version. bK.
126!> @todo slight differences in passive scalar and chem spcs when chem reactions
127!>       turned off. Need to be fixed. bK
128!> @todo test nesting for chem spcs, was implemented by suehring (kanani)
129!> @todo subroutine set_const_initial_values to be taken out from chemistry_model_mod !bK.
130!> @todo chemistry error messages
131!> @todo Format this module according to PALM coding standard (see e.g. module
132!>       template under http://palm.muk.uni-hannover.de/mosaik/downloads/8 or
133!>       D3_coding_standard.pdf under https://palm.muk.uni-hannover.de/trac/downloads/16)
134!
135!------------------------------------------------------------------------------!
136
137MODULE chemistry_model_mod
138
139   USE kinds,              ONLY: wp, iwp
140   USE indices,            ONLY: nz, nzb,nzt,nysg,nyng,nxlg,nxrg,nys,nyn
141   USE pegrid,             ONLY: myid, threads_per_task
142   USE control_parameters, ONLY: dt_3d, ws_scheme_sca, initializing_actions, message_string, &
143                                 omega, tsc, intermediate_timestep_count,      &
144                                 intermediate_timestep_count_max,              &
145                                 timestep_scheme, use_prescribed_profile_data 
146   USE arrays_3d,          ONLY: hyp, pt, rdf_sc, tend, zu                     
147   USE chem_gasphase_mod,  ONLY: nspec, spc_names, nkppctrl, nmaxfixsteps,     &
148                                 t_steps, fill_temp, chem_gasphase_integrate,  &
149                                 nvar, atol, rtol, nphot, phot_names
150   USE cpulog,             ONLY: cpu_log, log_point
151
152   USE chem_modules
153 
154
155   IMPLICIT   NONE
156   PRIVATE
157   SAVE
158
159   LOGICAL ::  nest_chemistry = .TRUE. !< flag for nesting mode of chemical species, independent on parent or not
160!
161!- Define chemical variables
162   TYPE   species_def
163      CHARACTER(LEN=8)                                   :: name
164      CHARACTER(LEN=16)                                  :: unit
165      REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: conc
166      REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: conc_av
167      REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: conc_p
168      REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: tconc_m
169      REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:)           :: cssws_av           
170      REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:)           :: flux_s_cs
171      REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:)           :: diss_s_cs
172      REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:)         :: flux_l_cs
173      REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:)         :: diss_l_cs
174      REAL(kind=wp),ALLOCATABLE,DIMENSION(:)             :: conc_pr_init
175   END TYPE species_def
176
177
178   TYPE   photols_def                                                           
179      CHARACTER(LEN=8)                                   :: name
180      CHARACTER(LEN=16)                                  :: unit
181      REAL(kind=wp),POINTER,DIMENSION(:,:,:)             :: freq
182   END TYPE photols_def
183
184
185   PUBLIC  species_def
186   PUBLIC  photols_def
187
188
189   TYPE(species_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC     :: chem_species
190   TYPE(photols_def),ALLOCATABLE,DIMENSION(:),TARGET, PUBLIC     :: phot_frequen 
191
192   REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: spec_conc_1
193   REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: spec_conc_2
194   REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: spec_conc_3
195   REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: spec_conc_av       
196   REAL(kind=wp),ALLOCATABLE,DIMENSION(:,:,:,:),TARGET   :: freq_1
197
198   INTEGER,DIMENSION(nkppctrl)                           :: icntrl                            ! Fine tuning kpp
199   REAL(kind=wp),DIMENSION(nkppctrl)                     :: rcntrl                            ! Fine tuning kpp
200
201   CHARACTER(10), PUBLIC :: photolysis_scheme
202                                         ! 'constant',
203                                         ! 'simple' (Simple parameterisation from MCM, Saunders et al., 2003, Atmos. Chem. Phys., 3, 161-180
204                                         ! 'fastj'  (Wild et al., 2000, J. Atmos. Chem., 37, 245-282) STILL NOT IMPLEMENTED
205
206   PUBLIC nest_chemistry, nspec
207   PUBLIC nvar       
208   PUBLIC spc_names
209   PUBLIC spec_conc_2 
210
211!- Interface section
212  INTERFACE chem_boundary_conds
213      MODULE PROCEDURE chem_boundary_conds
214  END INTERFACE chem_boundary_conds
215
216   INTERFACE chem_init
217      MODULE PROCEDURE chem_init
218   END INTERFACE chem_init
219
220   INTERFACE chem_init_profiles
221      MODULE PROCEDURE chem_init_profiles
222   END INTERFACE chem_init_profiles
223
224   INTERFACE chem_parin
225      MODULE PROCEDURE chem_parin
226   END INTERFACE chem_parin
227
228   INTERFACE chem_integrate
229      MODULE PROCEDURE chem_integrate_ij
230   END INTERFACE chem_integrate
231
232   INTERFACE chem_swap_timelevel
233      MODULE PROCEDURE chem_swap_timelevel
234   END INTERFACE chem_swap_timelevel
235
236   INTERFACE chem_define_netcdf_grid
237      MODULE PROCEDURE chem_define_netcdf_grid
238   END INTERFACE chem_define_netcdf_grid
239
240   INTERFACE chem_data_output_3d
241      MODULE PROCEDURE chem_data_output_3d
242   END INTERFACE chem_data_output_3d
243
244   INTERFACE chem_check_data_output
245      MODULE PROCEDURE chem_check_data_output
246   END INTERFACE chem_check_data_output
247
248   INTERFACE chem_check_data_output_pr
249      MODULE PROCEDURE chem_check_data_output_pr
250   END INTERFACE chem_check_data_output_pr
251
252   INTERFACE chem_3d_data_averaging
253      MODULE PROCEDURE chem_3d_data_averaging 
254   END INTERFACE chem_3d_data_averaging
255
256   INTERFACE chem_wrd_local
257      MODULE PROCEDURE chem_wrd_local 
258   END INTERFACE chem_wrd_local
259
260   INTERFACE chem_rrd_local
261      MODULE PROCEDURE chem_rrd_local
262   END INTERFACE chem_rrd_local
263
264   INTERFACE chem_prognostic_equations
265      MODULE PROCEDURE chem_prognostic_equations
266      MODULE PROCEDURE chem_prognostic_equations_ij
267   END INTERFACE chem_prognostic_equations
268
269   INTERFACE chem_header
270      MODULE PROCEDURE chem_header
271   END INTERFACE chem_header
272
273   INTERFACE chem_emissions
274      MODULE PROCEDURE chem_emissions
275   END INTERFACE chem_emissions
276
277!    INTERFACE chem_wrd_global
278!       MODULE PROCEDURE chem_wrd_global
279!    END INTERFACE chem_wrd_global
280!
281!    INTERFACE chem_rrd_global
282!       MODULE PROCEDURE chem_rrd_global
283!    END INTERFACE chem_rrd_global
284
285
286   PUBLIC chem_3d_data_averaging, chem_boundary_conds, chem_check_data_output, &
287          chem_check_data_output_pr, chem_data_output_3d,                      &
288          chem_define_netcdf_grid, chem_emissions, chem_header, chem_init,     &
289          chem_init_profiles, chem_integrate, chem_wrd_local,                  &
290          chem_parin, chem_prognostic_equations,                               &
291          chem_rrd_local, chem_swap_timelevel
292         
293
294 CONTAINS
295
296!------------------------------------------------------------------------------!
297!
298! Description:
299! ------------
300!> Subroutine to initialize and set all boundary conditions for chemical species
301!------------------------------------------------------------------------------!
302
303 SUBROUTINE chem_boundary_conds( mode )                                           
304                                                                                 
305    USE control_parameters,                                                    & 
306        ONLY:  bc_radiation_l, bc_radiation_n, bc_radiation_r, bc_radiation_s
307
308    USE indices,                                                               & 
309        ONLY:  nxl, nxr, nzt
310                                                                                 
311    USE arrays_3d,                                                             &     
312        ONLY:  dzu                                               
313
314    USE surface_mod,                                                           &
315        ONLY:  bc_h                                                           
316
317    CHARACTER (len=*), INTENT(IN) :: mode
318    INTEGER(iwp) ::  i                                                            !< grid index x direction.
319    INTEGER(iwp) ::  j                                                            !< grid index y direction.
320    INTEGER(iwp) ::  k                                                            !< grid index z direction.
321    INTEGER(iwp) ::  kb                                                           !< variable to set respective boundary value, depends on facing.
322    INTEGER(iwp) ::  l                                                            !< running index boundary type, for up- and downward-facing walls.
323    INTEGER(iwp) ::  m                                                            !< running index surface elements.
324    INTEGER(iwp) ::  lsp                                                          !< running index for chem spcs.
325
326
327    SELECT CASE ( TRIM( mode ) )       
328       CASE ( 'init' )     
329          DO lsp = 1, nspec           
330            IF ( surface_csflux(lsp) == 9999999.9_wp )  THEN
331                 constant_csflux(lsp) = .FALSE.           
332            ENDIF
333          ENDDO
334
335          IF ( bc_cs_b == 'dirichlet' )    THEN
336             ibc_cs_b = 0 
337          ELSEIF ( bc_cs_b == 'neumann' )  THEN
338             ibc_cs_b = 1 
339          ELSE
340!             message_string = 'unknown boundary condition: bc_cs_b ="' // TRIM( bc_cs_b ) // '"'  ! bK commented
341             CALL message( 'chem_boundary_conds', 'CM0010', 1, 2, 0, 6, 0 )     !< chemistry_model_mod should have special error numbers --> "CHEM###",
342          ENDIF                                                                 
343!
344!--       Set Integer flags and check for possible erroneous settings for top
345!--       boundary condition. bK added *_cs_* here.
346          IF ( bc_cs_t == 'dirichlet' )             THEN
347             ibc_cs_t = 0 
348          ELSEIF ( bc_cs_t == 'neumann' )           THEN
349             ibc_cs_t = 1
350          ELSEIF ( bc_cs_t == 'initial_gradient' )  THEN
351             ibc_cs_t = 2
352          ELSEIF ( bc_cs_t == 'nested' )            THEN         
353             ibc_cs_t = 3
354          ELSE
355!            message_string = 'unknown boundary condition: bc_c_t ="' // TRIM( bc_cs_t ) // '"' 
356             CALL message( 'check_parameters', 'CM0011', 1, 2, 0, 6, 0 )
357          ENDIF
358
359     
360       CASE ( 'set_bc_bottomtop' )                   
361!--      Bottom boundary condtions for chemical species     
362          DO lsp = 1, nspec                                                     
363             IF ( ibc_cs_b == 0 ) THEN                   
364                DO l = 0, 1 
365!--             Set index kb: For upward-facing surfaces (l=0), kb=-1, i.e.
366!--             the chem_species(nsp)%conc_p value at the topography top (k-1)
367!--             is set; for downward-facing surfaces (l=1), kb=1, i.e. the
368!--             value at the topography bottom (k+1) is set.
369
370                   kb = MERGE( -1, 1, l == 0 )
371                   !$OMP PARALLEL DO PRIVATE( i, j, k )
372                   DO m = 1, bc_h(l)%ns
373                      i = bc_h(l)%i(m)           
374                      j = bc_h(l)%j(m)
375                      k = bc_h(l)%k(m)
376                      chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc(k+kb,j,i) 
377                   ENDDO                                       
378                ENDDO                                       
379         
380             ELSEIF ( ibc_cs_b == 1 ) THEN
381         ! in boundary_conds there is som extra loop over m here for passive tracer
382                DO l = 0, 1
383                   kb = MERGE( -1, 1, l == 0 )
384                   !$OMP PARALLEL DO PRIVATE( i, j, k )                                           
385                   DO m = 1, bc_h(l)%ns
386                      i = bc_h(l)%i(m)           
387                      j = bc_h(l)%j(m)
388                      k = bc_h(l)%k(m)
389                      chem_species(lsp)%conc_p(k+kb,j,i) = chem_species(lsp)%conc_p(k,j,i)
390
391!< @todo: chem_species(nsp)%conc_p(k+kb,j,i) = chem_species(nsp)%conc(k+kb,j,i),    &
392!<  pls loop over nsp=1, NSPEC.
393!< @todo: We should also think about the possibility to have &
394!< individual boundary conditions for each species? This means, bc_cs_b,       &
395!< bc_cs_t, ibc_cs_b, ibc_cs_t would need to be added to TYPE chem_species(nsp)%,   &
396!< and the loop over nsp would have to be put outside of this IF-clause.
397!< i think its better we keep the same bonundary cond i.e. dirichlet or neumann
398!< for all chem spcs. ... !bK
399
400                   ENDDO
401                ENDDO
402             ENDIF
403          ENDDO    ! end lsp loop 
404
405!--    Top boundary conditions for chemical species - Should this not be done for all species?
406          IF ( ibc_cs_t == 0 )  THEN                   
407             DO lsp = 1, nspec
408                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc(nzt+1,:,:)       
409             ENDDO
410          ELSEIF ( ibc_cs_t == 1 )  THEN
411             DO lsp = 1, nspec
412                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:)
413             ENDDO
414          ELSEIF ( ibc_cs_t == 2 )  THEN
415             DO lsp = 1, nspec
416                chem_species(lsp)%conc_p(nzt+1,:,:) = chem_species(lsp)%conc_p(nzt,:,:) + bc_cs_t_val(lsp) * dzu(nzt+1)
417             ENDDO
418                                    !@todo: bc_cs_t_val needs to be calculated,   
419                                    !see bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
420                                    !(in time_integration). pt_init would be the counterpart to
421                                    !chem_species(i)%conc_pr_init (see kchem_driver_FKa1408.f90 of my
422                                    !"Hints: prescribing initial concentration.
423          ENDIF
424!
425       CASE ( 'set_bc_lateral' )                       ! bK commented it
426!--          Lateral boundary conditions for chem species at inflow boundary
427!--          are automatically set when chem_species concentration is
428!--          initialized. The initially set value at the inflow boundary is not
429!--          touched during time integration, hence, this boundary value remains
430!--          at a constant value, which is the concentration that flows into the
431!--          domain.                                                           
432!--          Lateral boundary conditions for chem species at outflow boundary
433
434          IF ( bc_radiation_s )  THEN
435             DO lsp = 1, nspec
436                chem_species(lsp)%conc_p(:,nys-1,:) = chem_species(lsp)%conc_p(:,nys,:)
437             ENDDO
438          ELSEIF ( bc_radiation_n )  THEN
439             DO lsp = 1, nspec
440                chem_species(lsp)%conc_p(:,nyn+1,:) = chem_species(lsp)%conc_p(:,nyn,:)
441             ENDDO
442          ELSEIF ( bc_radiation_l )  THEN
443             DO lsp = 1, nspec
444                chem_species(lsp)%conc_p(:,:,nxl-1) = chem_species(lsp)%conc_p(:,:,nxl)
445             ENDDO
446          ELSEIF ( bc_radiation_r )  THEN
447             DO lsp = 1, nspec
448                chem_species(lsp)%conc_p(:,:,nxr+1) = chem_species(lsp)%conc_p(:,:,nxr)
449             ENDDO
450          ENDIF
451         
452    END SELECT
453
454 END SUBROUTINE chem_boundary_conds
455!
456!------------------------------------------------------------------------------!
457!
458! Description:
459! ------------
460!> Subroutine defining initial vertical profiles of chemical species (given by
461!> namelist parameters chem_profiles and chem_heights)  --> which should work
462!> analogue to parameters u_profile, v_profile and uv_heights)
463!------------------------------------------------------------------------------!
464 SUBROUTINE chem_init_profiles              !< SUBROUTINE is called from chem_init in case of
465                                            !< TRIM( initializing_actions ) /= 'read_restart_data'
466                                            !< We still need to see what has to be done in case of restart run
467    USE chem_modules
468
469    IMPLICIT NONE
470
471!-- Local variables
472    INTEGER ::  lsp        !< running index for number of species in derived data type species_def
473    INTEGER ::  lsp_pr     !< running index for number of species in cs_names, cs_profiles etc
474    INTEGER ::  lpr_lev    !< running index for profile level for each chem spcs.
475!    INTEGER ::  npr_lev    !< the next available profile lev
476
477!-----------------
478!-- To prescribe/initialize a vertically constant  'cs_profile', another parameter
479!-- "cs_surface" is introduced. If "cs_profile" and "cs_heights" are prescribed,
480!-- their values will override the constant profile given by "cs_surface".
481
482!     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
483       lsp_pr = 1
484       DO  WHILE ( TRIM( cs_name( lsp_pr ) ) /= 'novalue' )   !'novalue' is the default
485          DO  lsp = 1, nspec                                !
486!--          Create initial profile (conc_pr_init) for each chemical species
487             IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_pr) ) )  THEN   !
488!                IF ( cs_profile(1,1) == 9999999.9_wp ) THEN
489!--               Set a vertically constant profile based on the surface conc (cs_surface(lsp_pr)) of each species
490                  DO lpr_lev = 0, nzt+1
491                     chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_pr)
492                  ENDDO
493
494!                ELSE
495!                   IF ( cs_heights(lsp,1) /= 0.0_wp )  THEN
496!                      message_string = 'cs_heights(1,1) must be 0.0'
497!                      CALL message( 'chem_check_parameters', 'CM0012', 1, 2, 0, 6, 0 )
498!                   ENDIF
499
500!                   IF ( omega /= 0.0_wp )  THEN
501!                      message_string = 'Coriolis force must be switched off (by setting omega=0.0)' //  &
502!                                       ' when prescribing the forcing by u_profile and v_profile'
503!                      CALL message( 'check_parameters', 'PA0347', 1, 2, 0, 6, 0 )
504!                   ENDIF
505
506!                   use_prescribed_profile_data = .TRUE.
507!
508!                   npr_lev = 1
509! !                  chem_sddpecies(lsp)%conc_pr_init(0) = 0.0_wp
510!                   DO  lpr_lev = 1, nz+1
511!                      IF ( npr_lev < 100 )  THEN
512!                         DO  WHILE ( cs_heights(lsp, npr_lev+1) <= zu(lpr_lev) )
513!                            npr_lev = npr_lev + 1
514!                            IF ( npr_lev == 100 )  EXIT
515!                         ENDDO
516!                      ENDIF
517
518!                      IF ( npr_lev < 100  .AND.  cs_heights(lsp, npr_lev + 1) /= 9999999.9_wp )  THEN
519!                         chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp, npr_lev) +         &
520!                                                 ( zu(lpr_lev) - cs_heights(lsp, npr_lev) ) /         &
521!                                                 ( cs_heights(lsp, npr_lev + 1) - cs_heights(lsp, npr_lev ) ) * &
522!                                                 ( cs_profile(lsp, npr_lev + 1) - cs_profile(lsp, npr_lev ) )
523!                      ELSE
524!                         chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp, npr_lev)
525!                      ENDIF
526!                   ENDDO
527!                ENDIF
528!-- If a profile is prescribed explicity using cs_profiles and cs_heights, we then have to fill the
529!-- chem_species(lsp)%conc_pr_init for the specific "lsp" based on the cs_profiles(lsp_pr,:)
530!-- and cs_heights(lsp_pr,:).
531             ENDIF
532          ENDDO
533          lsp_pr = lsp_pr + 1
534       ENDDO
535!     ELSE
536!       IF (chem_debug1 ) print*,'code to be added for initializing_actions == read_restart_data'   !bK
537!     ENDIF
538
539!-- Now, go back to chem_init and use the contents of chem_species(lsp)%conc_pr_init to
540!-- initialize the 3D conc arrays, as it is currently taken care of in chem_set_constant_values.
541!-- After initializing the 3D arrays, these can be used to set the boundary conditions in the
542!-- subroutine kchem_boundary_conds, which should be called from boundary_conds.f90.
543
544
545 END SUBROUTINE chem_init_profiles
546!
547!------------------------------------------------------------------------------!
548!
549! Description:
550! ------------
551!> Subroutine initializating chemistry_model_mod
552!------------------------------------------------------------------------------!
553   SUBROUTINE chem_init
554
555      USE pegrid
556
557    IMPLICIT   none
558!-- local variables
559    INTEGER                 :: i,j               !< running index for for horiz numerical grid points
560    INTEGER                 :: lsp               !< running index for chem spcs
561    INTEGER                 :: lpr_lev           !< running index for chem spcs profile level
562!
563!-- NOPOINTER version not implemented yet
564! #if defined( __nopointer )
565!     message_string = 'The chemistry module only runs with POINTER version'
566!     CALL message( 'chemistry_model_mod', 'CM0001', 1, 2, 0, 6, 0 )     
567! #endif
568!
569!-- Allocate memory for chemical species
570    ALLOCATE( chem_species(nspec) )
571    ALLOCATE( spec_conc_1 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
572    ALLOCATE( spec_conc_2 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
573    ALLOCATE( spec_conc_3 (nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) )
574    ALLOCATE( spec_conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nspec) ) 
575    ALLOCATE( phot_frequen(nphot) ) 
576    ALLOCATE( freq_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg,nphot) )
577    ALLOCATE( bc_cs_t_val(nspec) )
578!
579!-- Initialize arrays
580    spec_conc_1 (:,:,:,:) = 0.0_wp
581    spec_conc_2 (:,:,:,:) = 0.0_wp
582    spec_conc_3 (:,:,:,:) = 0.0_wp
583    spec_conc_av(:,:,:,:) = 0.0_wp
584
585
586    DO lsp = 1, nspec
587       chem_species(lsp)%name    = spc_names(lsp)
588
589       chem_species(lsp)%conc   (nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_1 (:,:,:,lsp)
590       chem_species(lsp)%conc_p (nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_2 (:,:,:,lsp)
591       chem_species(lsp)%tconc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_3 (:,:,:,lsp)
592       chem_species(lsp)%conc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg)       => spec_conc_av(:,:,:,lsp)     
593
594       ALLOCATE (chem_species(lsp)%cssws_av(nysg:nyng,nxlg:nxrg))                   
595       chem_species(lsp)%cssws_av    = 0.0_wp
596!
597!--    (todo (FK): This needs to be revised. This block must go somewhere else)                                         
598!      IF ( ws_scheme_sca )  THEN                                               
599          ALLOCATE (chem_species(lsp)%flux_s_cs(nzb+1:nzt,0:threads_per_task-1))   
600          ALLOCATE (chem_species(lsp)%diss_s_cs(nzb+1:nzt,0:threads_per_task-1))   
601          ALLOCATE (chem_species(lsp)%flux_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1)) 
602          ALLOCATE (chem_species(lsp)%diss_l_cs(nzb+1:nzt,nys:nyn,0:threads_per_task-1))   
603          chem_species(lsp)%flux_s_cs = 0.0_wp                                     
604          chem_species(lsp)%flux_l_cs = 0.0_wp                                     
605          chem_species(lsp)%diss_s_cs = 0.0_wp                                     
606          chem_species(lsp)%diss_l_cs = 0.0_wp                                     
607!      ENDIF         
608!
609!--   Allocate memory for initial concentration profiles
610!--   (concentration values come from namelist)
611!--   (todo (FK): Because of this, chem_init is called in palm before
612!--               check_parameters, since conc_pr_init is used there.
613!--               We have to find another solution since chem_init should
614!--               eventually be called from init_3d_model!!)
615      ALLOCATE ( chem_species(lsp)%conc_pr_init(0:nz+1) )
616      chem_species(lsp)%conc_pr_init(:) = 0.0_wp
617
618    ENDDO
619
620!
621!-- Set initial concentration of profiles prescribed by parameters cs_profile
622!-- and cs_heights in the namelist &chemistry_parameters
623!-- (todo (FK): chem_init_profiles not ready yet, has some bugs)
624!     CALL chem_init_profiles
625!
626!-- Initialize model variables
627
628
629    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
630         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
631
632
633!--    First model run of a possible job queue.
634!--    Initial profiles of the variables must be computes.
635       IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
636!            CALL location_message( 'initializing with 1D model profiles', .FALSE. )
637!
638!          CALL init_1d_model    ...... decide to call or not later     !bK
639
640!--        Transfer initial profiles to the arrays of the 3D model
641          DO lsp = 1, nspec
642             DO  i = nxlg, nxrg
643                DO  j = nysg, nyng
644                   DO lpr_lev = 1, nz + 1
645                      chem_species(lsp)%conc(lpr_lev,j,i) = chem_species(lsp)%conc_pr_init(lpr_lev)
646                   ENDDO
647                ENDDO
648             ENDDO   
649          ENDDO
650       
651       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 )    &
652       THEN
653!           CALL location_message( 'initializing with constant profiles', .FALSE. )
654
655
656
657!--       Set initial horizontal velocities at the lowest computational grid
658!--       levels to zero in order to avoid too small time steps caused by the
659!--       diffusion limit in the initial phase of a run (at k=1, dz/2 occurs
660!--       in the limiting formula!).
661
662          DO lsp = 1, nspec 
663             DO  i = nxlg, nxrg
664                DO  j = nysg, nyng
665                   chem_species(lsp)%conc(:,j,i) = chem_species(lsp)%conc_pr_init    !ITS THERE bK
666                ENDDO
667             ENDDO
668          ENDDO
669
670!        ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 )                  &
671!        THEN
672!           CALL location_message( 'initializing by user', .FALSE. )
673!
674!--       Initialization will completely be done by the user
675!--       (FK: This should be called only once, in init_3d_model, i.e. remove it here)
676!           CALL user_init_3d_model
677!           CALL location_message( 'finished', .TRUE. )
678
679       ENDIF
680
681!-- Store initial chem spcs profile
682!         DO lsp = 1, nvar
683!          hom_cs(:,1,115,:) = SPREAD(  chem_species(lsp)%conc(:,nys,nxl), 2, statistic_regions+1 )
684!         ENDDO
685!
686!-- If required, change the surface chem spcs at the start of the 3D run
687       IF ( cs_surface_initial_change(1) /= 0.0_wp ) THEN           
688          DO lsp = 1, nspec 
689             chem_species(lsp)%conc(nzb,:,:) = chem_species(lsp)%conc(nzb,:,:) +  &
690                                               cs_surface_initial_change(lsp)
691          ENDDO
692       ENDIF 
693!
694!-- Initiale old and new time levels.
695       DO lsp = 1, nvar
696          chem_species(lsp)%tconc_m = 0.0_wp                     
697          chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
698       ENDDO
699
700    ENDIF
701
702
703
704!--- new code add above this line
705    DO lsp = 1, nphot
706       phot_frequen(lsp)%name = phot_names(lsp)
707!        IF( myid == 0 )  THEN
708!           WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,trim(phot_names(lsp))
709!        ENDIF
710         phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => freq_1(:,:,:,lsp)
711    ENDDO
712
713!--   Set initial values
714!    Not required any more, this can now be done with the namelist  by setting cs_surface
715!    and cs_name without specifying cs_profile (Nevertheless, we still want to keep it for a while)
716!   IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN             
717!         CALL set_const_initial_values
718!   ENDIF
719
720    RETURN
721
722    CONTAINS
723!------------------------------------------------------------------------------!
724!
725! Description:
726! ------------
727!> Subroutine setting constant initial values of chemical species
728!------------------------------------------------------------------------------!
729     SUBROUTINE set_const_initial_values
730!    Not required any more, this can now be done with the namelist  by setting cs_surface
731!    and cs_name without specifying cs_profile (Nevertheless, we still want to keep it for a while)
732         IMPLICIT   none
733
734!--      local variables
735         INTEGER                  :: lsp
736
737         IF(myid == 0)  THEN
738            write(6,'(/,a,/)')  ' chemics >>>> Set constant Initial Values: '
739         ENDIF
740
741!        Default values are taken from smog.def from supplied kpp example
742         DO lsp = 1, nspec
743            IF(trim(chem_species(lsp)%name) == 'NO')           THEN
744!              chem_species(lsp)%conc = 8.725*1.0E+08
745!              chem_species(lsp)%conc = 0.05_wp                          !added by bK
746               chem_species(lsp)%conc = 0.01_wp                          !added by RFo
747            ELSEIF (trim(chem_species(lsp)%name) == 'NO2')     THEN
748!              chem_species(lsp)%conc = 2.240*1.0E+08             
749!              chem_species(lsp)%conc = 0.01_wp                          !added by bK
750               chem_species(lsp)%conc = 0.05_wp                          !added by RFo
751            ELSEIF( trim( chem_species(lsp)%name ) == 'O3' )   THEN
752               chem_species(lsp)%conc = 0.05_wp                          !added by bK
753            ELSEIF(trim(chem_species(lsp)%name) == 'H2O')      THEN
754!              chem_species(lsp)%conc = 5.326*1.0E+11
755               chem_species(lsp)%conc = 1.30*1.0E+4_wp                   !added by bK
756            ELSEIF(trim(chem_species(lsp)%name) == 'O2')       THEN
757               chem_species(lsp)%conc = 2.0*1.0E+5_wp                    !added by bK
758            ELSEIF(trim(chem_species(lsp)%name) == 'RH')       THEN
759               chem_species(lsp)%conc = 0.001_wp                         !added by RFo
760            ELSEIF(trim(chem_species(lsp)%name) == 'CO')       THEN
761               chem_species(lsp)%conc = 0.5_wp                           !added by RFo
762            ELSEIF(trim(chem_species(lsp)%name) == 'RCHO')     THEN
763!              chem_species(lsp)%conc = 2.0_wp                           !added by bK
764               chem_species(lsp)%conc = 0.01_wp                          !added by RFo
765!           ELSEIF(trim(chem_species(lsp)%name) == 'OH')       THEN
766!              chem_species(lsp)%conc = 1.0*1.0E-07_wp                   !added by bK
767!           ELSEIF(trim(chem_species(lsp)%name) == 'HO2')      THEN
768!              chem_species(lsp)%conc = 1*1.0E-7_wp                      !added by bK
769!           ELSEIF(trim(chem_species(lsp)%name) == 'RCOO2')    THEN      ! corrected RFo
770!              chem_species(lsp)%conc = 1.0*1.0E-7_wp                    !added by bK
771!           ELSEIF(trim(chem_species(lsp)%name) == 'RCOO2NO2') THEN
772!              chem_species(lsp)%conc = 1.0*1.0E-7_wp                    !added by bK
773            ELSE
774!              H2O = 2.0e+04;
775               chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) = 0.0_wp
776            ENDIF
777
778            IF(myid == 0)  THEN
779               WRITE(6,'(a,3x,a,3x,a,e12.4)')  '   Species:     ',chem_species(lsp)%name(1:7), & 
780                                        'Initial Value = ',chem_species(lsp)%conc(nzb,nysg,nxlg)
781            ENDIF
782         ENDDO
783
784!   #if defined( __nopointer )
785!   !kk      Hier mit message abbrechen
786!            if(myid == 0)  then
787!               write(6,*)  '   KPP does only run with POINTER Version'
788!            end if
789!            stop 'error'
790!   #endif
791
792         RETURN
793      END SUBROUTINE set_const_initial_values
794
795
796   END SUBROUTINE chem_init
797!
798!------------------------------------------------------------------------------!
799!
800! Description:
801! ------------
802!> Subroutine defining parin for &chemistry_parameters for chemistry model
803!------------------------------------------------------------------------------!
804   SUBROUTINE chem_parin
805   
806      USE control_parameters,                                                  &
807          ONLY: air_chemistry
808         
809      USE chem_modules
810     
811      USE kinds
812
813      IMPLICIT   none
814
815      CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
816     
817      REAL(wp), DIMENSION(nmaxfixsteps) ::   my_steps   !< List of fixed timesteps   my_step(1) = 0.0 automatic stepping
818
819      NAMELIST /chemistry_parameters/   bc_cs_b,                          &
820                                        bc_cs_t,                          &
821                                        call_chem_at_all_substeps,        &
822                                        chem_debug0,                      &
823                                        chem_debug1,                      &
824                                        chem_debug2,                      &
825                                        chem_gasphase_on,                 &
826                                        cs_heights,                       &
827                                        cs_name,                          &
828                                        cs_profile,                       &
829                                        cs_profile_name,                  & 
830                                        cs_surface,                       &
831                                        emiss_factor_main,                &
832                                        emiss_factor_side,                &                     
833                                        icntrl,                           &
834                                        main_street_id,                   &
835                                        max_street_id,                    &
836                                        my_steps,                         &
837                                        nest_chemistry,                   &
838                                        rcntrl,                           &
839                                        side_street_id,                   &
840                                        photolysis_scheme,                &
841                                        wall_csflux,                      &
842                                        cs_vertical_gradient,             &
843                                        top_csflux,                       & 
844                                        surface_csflux,                   &
845                                        surface_csflux_name,              &
846                                        cs_surface_initial_change,        &
847                                        cs_vertical_gradient_level
848                             
849!-- analogue to chem_names(nspj) we could invent chem_surfaceflux(nspj) and chem_topflux(nspj)
850!-- so this way we could prescribe a specific flux value for each species
851!>  chemistry_parameters for initial profiles
852!>  cs_names = 'O3', 'NO2', 'NO', ...   to set initial profiles)
853!>  cs_heights(1,:) = 0.0, 100.0, 500.0, 2000.0, .... (height levels where concs will be prescribed for O3)
854!>  cs_heights(2,:) = 0.0, 200.0, 400.0, 1000.0, .... (same for NO2 etc.)
855!>  cs_profiles(1,:) = 10.0, 20.0, 20.0, 30.0, .....  (chem spcs conc at height lvls chem_heights(1,:)) etc.
856!>  If the respective concentration profile should be constant with height, then use "cs_surface( number of spcs)"
857!>  then write these cs_surface values to chem_species(lsp)%conc_pr_init(:)
858
859!
860!--   Read chem namelist   
861!--   (todo: initialize these parameters in declaration part, do this for
862!--          all chemistry_parameters namelist parameters)
863      icntrl    = 0
864      rcntrl    = 0.0_wp
865      my_steps  = 0.0_wp
866      icntrl(2) = 1                                 
867      photolysis_scheme = 'simple'
868      atol = 1.0_wp
869      rtol = 0.01_wp
870!
871!--   Try to find chemistry package
872      REWIND ( 11 )
873      line = ' '
874      DO WHILE ( INDEX( line, '&chemistry_parameters' ) == 0 )
875         READ ( 11, '(A)', END=20 )  line
876      ENDDO
877      BACKSPACE ( 11 )
878!
879!--   Read chemistry namelist
880      READ ( 11, chemistry_parameters, ERR = 10, END = 20 )
881!
882!--   Enable chemistry model
883      air_chemistry = .TRUE.
884      GOTO 20
885
886 10   BACKSPACE( 11 )
887      READ( 11 , '(A)') line
888      CALL parin_fail_message( 'chemistry_parameters', line )
889
890 20   CONTINUE
891
892      t_steps = my_steps          !(todo: Why not directly make t_steps a
893                                  !       namelist parameter?)
894
895   END SUBROUTINE chem_parin
896
897 !
898!------------------------------------------------------------------------------!
899!
900! Description:
901! ------------
902!> Subroutine to integrate chemical species in the given chemical mechanism
903!------------------------------------------------------------------------------!
904
905   SUBROUTINE chem_integrate_ij (i, j)
906
907      USE cpulog,                                                              &
908           ONLY: cpu_log, log_point
909      USE statistics,                                                          &   ! ## RFo
910           ONLY:  weight_pres
911       USE control_parameters,                                                 &   ! ## RFo
912           ONLY:  dt_3d, intermediate_timestep_count
913
914      IMPLICIT   none
915      INTEGER,INTENT(IN)       :: i,j
916
917!--   local variables
918      INTEGER(iwp) ::  lsp                                                          !< running index for chem spcs.
919      INTEGER(iwp) ::  lph                                                          !< running index for photolysis frequencies
920      INTEGER                  :: istatf
921      INTEGER,dimension(20)    :: istatus
922      REAL(kind=wp),dimension(nzb+1:nzt,nspec)                :: tmp_conc           
923      REAL(kind=wp),dimension(nzb+1:nzt)                      :: tmp_temp
924      REAL(kind=wp),dimension(nzb+1:nzt,nphot)                :: tmp_phot
925      REAL(kind=wp),dimension(nzb+1:nzt)                      :: tmp_fact   
926     REAL(kind=wp),dimension(nzb+1:nzt)                      :: tmp_fact_i    !< conversion factor between molecules cm^{-3} and ppm
927     REAL(wp)                         ::  conv                                !< conversion factor
928     REAL(wp), PARAMETER              ::  ppm2fr  = 1.0e-6_wp                 !< Conversion factor ppm to fraction
929     REAL(wp), PARAMETER              ::  t_std   = 273.15_wp                 !< standard pressure (Pa)
930     REAL(wp), PARAMETER              ::  p_std   = 101325.0_wp               !< standard pressure (Pa)
931     REAL(wp), PARAMETER              ::  vmolcm  = 22.414e3_wp               !< Mole volume (22.414 l) in cm^{-3}
932     REAL(wp), PARAMETER              ::  xna     = 6.022e23_wp               !< Avogadro number (molecules/mol)
933
934
935      REAL(kind=wp)  :: dt_chem                                             
936
937       CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'start' )
938!<     Set chem_gasphase_on to .FALSE. if you want to skip computation of gas phase chemistry
939    IF (chem_gasphase_on) THEN
940
941       tmp_temp(:) = pt(:,j,i) * ( hyp(nzb+1:nzt) / 100000.0_wp )**0.286_wp
942! ppm to molecules/cm**3
943!      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 )
944       conv = ppm2fr * xna / vmolcm
945       tmp_fact(:) = conv * t_std * hyp(nzb+1:nzt) / (tmp_temp(:) * p_std)
946       tmp_fact_i = 1.0_wp/tmp_fact
947
948       CALL fill_temp (istatf, tmp_temp)                                      !< Load constant temperature into kpp context
949!      CALL fill_temp (istatf, pt(nzb+1:nzt,j,i))                             !< Load temperature into kpp context
950
951       DO lsp = 1,nspec
952          tmp_conc(:,lsp) = chem_species(lsp)%conc(nzb+1:nzt,j,i) * tmp_fact(:) 
953       ENDDO
954
955       DO lph = 1,nphot
956          tmp_phot(:,lph) = phot_frequen(lph)%freq(nzb+1:nzt,j,i)               
957       ENDDO
958
959       IF(myid == 0 .AND. chem_debug0 ) THEN
960          IF (i == 10 .and. j == 10) WRITE(0,*) 'begin chemics step ',dt_3d
961       ENDIF
962
963!--    Compute length of time step     # RFo
964       IF ( call_chem_at_all_substeps )  THEN
965          dt_chem = dt_3d * weight_pres(intermediate_timestep_count)
966       ELSE
967          dt_chem = dt_3d
968       ENDIF
969
970       CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'start' )
971
972
973       CALL chem_gasphase_integrate (dt_chem, tmp_conc, tmp_temp, tmp_phot, istatus=istatus)
974
975
976       CALL cpu_log( log_point(81), '{chem_gasphase_integrate}', 'stop' )
977
978       DO lsp = 1,nspec
979          chem_species(lsp)%conc (nzb+1:nzt,j,i) = tmp_conc(:,lsp) * tmp_fact_i(:)  ! RFo
980       ENDDO
981
982!      IF (myid == 0 .AND. chem_debug2 ) THEN
983!         IF (i == 10 .and. j == 10)   WRITE(6,'(a,8i7)') ' KPP Status ',istatus(1:8)
984!         write(6,'(a,8i7)') ' KPP Status ',istatus(1:8)
985!      ENDIF
986
987    ENDIF
988       CALL cpu_log( log_point(80), '[chem_integrate_ij]', 'stop' )
989
990       RETURN
991   END SUBROUTINE chem_integrate_ij
992!
993!------------------------------------------------------------------------------!
994!
995! Description:
996! ------------
997!> Subroutine for swapping of timelevels for chemical species
998!> called out from subroutine swap_timelevel
999!------------------------------------------------------------------------------!
1000
1001   SUBROUTINE chem_swap_timelevel (level)
1002      IMPLICIT   none
1003
1004      INTEGER,INTENT(IN)                  :: level
1005
1006!--   local variables
1007
1008      INTEGER               :: lsp
1009
1010!        print*,' *** entering chem_swap_timelevel ***) '
1011      if(level == 0)  then
1012         do lsp=1, nvar                                       
1013            chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_1(:,:,:,lsp)
1014            chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_2(:,:,:,lsp)
1015         end do
1016      else
1017         do lsp=1, nvar                                       
1018            chem_species(lsp)%conc(nzb:nzt+1,nysg:nyng,nxlg:nxrg)    => spec_conc_2(:,:,:,lsp)
1019            chem_species(lsp)%conc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  => spec_conc_1(:,:,:,lsp)
1020         end do
1021      end if
1022
1023      RETURN
1024   END SUBROUTINE chem_swap_timelevel
1025
1026!------------------------------------------------------------------------------!
1027!
1028! Description:
1029! ------------
1030!> Subroutine defining appropriate grid for netcdf variables.
1031!> It is called out from subroutine netcdf.
1032!------------------------------------------------------------------------------!
1033   SUBROUTINE chem_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
1034
1035      IMPLICIT NONE
1036
1037      CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
1038      LOGICAL, INTENT(OUT)           ::  found       !<
1039      CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
1040      CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
1041      CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
1042
1043      found  = .TRUE.
1044
1045      if(var(1:3) == 'kc_')   then                    !< always the same grid for chemistry variables
1046            grid_x = 'x'
1047            grid_y = 'y'
1048            grid_z = 'zu'                             !< kk Use same z axis as u variables. Has to be checked if OK
1049      else
1050            found  = .FALSE.
1051            grid_x = 'none'
1052            grid_y = 'none'
1053            grid_z = 'none'
1054      end if
1055
1056!     write(6,*) 'chem_define_netcdf_grid ',TRIM(var),' ',trim(grid_x),' ',found
1057
1058   END SUBROUTINE chem_define_netcdf_grid
1059!
1060!------------------------------------------------------------------------------!
1061!
1062! Description:
1063! ------------
1064!> Subroutine for checking data output for chemical species
1065!------------------------------------------------------------------------------!
1066
1067   SUBROUTINE chem_check_data_output( var, unit, i, ilen, k )
1068
1069      IMPLICIT NONE
1070
1071      CHARACTER (LEN=*) ::  unit     !<
1072      CHARACTER (LEN=*) ::  var      !<
1073
1074      INTEGER(iwp) ::  i, ilen, k, lsp
1075
1076      CHARACTER(len=16)    ::  spec_name
1077
1078      unit = 'illegal'
1079
1080      spec_name = TRIM(var(4:))             !< var 1:3 is 'kc_'.
1081
1082       DO lsp=1,nspec
1083         IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
1084            unit = 'ppm'
1085         ENDIF
1086         ! It is possible  to plant PM10 and PM25 into the gasphase chemistry code
1087!        ! as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
1088!        ! set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
1089         IF (spec_name(1:2) == 'PM')   THEN 
1090            unit = 'ug m-3'
1091         ENDIF
1092       ENDDO
1093
1094       DO lsp=1,nphot                                                     
1095         IF (TRIM(spec_name) == TRIM(phot_frequen(lsp)%name))   THEN
1096            unit = 'sec-1'
1097         ENDIF
1098       ENDDO
1099
1100
1101       RETURN
1102   END SUBROUTINE chem_check_data_output
1103!
1104!------------------------------------------------------------------------------!
1105!
1106! Description:
1107! ------------
1108!> Subroutine for checking data output of profiles for chemistry model
1109!------------------------------------------------------------------------------!
1110
1111   SUBROUTINE chem_check_data_output_pr( variable, var_count, unit, dopr_unit )
1112
1113
1114    USE arrays_3d,                                                             &
1115        ONLY: zu
1116
1117    USE control_parameters,                                                    &
1118        ONLY: data_output_pr, message_string, air_chemistry
1119
1120    USE indices
1121
1122    USE profil_parameter
1123
1124    USE statistics
1125
1126
1127    IMPLICIT NONE
1128
1129    CHARACTER (LEN=*) ::  unit     !<
1130    CHARACTER (LEN=*) ::  variable !<
1131    CHARACTER (LEN=*) ::  dopr_unit
1132    CHARACTER(len=16) ::  spec_name
1133 
1134    INTEGER(iwp) ::  var_count, lsp  !<
1135   
1136
1137    spec_name = TRIM(variable(4:))   
1138!             write(9,*) 'fm #32 .. air_chemistry ', air_chemistry
1139
1140          IF (  .NOT.  air_chemistry )  THEN
1141                 message_string = 'data_output_pr = ' //                        &
1142                 TRIM( data_output_pr(var_count) ) // ' is not imp' // &
1143                          'lemented for air_chemistry = .FALSE.'
1144!                CALL message( 'check_parameters', 'PA0XXX', 1, 2, 0, 6, 0 )
1145          ELSE
1146             DO lsp = 1, nspec
1147                IF (TRIM( spec_name ) == TRIM( chem_species(lsp)%name ) ) THEN
1148                    dopr_index(var_count) = 900 
1149                    dopr_unit  = 'ppm'
1150                    hom(:,2,900,:) = SPREAD( zu, 2, statistic_regions+1 )
1151                ENDIF
1152             ENDDO
1153          ENDIF
1154
1155   END SUBROUTINE chem_check_data_output_pr
1156!
1157!------------------------------------------------------------------------------!
1158!
1159! Description:
1160! ------------
1161!> Subroutine defining 3D output variables for chemical species
1162!------------------------------------------------------------------------------!
1163
1164
1165   SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
1166
1167
1168      USE indices
1169
1170      USE kinds
1171
1172
1173      IMPLICIT NONE
1174
1175      CHARACTER (LEN=*) ::  variable !<
1176
1177      INTEGER(iwp) ::  av    !<
1178      INTEGER(iwp) ::  nzb_do !< lower limit of the data output (usually 0)
1179      INTEGER(iwp) ::  nzt_do !< vertical upper limit of the data output (usually nz_do3d)
1180
1181      LOGICAL      ::  found !<
1182
1183      REAL(wp) ::  fill_value !<
1184      REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf 
1185
1186
1187      !-- local variables
1188
1189      INTEGER              ::  i, j, k, lsp
1190      CHARACTER(len=16)    ::  spec_name
1191
1192
1193      found = .FALSE.
1194
1195      spec_name = TRIM(variable(4:))
1196!av == 0
1197
1198       DO lsp=1,nspec
1199          IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
1200             IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM(variable)  //  &
1201                                                          TRIM(chem_species(lsp)%name)       
1202             
1203             IF (av == 0) THEN
1204                DO  i = nxl, nxr
1205                   DO  j = nys, nyn
1206                      DO  k = nzb_do, nzt_do
1207                          local_pf(i,j,k) = MERGE(                             &
1208                                              chem_species(lsp)%conc(k,j,i),   &
1209                                              REAL( fill_value, KIND = wp ),   &
1210                                              BTEST( wall_flags_0(k,j,i), 0 ) )
1211                      ENDDO
1212                   ENDDO
1213                ENDDO
1214           
1215             ELSE
1216                DO  i = nxl, nxr
1217                   DO  j = nys, nyn
1218                      DO  k = nzb_do, nzt_do
1219                          local_pf(i,j,k) = MERGE(                             &
1220                                              chem_species(lsp)%conc_av(k,j,i),&
1221                                              REAL( fill_value, KIND = wp ),   &
1222                                              BTEST( wall_flags_0(k,j,i), 0 ) )
1223                      ENDDO
1224                   ENDDO
1225                ENDDO
1226             ENDIF
1227
1228            found = .TRUE.
1229          ENDIF
1230       ENDDO
1231
1232       RETURN
1233   END SUBROUTINE chem_data_output_3d
1234!
1235!------------------------------------------------------------------------------!
1236!
1237! Description:
1238! ------------
1239!> Subroutine for averaging 3D data of chemical species. Due to the fact that
1240!> the averaged chem arrays are allocated in chem_init, no if-query concerning
1241!> the allocation is required (in any mode). Attention: If you just specify an
1242!> averaged output quantity in the _p3dr file during restarts the first output
1243!> includes the time between the beginning of the restart run and the first
1244!> output time (not necessarily the whole averaging_interval you have
1245!> specified in your _p3d/_p3dr file )
1246!------------------------------------------------------------------------------!
1247
1248    SUBROUTINE chem_3d_data_averaging ( mode, variable )
1249
1250    USE control_parameters
1251
1252    USE indices
1253
1254    USE kinds
1255
1256    USE surface_mod,                                                         &
1257        ONLY: surf_def_h, surf_lsm_h, surf_usm_h   
1258 
1259    IMPLICIT NONE
1260 
1261    CHARACTER (LEN=*) ::  mode    !<
1262    CHARACTER (LEN=*) :: variable !<
1263 
1264    LOGICAL      ::  match_def !< flag indicating natural-type surface
1265    LOGICAL      ::  match_lsm !< flag indicating natural-type surface
1266    LOGICAL      ::  match_usm !< flag indicating urban-type surface
1267
1268    INTEGER(iwp) ::  i                  !< grid index x direction
1269    INTEGER(iwp) ::  j                  !< grid index y direction
1270    INTEGER(iwp) ::  k                  !< grid index z direction
1271    INTEGER(iwp) ::  m                  !< running index surface type
1272    INTEGER(iwp) :: lsp                 !< running index for chem spcs
1273 
1274    IF ( mode == 'allocate' )  THEN
1275       DO lsp = 1, nspec
1276          IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
1277             chem_species(lsp)%conc_av = 0.0_wp
1278          ENDIF
1279       ENDDO
1280
1281    ELSEIF ( mode == 'sum' )  THEN
1282   
1283       DO lsp = 1, nspec
1284          IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
1285             DO  i = nxlg, nxrg
1286                DO  j = nysg, nyng
1287                   DO  k = nzb, nzt+1
1288                       chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) + &
1289                                                          chem_species(lsp)%conc(k,j,i)
1290                   ENDDO
1291                ENDDO
1292             ENDDO
1293          ELSEIF ( TRIM(variable(4:)) == TRIM('cssws*') )        THEN
1294             DO  i = nxl, nxr
1295                DO  j = nys, nyn
1296                   match_def = surf_def_h(0)%start_index(j,i) <=               &
1297                               surf_def_h(0)%end_index(j,i)
1298                   match_lsm = surf_lsm_h%start_index(j,i) <=                  &
1299                               surf_lsm_h%end_index(j,i)
1300                   match_usm = surf_usm_h%start_index(j,i) <=                  &
1301                               surf_usm_h%end_index(j,i)
1302
1303                   IF ( match_def )  THEN
1304                      m = surf_def_h(0)%end_index(j,i)
1305                      chem_species(lsp)%cssws_av(j,i) =                        &
1306                                             chem_species(lsp)%cssws_av(j,i) + &
1307                                             surf_def_h(0)%cssws(lsp,m)
1308                   ELSEIF ( match_lsm  .AND.  .NOT. match_usm )  THEN
1309                      m = surf_lsm_h%end_index(j,i)
1310                      chem_species(lsp)%cssws_av(j,i) =                        &
1311                                             chem_species(lsp)%cssws_av(j,i) + &
1312                                             surf_lsm_h%cssws(lsp,m)
1313                   ELSEIF ( match_usm )  THEN
1314                      m = surf_usm_h%end_index(j,i)
1315                      chem_species(lsp)%cssws_av(j,i) =                        &
1316                                             chem_species(lsp)%cssws_av(j,i) + &
1317                                             surf_usm_h%cssws(lsp,m)
1318                   ENDIF
1319                ENDDO
1320             ENDDO
1321          ENDIF
1322       ENDDO
1323 
1324    ELSEIF ( mode == 'average' )  THEN
1325       DO lsp = 1, nspec
1326          IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
1327             DO  i = nxlg, nxrg
1328                DO  j = nysg, nyng
1329                   DO  k = nzb, nzt+1
1330                       chem_species(lsp)%conc_av(k,j,i) = chem_species(lsp)%conc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
1331                   ENDDO
1332                ENDDO
1333             ENDDO
1334
1335          ELSEIF (TRIM(variable(4:)) == TRIM('cssws*') )            THEN
1336             DO i = nxlg, nxrg
1337                DO  j = nysg, nyng
1338                     chem_species(lsp)%cssws_av(j,i) = chem_species(lsp)%cssws_av(j,i) / REAL( average_count_3d, KIND=wp )
1339                ENDDO
1340             ENDDO
1341                CALL exchange_horiz_2d( chem_species(lsp)%cssws_av, nbgp )                 
1342          ENDIF
1343       ENDDO
1344       
1345    ENDIF     
1346
1347    END SUBROUTINE chem_3d_data_averaging
1348
1349!------------------------------------------------------------------------------!
1350!
1351! Description:
1352! ------------
1353!> Subroutine to write restart data for chemistry model
1354!------------------------------------------------------------------------------!
1355 SUBROUTINE chem_wrd_local
1356
1357
1358    IMPLICIT NONE
1359
1360    INTEGER(iwp) ::  lsp !<
1361
1362!      REAL(kind=wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg)   :: chems_conc
1363
1364
1365    DO  lsp = 1, nspec
1366
1367       CALL wrd_write_string( TRIM( chem_species(lsp)%name ))
1368       WRITE ( 14 )  chem_species(lsp)%conc
1369
1370       CALL wrd_write_string( TRIM( chem_species(lsp)%name )//'_av' )
1371       WRITE ( 14 )  chem_species(lsp)%conc_av
1372
1373    ENDDO
1374
1375
1376 END SUBROUTINE chem_wrd_local
1377
1378!------------------------------------------------------------------------------!
1379!
1380! Description:
1381! ------------
1382!> Subroutine to read restart data of chemical species
1383!------------------------------------------------------------------------------!
1384
1385 SUBROUTINE chem_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,         &
1386                            nxr_on_file, nynf, nync, nyn_on_file, nysf, nysc,  &
1387                            nys_on_file, tmp_3d, found )   
1388                                     
1389    USE control_parameters
1390           
1391    USE indices
1392   
1393    USE pegrid
1394
1395    IMPLICIT NONE
1396
1397    CHARACTER (LEN=20) :: spc_name_av !<   
1398       
1399    INTEGER(iwp) ::  i, lsp          !<
1400    INTEGER(iwp) ::  k               !<
1401    INTEGER(iwp) ::  nxlc            !<
1402    INTEGER(iwp) ::  nxlf            !<
1403    INTEGER(iwp) ::  nxl_on_file     !<   
1404    INTEGER(iwp) ::  nxrc            !<
1405    INTEGER(iwp) ::  nxrf            !<
1406    INTEGER(iwp) ::  nxr_on_file     !<   
1407    INTEGER(iwp) ::  nync            !<
1408    INTEGER(iwp) ::  nynf            !<
1409    INTEGER(iwp) ::  nyn_on_file     !<   
1410    INTEGER(iwp) ::  nysc            !<
1411    INTEGER(iwp) ::  nysf            !<
1412    INTEGER(iwp) ::  nys_on_file     !<   
1413   
1414    LOGICAL, INTENT(OUT)  :: found 
1415
1416    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
1417
1418
1419    found = .FALSE. 
1420
1421
1422       IF ( ALLOCATED(chem_species) )  THEN
1423
1424          DO  lsp = 1, nspec
1425
1426              !< for time-averaged chemical conc.
1427             spc_name_av  =  TRIM(chem_species(lsp)%name)//'_av'
1428
1429             IF (restart_string(1:length) == TRIM(chem_species(lsp)%name) )    &
1430             THEN
1431                !< read data into tmp_3d
1432                IF ( k == 1 )  READ ( 13 )  tmp_3d 
1433                !< fill ..%conc in the restart run   
1434                chem_species(lsp)%conc(:,nysc-nbgp:nync+nbgp,                  &
1435                                       nxlc-nbgp:nxrc+nbgp) =                  & 
1436                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1437                found = .TRUE.
1438             ELSEIF (restart_string(1:length) == spc_name_av )  THEN
1439                IF ( k == 1 )  READ ( 13 )  tmp_3d
1440                chem_species(lsp)%conc_av(:,nysc-nbgp:nync+nbgp,               &
1441                                          nxlc-nbgp:nxrc+nbgp) =               &
1442                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1443                found = .TRUE.
1444             ENDIF
1445
1446          ENDDO
1447
1448       ENDIF
1449
1450
1451 END SUBROUTINE chem_rrd_local
1452
1453
1454!------------------------------------------------------------------------------!
1455!
1456! Description:
1457! ------------
1458!> Subroutine calculating prognostic equations for chemical species
1459!> (cache-optimized).
1460!> Routine is called separately for each chemical species over a loop from
1461!> prognostic_equations.
1462!------------------------------------------------------------------------------!
1463 SUBROUTINE chem_prognostic_equations_ij ( cs_scalar_p, cs_scalar, tcs_scalar_m, pr_init_cs,  &
1464                            i, j, i_omp_start, tn, ilsp, flux_s_cs, diss_s_cs, &
1465                            flux_l_cs, diss_l_cs )
1466    USE pegrid         
1467    USE advec_ws,        ONLY:  advec_s_ws 
1468    USE advec_s_pw_mod,  ONLY:  advec_s_pw
1469    USE advec_s_up_mod,  ONLY:  advec_s_up
1470    USE diffusion_s_mod, ONLY:  diffusion_s
1471    USE indices,         ONLY:  wall_flags_0
1472    USE surface_mod,     ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
1473                                surf_usm_v
1474
1475
1476    IMPLICIT NONE
1477
1478    REAL(wp), DIMENSION(:,:,:), POINTER   :: cs_scalar_p, cs_scalar, tcs_scalar_m
1479
1480    INTEGER(iwp),INTENT(IN) :: i, j, i_omp_start, tn, ilsp
1481    REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1)         :: flux_s_cs   !<
1482    REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1)         :: diss_s_cs   !<
1483    REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: flux_l_cs   !<
1484    REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: diss_l_cs   !<
1485    REAL(wp), DIMENSION(0:nz+1)                                 :: pr_init_cs  !<
1486
1487!-- local variables
1488
1489    INTEGER :: k
1490    !
1491    !--    Tendency-terms for chem spcs.
1492    tend(:,j,i) = 0.0_wp
1493!   
1494!-- Advection terms
1495    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1496       IF ( ws_scheme_sca )  THEN
1497          CALL advec_s_ws( i, j, cs_scalar, 'kc', flux_s_cs, diss_s_cs,                &
1498             flux_l_cs, diss_l_cs, i_omp_start, tn )
1499       ELSE
1500          CALL advec_s_pw( i, j, cs_scalar )
1501       ENDIF
1502    ELSE
1503         CALL advec_s_up( i, j, cs_scalar )
1504    ENDIF
1505
1506!
1507
1508!-- Diffusion terms (the last three arguments are zero)
1509
1510      CALL diffusion_s( i, j, cs_scalar,                                                 &
1511                        surf_def_h(0)%cssws(ilsp,:), surf_def_h(1)%cssws(ilsp,:),        &
1512                        surf_def_h(2)%cssws(ilsp,:),                                     &
1513                        surf_lsm_h%cssws(ilsp,:), surf_usm_h%cssws(ilsp,:),              &
1514                        surf_def_v(0)%cssws(ilsp,:), surf_def_v(1)%cssws(ilsp,:),        &
1515                        surf_def_v(2)%cssws(ilsp,:), surf_def_v(3)%cssws(ilsp,:),        &
1516                        surf_lsm_v(0)%cssws(ilsp,:), surf_lsm_v(1)%cssws(ilsp,:),        &
1517                        surf_lsm_v(2)%cssws(ilsp,:), surf_lsm_v(3)%cssws(ilsp,:),        &
1518                        surf_usm_v(0)%cssws(ilsp,:), surf_usm_v(1)%cssws(ilsp,:),        &
1519                        surf_usm_v(2)%cssws(ilsp,:), surf_usm_v(3)%cssws(ilsp,:) )
1520 
1521!   
1522!-- Prognostic equation for chem spcs
1523    DO k = nzb+1, nzt
1524       cs_scalar_p(k,j,i) = cs_scalar(k,j,i) + ( dt_3d  *                       &
1525                                               ( tsc(2) * tend(k,j,i) +         &
1526                                                 tsc(3) * tcs_scalar_m(k,j,i) ) & 
1527                                               - tsc(5) * rdf_sc(k)             &
1528                                                        * ( cs_scalar(k,j,i) - pr_init_cs(k) )    &   
1529                                               )                                                  &
1530                                                        * MERGE( 1.0_wp, 0.0_wp,                  &     
1531                                                                BTEST( wall_flags_0(k,j,i), 0 )   &             
1532                                                                 )       
1533
1534       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
1535    ENDDO
1536
1537!
1538!-- Calculate tendencies for the next Runge-Kutta step
1539    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1540       IF ( intermediate_timestep_count == 1 )  THEN
1541          DO  k = nzb+1, nzt
1542             tcs_scalar_m(k,j,i) = tend(k,j,i)
1543          ENDDO
1544       ELSEIF ( intermediate_timestep_count < &
1545          intermediate_timestep_count_max )  THEN
1546          DO  k = nzb+1, nzt
1547             tcs_scalar_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1548                5.3125_wp * tcs_scalar_m(k,j,i)
1549          ENDDO
1550       ENDIF
1551    ENDIF
1552
1553 END SUBROUTINE chem_prognostic_equations_ij
1554
1555
1556!------------------------------------------------------------------------------!
1557!
1558! Description:
1559! ------------
1560!> Subroutine calculating prognostic equations for chemical species
1561!> (vector-optimized).
1562!> Routine is called separately for each chemical species over a loop from
1563!> prognostic_equations.
1564!------------------------------------------------------------------------------!
1565 SUBROUTINE chem_prognostic_equations ( cs_scalar_p, cs_scalar, tcs_scalar_m,  &
1566                                        pr_init_cs, ilsp )
1567
1568    USE advec_s_pw_mod,                                                        &
1569        ONLY:  advec_s_pw
1570
1571    USE advec_s_up_mod,                                                        &
1572        ONLY:  advec_s_up
1573
1574    USE advec_ws,                                                              &
1575        ONLY:  advec_s_ws 
1576
1577    USE diffusion_s_mod,                                                       &
1578        ONLY:  diffusion_s
1579
1580    USE indices,                                                               &
1581        ONLY:  nxl, nxr, nyn, nys, wall_flags_0
1582
1583    USE pegrid
1584
1585    USE surface_mod,                                                           &
1586        ONLY:  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,     &
1587               surf_usm_v
1588
1589    IMPLICIT NONE
1590
1591    INTEGER ::  i   !< running index
1592    INTEGER ::  j   !< running index
1593    INTEGER ::  k   !< running index
1594
1595    INTEGER(iwp),INTENT(IN) ::  ilsp          !<
1596
1597    REAL(wp), DIMENSION(0:nz+1) ::  pr_init_cs   !<
1598
1599    REAL(wp), DIMENSION(:,:,:), POINTER ::  cs_scalar      !<
1600    REAL(wp), DIMENSION(:,:,:), POINTER ::  cs_scalar_p    !<
1601    REAL(wp), DIMENSION(:,:,:), POINTER ::  tcs_scalar_m   !<
1602
1603
1604!
1605!-- Tendency terms for chemical species
1606    tend = 0.0_wp
1607!   
1608!-- Advection terms
1609    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1610       IF ( ws_scheme_sca )  THEN
1611          CALL advec_s_ws( cs_scalar, 'kc' )
1612       ELSE
1613          CALL advec_s_pw( cs_scalar )
1614       ENDIF
1615    ELSE
1616         CALL advec_s_up( cs_scalar )
1617    ENDIF
1618!
1619!-- Diffusion terms  (the last three arguments are zero)
1620    CALL diffusion_s( cs_scalar,                                               &
1621                      surf_def_h(0)%cssws(ilsp,:),                             &
1622                      surf_def_h(1)%cssws(ilsp,:),                             &
1623                      surf_def_h(2)%cssws(ilsp,:),                             &
1624                      surf_lsm_h%cssws(ilsp,:),                                &
1625                      surf_usm_h%cssws(ilsp,:),                                &
1626                      surf_def_v(0)%cssws(ilsp,:),                             &
1627                      surf_def_v(1)%cssws(ilsp,:),                             &
1628                      surf_def_v(2)%cssws(ilsp,:),                             &
1629                      surf_def_v(3)%cssws(ilsp,:),                             &
1630                      surf_lsm_v(0)%cssws(ilsp,:),                             &
1631                      surf_lsm_v(1)%cssws(ilsp,:),                             &
1632                      surf_lsm_v(2)%cssws(ilsp,:),                             &
1633                      surf_lsm_v(3)%cssws(ilsp,:),                             &
1634                      surf_usm_v(0)%cssws(ilsp,:),                             &
1635                      surf_usm_v(1)%cssws(ilsp,:),                             &
1636                      surf_usm_v(2)%cssws(ilsp,:),                             &
1637                      surf_usm_v(3)%cssws(ilsp,:) )
1638!   
1639!-- Prognostic equation for chemical species
1640    DO  i = nxl, nxr
1641       DO  j = nys, nyn   
1642          DO  k = nzb+1, nzt
1643             cs_scalar_p(k,j,i) =   cs_scalar(k,j,i)                           &
1644                                  + ( dt_3d  *                                 &
1645                                      (   tsc(2) * tend(k,j,i)                 &
1646                                        + tsc(3) * tcs_scalar_m(k,j,i)         &
1647                                      )                                        & 
1648                                    - tsc(5) * rdf_sc(k)                       &
1649                                             * ( cs_scalar(k,j,i) - pr_init_cs(k) )    &   
1650                                    )                                          &
1651                                    * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )       
1652
1653             IF ( cs_scalar_p(k,j,i) < 0.0_wp )  cs_scalar_p(k,j,i) = 0.1_wp * cs_scalar(k,j,i)
1654          ENDDO
1655       ENDDO
1656    ENDDO
1657!
1658!-- Calculate tendencies for the next Runge-Kutta step
1659    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1660       IF ( intermediate_timestep_count == 1 )  THEN
1661          DO  i = nxl, nxr
1662             DO  j = nys, nyn   
1663                DO  k = nzb+1, nzt
1664                   tcs_scalar_m(k,j,i) = tend(k,j,i)
1665                ENDDO
1666             ENDDO
1667          ENDDO
1668       ELSEIF ( intermediate_timestep_count < &
1669          intermediate_timestep_count_max )  THEN
1670          DO  i = nxl, nxr
1671             DO  j = nys, nyn
1672                DO  k = nzb+1, nzt
1673                   tcs_scalar_m(k,j,i) = - 9.5625_wp * tend(k,j,i)             &
1674                                         + 5.3125_wp * tcs_scalar_m(k,j,i)
1675                ENDDO
1676             ENDDO
1677          ENDDO
1678       ENDIF
1679    ENDIF
1680
1681 END SUBROUTINE chem_prognostic_equations
1682
1683
1684!------------------------------------------------------------------------------!
1685!
1686! Description:
1687! ------------
1688!> Subroutine defining header output for chemistry model
1689!------------------------------------------------------------------------------!
1690 SUBROUTINE chem_header ( io )
1691       
1692    IMPLICIT NONE
1693 
1694       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
1695
1696!      print*,'the header subroutine is still not operational'     
1697!!
1698!!--    Write chemistry model header
1699!       WRITE( io, 3 )
1700!
1701!       IF ( radiation_scheme == "constant" )  THEN
1702!          WRITE( io, 4 ) net_radiation
1703!       ELSEIF ( radiation_scheme == "clear-sky" )  THEN
1704!          WRITE( io, 5 )
1705!       ELSEIF ( radiation_scheme == "rrtmg" )  THEN
1706!          WRITE( io, 6 )
1707!          IF ( .NOT. lw_radiation )  WRITE( io, 10 )
1708!          IF ( .NOT. sw_radiation )  WRITE( io, 11 )
1709!       ENDIF
1710!
1711!       IF ( albedo_type == 0 )  THEN
1712!          WRITE( io, 7 ) albedo
1713!       ELSE
1714!          WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
1715!       ENDIF
1716!       IF ( constant_albedo )  THEN
1717!          WRITE( io, 9 )
1718!       ENDIF
1719!     
1720!       IF ( radiation .AND. radiation_scheme /= 'constant' )  THEN
1721!          WRITE ( io, 1 )  lambda
1722!          WRITE ( io, 2 )  day_init, time_utc_init
1723!       ENDIF
1724!
1725!       WRITE( io, 12 ) dt_radiation
1726!
1727! 1 FORMAT ('    Geograph. longitude            :   lambda = ',F4.1,' degr')
1728! 2 FORMAT ('    Day of the year at model start :   day_init = ',I3   &
1729!            /'    UTC time at model start        :   time_utc_init = ',F7.1' s')
1730! 3 FORMAT (//' Radiation model information:'/                                  &
1731!              ' ----------------------------'/)
1732! 4 FORMAT ('    --> Using constant net radiation: net_radiation = ', F6.2,        &
1733!           // 'W/m**2')
1734! 5 FORMAT ('    --> Simple radiation scheme for clear sky is used (no clouds,',   &
1735!                   ' default)')
1736! 6 FORMAT ('    --> RRTMG scheme is used')
1737! 7 FORMAT (/'    User-specific surface albedo: albedo =', F6.3)
1738! 8 FORMAT (/'    Albedo is set for land surface type: ', A)
1739! 9 FORMAT (/'    --> Albedo is fixed during the run')
1740!10 FORMAT (/'    --> Longwave radiation is disabled')
1741!11 FORMAT (/'    --> Shortwave radiation is disabled.')
1742!12 FORMAT  ('    Timestep: dt_radiation = ', F6.2, '  s')
1743!
1744!
1745 END SUBROUTINE chem_header
1746
1747!------------------------------------------------------------------------------
1748! Description:
1749! ------------
1750!> Subroutine reading restart data for chemistry model input parameters
1751!  (FK: To make restarts work, I had to comment this routine. We actually
1752!       don't need it, since the namelist parameters are always read in,
1753!       also in case of a restart run)
1754!------------------------------------------------------------------------------
1755!  SUBROUTINE chem_rrd_global
1756!
1757!    USE chem_modules
1758!
1759!    USE control_parameters,                                                   &
1760!        ONLY: length, message_string, restart_string
1761!
1762!     
1763!     IMPLICIT NONE
1764!     
1765!     
1766!     
1767!     DO
1768!
1769!        SELECT CASE ( restart_string(1:length) )
1770!       
1771!           CASE ( 'bc_cs_b' )
1772!              READ ( 13 )  bc_cs_b
1773!
1774!           CASE DEFAULT
1775!
1776!              EXIT
1777!           
1778!        END SELECT
1779!       
1780!!
1781!!--     Read next string and its length
1782!        READ ( 13 )  length
1783!        READ ( 13 )  restart_string(1:length)
1784!       
1785!     ENDDO
1786!     
1787!  END SUBROUTINE chem_rrd_global
1788
1789
1790!------------------------------------------------------------------------------!
1791!
1792! Description:
1793! ------------
1794!> Subroutine writing restart data for chemistry model input parameters
1795!  (FK: To make restarts work, I had to comment this routine. We actually
1796!       don't need it, since the namelist parameters are always read in,
1797!       also in case of a restart run)
1798!------------------------------------------------------------------------------!
1799!  SUBROUTINE chem_wrd_global
1800!
1801!     USE chem_modules
1802!     
1803!     USE kinds
1804!
1805!
1806!     IMPLICIT NONE
1807!
1808!     INTEGER(iwp) ::  lsp  !< running index for chem spcs
1809!
1810! !
1811! !-- Writing out input parameters that are not part of chemistry_parameters
1812! !-- namelist (namelist parameters are anyway read in again in case of restart)
1813!     DO lsp = 1, nvar
1814!        CALL wrd_write_string( 'conc_pr_init_'//chem_species(lsp)%name )
1815!        WRITE ( 14 )  chem_species(lsp)%conc_pr_init
1816!     ENDDO
1817!
1818!
1819!  END SUBROUTINE chem_wrd_global
1820
1821
1822!------------------------------------------------------------------------------!
1823!
1824! Description:
1825! ------------
1826!> Subroutine for emission
1827!------------------------------------------------------------------------------!
1828 SUBROUTINE chem_emissions
1829
1830    USE chem_modules
1831   
1832    USE netcdf_data_input_mod,                                                 &
1833        ONLY:  street_type_f
1834   
1835    USE surface_mod,                                                           &
1836        ONLY:  surf_lsm_h
1837
1838
1839    IMPLICIT NONE
1840
1841    INTEGER(iwp) ::  i    !< running index for grid in x-direction
1842    INTEGER(iwp) ::  j    !< running index for grid in y-direction
1843    INTEGER(iwp) ::  m    !< running index for horizontal surfaces
1844    INTEGER(iwp) ::  lsp  !< running index for chem spcs
1845
1846!
1847!-- Comment??? (todo)
1848    IF ( street_type_f%from_file )  THEN
1849!
1850!--    Streets are lsm surfaces, hence, no usm surface treatment required
1851       DO  m = 1, surf_lsm_h%ns
1852          i = surf_lsm_h%i(m)
1853          j = surf_lsm_h%j(m)
1854         
1855          IF ( street_type_f%var(j,i) >= main_street_id  .AND.                 &
1856               street_type_f%var(j,i) < max_street_id )  THEN
1857             DO  lsp = 1, nvar
1858                surf_lsm_h%cssws(lsp,m) = emiss_factor_main * surface_csflux(lsp)
1859             ENDDO
1860          ELSEIF ( street_type_f%var(j,i) >= side_street_id  .AND.             &
1861                   street_type_f%var(j,i) < main_street_id )  THEN
1862             DO  lsp = 1, nvar
1863                surf_lsm_h%cssws(lsp,m) = emiss_factor_side * surface_csflux(lsp)
1864             ENDDO
1865          ELSE
1866             surf_lsm_h%cssws(:,m) = 0.0_wp
1867          ENDIF
1868       ENDDO
1869       
1870    ENDIF
1871
1872 END SUBROUTINE chem_emissions
1873
1874
1875 END MODULE chemistry_model_mod
1876
Note: See TracBrowser for help on using the repository browser.