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

Last change on this file since 2766 was 2766, checked in by kanani, 4 years ago

Removal of chem directive, plus minor changes

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