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

Last change on this file since 2742 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

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