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

Last change on this file since 3226 was 3215, checked in by suehring, 6 years ago

changes for commit 3209 documented

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