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

Last change on this file since 2768 was 2768, checked in by kanani, 5 years ago

Added parameter check, reduced line length, some formatting

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